From 772663a03c4eaabd247ecc9d445ebdbb785868f5 Mon Sep 17 00:00:00 2001 From: tangchao0503 <735056338@qq.com> Date: Mon, 6 Jan 2025 10:18:08 +0800 Subject: [PATCH] =?UTF-8?q?=E5=AE=9E=E7=8E=B0=E6=B0=B4=E4=BD=93=E5=90=AB?= =?UTF-8?q?=E9=87=8F=E5=8F=8D=E6=BC=94=E7=A8=8B=E5=BA=8F=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- deglint.py | 58 ++++ deglint_oxygen_absorption_valley.py | 62 ++++ deglint_regression_slope.py | 212 +++++++++++++ deglint_subtract_nir.py | 35 +++ extract_water_area.py | 142 +++++++++ find_severe_glint_area.py | 118 +++++++ get_coor_base_interval.py | 81 +++++ get_spectral_in_coor.py | 54 ++++ interpolate.py | 77 +++++ model_correction.py | 461 ++++++++++++++++++++++++++++ readme.txt | 214 +++++++++++++ requirements.txt | Bin 0 -> 302 bytes retrieval.py | 169 ++++++++++ type_define.py | 22 ++ util.py | 174 +++++++++++ 算法步骤.txt | 24 ++ 16 files changed, 1903 insertions(+) create mode 100644 deglint.py create mode 100644 deglint_oxygen_absorption_valley.py create mode 100644 deglint_regression_slope.py create mode 100644 deglint_subtract_nir.py create mode 100644 extract_water_area.py create mode 100644 find_severe_glint_area.py create mode 100644 get_coor_base_interval.py create mode 100644 get_spectral_in_coor.py create mode 100644 interpolate.py create mode 100644 model_correction.py create mode 100644 readme.txt create mode 100644 requirements.txt create mode 100644 retrieval.py create mode 100644 type_define.py create mode 100644 util.py create mode 100644 算法步骤.txt diff --git a/deglint.py b/deglint.py new file mode 100644 index 0000000..7006897 --- /dev/null +++ b/deglint.py @@ -0,0 +1,58 @@ +import argparse + +from deglint_subtract_nir import deglint_subtract_nir +from deglint_oxygen_absorption_valley import deglint_oxygen_absorption_valley +from deglint_regression_slope import deglint_regression_slope + + +def main(): + parser = argparse.ArgumentParser(description="此程序用于去除水域数据的耀斑，暂时支持3种算法。") + + # parser.add_argument("--global_arg", type=str, help="A global argument for all modes", required=True) + + # 创建子命令解析器 + subparsers = parser.add_subparsers(dest="algorithm", required=True, help="Choose a mode") + + subnir = subparsers.add_parser("subnir", help="Mode 1 description") + subnir.add_argument('-i1', '--input', type=str, required=True, help='输入影像文件的路径') + subnir.add_argument('-i2', '--input_water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + subnir.add_argument('-o', '--output', type=str, required=True, help='输出文件的路径') + subnir.add_argument('-s', '--start_wave', type=float, required=True, help='nir开始波段') + subnir.add_argument('-e', '--end_wave', type=float, required=True, help='nir结束波段') + subnir.set_defaults(func=deglint_subtract_nir) + + regression_slope = subparsers.add_parser("regression_slope", help="Mode 2 description") + regression_slope.add_argument('-i1', '--input', type=str, required=True, help='输入影像文件的路径') + regression_slope.add_argument('-i2', '--input_water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + regression_slope.add_argument('-o', '--output', type=str, required=True, help='输出文件的路径') + regression_slope.add_argument('-s', '--start_wave', type=float, required=True, help='nir开始波段') + regression_slope.add_argument('-e', '--end_wave', type=float, required=True, help='nir结束波段') + regression_slope.add_argument('-r', '--roi', type=str, required=True, help='输入roi文件的路径') + regression_slope.set_defaults(func=deglint_regression_slope) + + oxygen_absorption = subparsers.add_parser("oxygen_absorption", help="Mode 3 description") + oxygen_absorption.add_argument('-i1', '--input', type=str, required=True, help='输入影像文件的路径') + oxygen_absorption.add_argument('-i2', '--input_water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + oxygen_absorption.add_argument('-o', '--output', type=str, required=True, help='输出文件的路径') + oxygen_absorption.add_argument('-l', '--left_shoulder_wave', type=float, required=True, help='氧气吸收谷左肩波长') + oxygen_absorption.add_argument('-v', '--valley_wave', type=float, required=True, help='氧气吸收谷底波长') + oxygen_absorption.add_argument('-r', '--right_shoulder_wave', type=float, required=True, help='氧气吸收谷右肩波长') + + oxygen_absorption.add_argument('-sw', '--shoulder_window', type=float, required=True, help='氧气吸收谷左右肩平均窗口，半径') + oxygen_absorption.add_argument('-vw', '--valley_window', type=float, required=True, help='氧气吸收谷底平均窗口，半径') + oxygen_absorption.set_defaults(func=deglint_oxygen_absorption_valley) + + # 解析参数 + args = parser.parse_args() + if args.algorithm == "subnir": # 处理文件D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510时，730-750效果是最好的 + args.func(args.input, args.input_water_mask, args.output, args.start_wave, args.end_wave) + elif args.algorithm == "regression_slope": + args.func(args.input, args.input_water_mask, args.output, args.start_wave, args.end_wave, args.roi) + elif args.algorithm == "oxygen_absorption": + args.func(args.input, args.input_water_mask, args.output, args.left_shoulder_wave, args.valley_wave, args.right_shoulder_wave, + args.shoulder_window, args.valley_window) + + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + main() diff --git a/deglint_oxygen_absorption_valley.py b/deglint_oxygen_absorption_valley.py new file mode 100644 index 0000000..58b1146 --- /dev/null +++ b/deglint_oxygen_absorption_valley.py @@ -0,0 +1,62 @@ +import numpy as np + +from util import * +from osgeo import gdal + + +def calculate_glint_spatial_distribution(imgpath, left_shoulder_wave, valley_wave, right_shoulder_wave, shoulder_window=1, valley_window=1): + left_shoulder_band = average_bands(left_shoulder_wave - shoulder_window, left_shoulder_wave + shoulder_window, imgpath) + valley_band = average_bands(valley_wave - valley_window, valley_wave + valley_window, imgpath) + right_shoulder_band = average_bands(right_shoulder_wave - shoulder_window, right_shoulder_wave + shoulder_window, imgpath) + + glint_spatial_distribution = (left_shoulder_band + right_shoulder_band) / 2 - valley_band + glint_spatial_distribution = glint_spatial_distribution / glint_spatial_distribution.max() + + glint_spatial_distribution[np.where(glint_spatial_distribution < 0)] = 0 + + write_bands(imgpath, r"D:\PycharmProjects\sun_glint\test_data\glint_spatial_distribution", glint_spatial_distribution) + + return glint_spatial_distribution + + +@timeit +def deglint_oxygen_absorption_valley(imgpath, water_mask, out_imgpath, left_shoulder_wave, valley_wave, right_shoulder_wave, shoulder_window=1, valley_window=1): + glint_spatial_distribution = calculate_glint_spatial_distribution(imgpath, left_shoulder_wave, valley_wave, + right_shoulder_wave, shoulder_window, + valley_window) + + dataset = gdal.Open(imgpath) + im_width = dataset.RasterXSize + im_height = dataset.RasterYSize + num_bands = dataset.RasterCount + geotransform = dataset.GetGeoTransform() + im_proj = dataset.GetProjection() + + # 通过data ignore value获取影像有效区域 + + format = "ENVI" + driver = gdal.GetDriverByName(format) + dst_ds = driver.Create(out_imgpath, im_width, im_height, num_bands, gdal.GDT_Int16, + options=["INTERLEAVE=BSQ"]) + dst_ds.SetGeoTransform(geotransform) + dst_ds.SetProjection(im_proj) + + dataset_water_mask = gdal.Open(water_mask) + data_water_mask = dataset_water_mask.GetRasterBand(1).ReadAsArray() + del dataset_water_mask + for i in range(num_bands): + band = dataset.GetRasterBand(i + 1).ReadAsArray() + + # 找波段最大最下值时，需要掩膜掉没有影像的区域（一般为0值的区域） + glint = glint_spatial_distribution * (band.max() - min(band[np.where(band > 0)])) * 1 + glint[np.where(data_water_mask == 0)] = 0 + band_deglint = band - glint + + # 确保值大于等于0 + band_deglint[np.where(band_deglint < 0)] = 0 + + dst_ds.GetRasterBand(i + 1).WriteArray(band_deglint) + + del dataset, dst_ds + + write_fields_to_hdrfile(get_hdr_file_path(imgpath), get_hdr_file_path(out_imgpath)) diff --git a/deglint_regression_slope.py b/deglint_regression_slope.py new file mode 100644 index 0000000..f46ccdc --- /dev/null +++ b/deglint_regression_slope.py @@ -0,0 +1,212 @@ +import numpy as np + +from util import * +from osgeo import gdal +import json +from shapely.geometry import Point, Polygon +import xml.etree.ElementTree as ET + + +def parse_xml(roifilepath_xml): + # 读取 XML 文件 + tree = ET.parse(roifilepath_xml) # 替换为您的 XML 文件路径 + root = tree.getroot() + + # 打印根标签 + # print("Root tag:", root.tag) + + # 遍历 XML 文件的所有子元素 + # for child in root: + # print("Tag:", child.tag, "| Attributes:", child.attrib) + + # 查找特定标签 + polygons = [] + for region in root.findall(".//Region"): + # print("Region name:", region.get("name")) + for polygon in region.findall(".//Polygon"): + coordinates_str = polygon.find(".//Coordinates").text.strip() + coordinates = list(map(float, coordinates_str.split())) + points = [(coordinates[i], coordinates[i + 1]) for i in range(0, len(coordinates), 2)] + polygon = Polygon(points) + polygons.append(polygon) + # print("Polygon coordinates:", coordinates) + + return CoorType.depend_on_image, polygons + + +def create_json(filepath_json): + """ + 用户创建sample json文件（parse_json函数需要能够解析），创建逆时针闭合不规则多边形 + :param filepath_json: + :return: + """ + polygons = { + "Available_Coordinates_type": ["depend_on_image", "pixel"], + "Coordinates_type": "pixel", + "polygons": [ + { + "id": 1, + "coordinates": [ + {"x": 120, "y": 300}, + {"x": 250, "y": 220}, + {"x": 340, "y": 400}, + {"x": 200, "y": 500}, + {"x": 120, "y": 300} # 闭合 + ] + }, + { + "id": 2, + "coordinates": [ + {"x": 400, "y": 100}, + {"x": 550, "y": 180}, + {"x": 480, "y": 350}, + {"x": 370, "y": 280}, + {"x": 400, "y": 100} # 闭合 + ] + } + ] + } + + # 保存为 JSON 文件 + with open(filepath_json, "w") as json_file: + json.dump(polygons, json_file, indent=4) + + +def parse_json(roifilepath_json): + with open(roifilepath_json, "r") as json_file: + loaded_data = json.load(json_file) + + coor_type_tmp = loaded_data["Coordinates_type"] + if coor_type_tmp == "pixel": + coor_type = CoorType.pixel + elif coor_type_tmp == "depend_on_image": + coor_type = CoorType.depend_on_image + else: + print("Error: unknown coordinates type!!!") + return + + # 构建多边形 + polygons = [] + for poly in loaded_data["polygons"]: + # 提取多边形坐标 + coordinates = [(point["x"], point["y"]) for point in poly["coordinates"]] + # 使用 Shapely 构造多边形 + polygons.append(Polygon(coordinates)) + + return coor_type, polygons + + +def parse_roifile(roifilepath): + extension = os.path.splitext(roifilepath)[-1].removeprefix(".") + + if extension == "json": + coor_type, polygons = parse_json(roifilepath) + elif extension == "xml": + coor_type, polygons = parse_xml(roifilepath) + else: + print("Error: unknown roi file type!!!") + return + + return coor_type, polygons + + +def get_coors_in_polygons(polygons): + """ + :param polygons: list类型，坐标类型必须为像素坐标 + :return: + """ + x = [] + y = [] + for i in range(len(polygons)): + polygon = polygons[i] + bounds = polygon.bounds + + for x_tmp in range(int(bounds[0]), int(bounds[2])): + for y_tmp in range(int(bounds[1]), int(bounds[3])): + if polygon.contains(Point(x_tmp, y_tmp)): + x.append(x_tmp) + y.append(y_tmp) + + return x, y + + +@timeit +def get_pixel_coors_in_polygons(roifilepath, dataset): + coor_type, polygons = parse_roifile(roifilepath) + + polygons_pixel_coors = [] + + # 如果polygons坐标不是像素坐标，则转换为像素坐标 + if coor_type == CoorType.depend_on_image: + geotransform_input = dataset.GetGeoTransform() + inv_geotransform_input = gdal.InvGeoTransform(geotransform_input) + + for polygon in polygons: + shifted_coords = [gdal.ApplyGeoTransform(inv_geotransform_input, x, y) for x, y in polygon.exterior.coords] + shifted_polygon = Polygon(shifted_coords) + polygons_pixel_coors.append(shifted_polygon) + elif coor_type == CoorType.pixel: + polygons_pixel_coors = polygons + + return get_coors_in_polygons(polygons_pixel_coors) + + +def get_valid_pixel_value(band, x, y): + value = band[y, x] + value = np.sort(value) + + return value + + +@timeit +def deglint_regression_slope(imgpath, water_mask, out_imgpath, start_nir_wave, end_nir_wave, filepath_json): + """ + :param imgpath: + :param out_imgpath: + :param start_nir_wave: + :param end_nir_wave: + :param args: + :param kwargs: + :return: + """ + nir = average_bands(start_nir_wave, end_nir_wave, imgpath) + + dataset = gdal.Open(imgpath) + im_width = dataset.RasterXSize + im_height = dataset.RasterYSize + num_bands = dataset.RasterCount + geotransform = dataset.GetGeoTransform() + im_proj = dataset.GetProjection() + + valid_pixel_coor1_x, valid_pixel_coor1_y = get_pixel_coors_in_polygons(filepath_json, dataset) + x = get_valid_pixel_value(nir, valid_pixel_coor1_x, valid_pixel_coor1_y) + + format = "ENVI" + driver = gdal.GetDriverByName(format) + dst_ds = driver.Create(out_imgpath, im_width, im_height, num_bands, gdal.GDT_Int16, + options=["INTERLEAVE=BSQ"]) + dst_ds.SetGeoTransform(geotransform) + dst_ds.SetProjection(im_proj) + + dataset_water_mask = gdal.Open(water_mask) + data_water_mask = dataset_water_mask.GetRasterBand(1).ReadAsArray() + del dataset_water_mask + for i in range(num_bands): + band = dataset.GetRasterBand(i + 1).ReadAsArray() + y = get_valid_pixel_value(band, valid_pixel_coor1_x, valid_pixel_coor1_y) + coefficients = np.polyfit(x, y, 1) + + glint = (nir - min(x)) * coefficients[0] + glint[np.where(glint < 0)] = 0 + glint[np.where(data_water_mask == 0)] = 0 + + band_deglint = band - glint + + # 确保值大于等于0 + band_deglint[np.where(band_deglint < 0)] = 0 + + dst_ds.GetRasterBand(i + 1).WriteArray(band_deglint) + + del dataset, dst_ds + + write_fields_to_hdrfile(get_hdr_file_path(imgpath), get_hdr_file_path(out_imgpath)) diff --git a/deglint_subtract_nir.py b/deglint_subtract_nir.py new file mode 100644 index 0000000..cbd09c4 --- /dev/null +++ b/deglint_subtract_nir.py @@ -0,0 +1,35 @@ +import numpy as np +from util import * +from osgeo import gdal + + +@timeit +def deglint_subtract_nir(imgpath, water_mask, out_imgpath, start_wave, end_wave): + glint = average_bands_in_mask(start_wave, end_wave, imgpath, water_mask) + + dataset = gdal.Open(imgpath) + im_width = dataset.RasterXSize + im_height = dataset.RasterYSize + num_bands = dataset.RasterCount + geotransform = dataset.GetGeoTransform() + im_proj = dataset.GetProjection() + + format = "ENVI" + driver = gdal.GetDriverByName(format) + dst_ds = driver.Create(out_imgpath, im_width, im_height, num_bands, gdal.GDT_Int16, + options=["INTERLEAVE=BSQ"]) + dst_ds.SetGeoTransform(geotransform) + dst_ds.SetProjection(im_proj) + + for i in range(num_bands): + band = dataset.GetRasterBand(i + 1).ReadAsArray() + band_deglint = band - glint + + # 确保值大于等于0 + band_deglint[np.where(band_deglint < 0)] = 0 + + dst_ds.GetRasterBand(i + 1).WriteArray(band_deglint) + + del dataset, dst_ds + + write_fields_to_hdrfile(get_hdr_file_path(imgpath), get_hdr_file_path(out_imgpath)) diff --git a/extract_water_area.py b/extract_water_area.py new file mode 100644 index 0000000..698163a --- /dev/null +++ b/extract_water_area.py @@ -0,0 +1,142 @@ +from util import * +from osgeo import gdal +import argparse + + +def xml2shp(): + pass + + +def rasterize_envi_xml(shp_filepath): + pass + + +@timeit +def rasterize_shp(shp_filepath, raster_fn_out, img_path, NoData_value=0): + dataset = gdal.Open(img_path) + im_width = dataset.RasterXSize + im_height = dataset.RasterYSize + geotransform = dataset.GetGeoTransform() + imgdata_in = dataset.GetRasterBand(1).ReadAsArray() + del dataset + + # Open the data source and read in the extent + source_ds = gdal.OpenEx(shp_filepath) + # about 25 metres(ish) use 0.001 if you want roughly 100m + pixel_size = geotransform[1] + raster_fn_out_tmp = append2filename(raster_fn_out, "_tmp_delete") + gdal.Rasterize(raster_fn_out_tmp, source_ds, format='envi', outputType=gdal.GDT_Byte, + noData=NoData_value, initValues=NoData_value, xRes=pixel_size, yRes=-pixel_size, allTouched=True, + burnValues=1) + + dataset_tmp = gdal.Open(raster_fn_out_tmp) + geotransform_tmp = dataset_tmp.GetGeoTransform() + inv_geotransform_tmp = gdal.InvGeoTransform(geotransform_tmp) + data_tmp = dataset_tmp.GetRasterBand(1).ReadAsArray() + del dataset_tmp + + # 创建和输入影像相同行列号、相同分辨率的水域掩膜，方便后续使用 + water_mask = np.zeros((im_height, im_width)) + for row in range(im_height): + for column in range(im_width): + coor = gdal.ApplyGeoTransform(geotransform, column, row) + + coor_pixel = gdal.ApplyGeoTransform(inv_geotransform_tmp, coor[0], coor[1]) + coor_pixel = [int(num) for num in coor_pixel] + + if coor_pixel[0] < 0 or coor_pixel[0] >= data_tmp.shape[1]: + continue + if coor_pixel[1] < 0 or coor_pixel[1] >= data_tmp.shape[0]: + continue + + if imgdata_in[row, column] == 0: # 当shp区域比影像区域大时，略过 + continue + + water_mask[row, column] = data_tmp[coor_pixel[1], coor_pixel[0]] + + write_bands(img_path, raster_fn_out, water_mask) + + os.remove(raster_fn_out_tmp) + + +def calculate_NDWI(green_bandnumber, nir_bandnumber, filename): + dataset = gdal.Open(filename) # 打开文件 + num_bands = dataset.RasterCount # 栅格矩阵的波段数 + im_geotrans = dataset.GetGeoTransform() # 仿射矩阵 + im_proj = dataset.GetProjection() # 地图投影信息 + + tmp = dataset.GetRasterBand(green_bandnumber + 1) # 波段计数从1开始 + band_green = tmp.ReadAsArray().astype(np.int16) + tmp = dataset.GetRasterBand(nir_bandnumber + 1) # 波段计数从1开始 + band_nir = tmp.ReadAsArray().astype(np.int16) + + ndwi = (band_green - band_nir) / (band_green + band_nir) + + del dataset + + return ndwi + + +def extract_water(ndwi, threshold=0.3, data_ignore_value=0): + water_region = np.where(ndwi > threshold, 1, data_ignore_value) + + return water_region + + +def ndwi(file_path, ndwi_threshold=0.4, output_path=None, data_ignore_value=0): + if output_path is None: + output_path = append2filename(file_path, "_waterarea") + + dataset_in = gdal.Open(file_path) + im_width_in = dataset_in.RasterXSize # 栅格矩阵的列数 + im_height_in = dataset_in.RasterYSize # 栅格矩阵的行数 + num_bands_in = dataset_in.RasterCount # 栅格矩阵的波段数 + geotrans_in = dataset_in.GetGeoTransform() # 仿射矩阵 + proj_in = dataset_in.GetProjection() # 地图投影信息 + del dataset_in + + green_wave = 552.19 + nir_wave = 809.2890 + green_band_number = find_band_number(green_wave, file_path) + nir_band_number = find_band_number(nir_wave, file_path) + + ndwi = calculate_NDWI(green_band_number, nir_band_number, file_path) + + water_binary = extract_water(ndwi, threshold=ndwi_threshold) # 0.4 + + write_bands(file_path, output_path, water_binary) + + return output_path + + +def main(): + parser = argparse.ArgumentParser(description="此程序用于提取水域区域，输出的水域栅格和输入的影像具有相同的行列数。") + + # parser.add_argument("--global_arg", type=str, help="A global argument for all modes", required=True) + + # 创建子命令解析器 + subparsers = parser.add_subparsers(dest="algorithm", required=True, help="Choose a mode") + + rasterize_shp_ = subparsers.add_parser("rasterize_shp", help="Mode 1 description") + rasterize_shp_.add_argument('-i1', '--img_path', type=str, required=True, help='输入影像文件的路径') + rasterize_shp_.add_argument('-i2', '--shp_path', type=str, required=True, help='输入shp文件的路径') + rasterize_shp_.add_argument('-o', '--water_mask_outpath', required=True, type=str, help='输出水体掩膜文件的路径') + rasterize_shp_.set_defaults(func=rasterize_shp) + + ndwi_ = subparsers.add_parser("ndwi", help="Mode 2 description") + ndwi_.add_argument('-i1', '--img_path', type=str, required=True, help='输入影像文件的路径') + ndwi_.add_argument('-i2', '--ndwi_threshold', type=float, required=True, help='输入ndwi水体阈值，大于此值的为水域') + ndwi_.add_argument('-o', '--water_mask_outpath', required=True, type=str, help='输出水体掩膜文件的路径') + ndwi_.set_defaults(func=ndwi) + + # 解析参数 + args = parser.parse_args() + if args.algorithm == "rasterize_shp": + args.func(args.shp_path, args.water_mask_outpath, args.img_path) + elif args.algorithm == "ndwi": + args.func(args.img_path, args.ndwi_threshold, args.water_mask_outpath) + + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + main() diff --git a/find_severe_glint_area.py b/find_severe_glint_area.py new file mode 100644 index 0000000..015e1ce --- /dev/null +++ b/find_severe_glint_area.py @@ -0,0 +1,118 @@ +from util import * +from osgeo import gdal +import argparse + + +@timeit +def otsu(img, max_value, data_water_mask, ignore_value=0, foreground=1, background=0): + height = img.shape[0] + width = img.shape[1] + + hist = np.zeros([max_value], np.float32) + + # 计算直方图 + invalid_counter = 0 + for i in range(height): + for j in range(width): + if img[i, j] == ignore_value or img[i, j] < 0 or data_water_mask[i, j] == 0: + invalid_counter = invalid_counter + 1 + continue + + hist[img[i, j]] += 1 + + hist /= (height * width - invalid_counter) + + threshold = 0 + deltaMax = 0 + # 遍历像素值，计算最大类间方差 + for i in range(max_value): + wA = 0 + wB = 0 + uAtmp = 0 + uBtmp = 0 + uA = 0 + uB = 0 + u = 0 + for j in range(max_value): + if j <= i: + wA += hist[j] + uAtmp += j * hist[j] + else: + wB += hist[j] + uBtmp += j * hist[j] + if wA == 0: + wA = 1e-10 + if wB == 0: + wB = 1e-10 + uA = uAtmp / wA + uB = uBtmp / wB + u = uAtmp + uBtmp + + # 计算类间方差 + deltaTmp = wA * ((uA - u)**2) + wB * ((uB - u)**2) + # 找出最大类间方差以及阈值 + if deltaTmp > deltaMax: + deltaMax = deltaTmp + threshold = i + + # 二值化 + det_img = img.copy() + det_img[img > threshold] = foreground + det_img[img <= threshold] = background + det_img[np.where(data_water_mask == 0)] = background + return det_img + + +def find_overexposure_area(img_path, threhold=4095): + # 第一步通过某个像素的光谱找到信号最强的波段 + + # 根据上步所得的波段号检测过曝区域 + pass + + +@timeit +def find_severe_glint_area(img_path, water_mask, glint_wave=750, output_path=None): + if output_path is None: + output_path = append2filename(img_path, "_severe_glint_area") + + glint_band_number = find_band_number(glint_wave, img_path) + + dataset = gdal.Open(img_path) + num_bands = dataset.RasterCount + im_geotrans = dataset.GetGeoTransform() + im_proj = dataset.GetProjection() + + tmp = dataset.GetRasterBand(glint_band_number + 1) # 波段计数从1开始 + band_flare = tmp.ReadAsArray().astype(np.int16) + band_flare_stretch = (band_flare / band_flare.max() * 255).astype(np.int32) + + dataset_water_mask = gdal.Open(water_mask) + data_water_mask = dataset_water_mask.GetRasterBand(1).ReadAsArray() + del dataset_water_mask + + # 不加1报错：IndexError: index 73 is out of bounds for axis 0 with size 73 + flare_binary = otsu(band_flare_stretch, band_flare_stretch.max() + 1, data_water_mask) + + write_bands(img_path, output_path, flare_binary) + + del dataset + + return output_path + + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + img_path = r"D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_bsq" + + parser = argparse.ArgumentParser(description="此程序通过大律法分割图像，提取耀斑最严重的区域，输出的栅格和输入的影像具有相同的行列数。") + + parser.add_argument('-i1', '--input', type=str, required=True, help='输入影像文件的路径') + parser.add_argument('-i2', '--input_water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + parser.add_argument('-gw', '--glint_wave', type=float, default=750.0, help='用于提取耀斑严重区域的波段.') + parser.add_argument('-o', '--output', type=str, help='输出文件的路径') + + parser.add_argument('-v', '--verbose', action='store_true', help='启用详细模式') + + args = parser.parse_args() + + find_severe_glint_area(args.input, args.input_water_mask, args.glint_wave, args.output) diff --git a/get_coor_base_interval.py b/get_coor_base_interval.py new file mode 100644 index 0000000..1cf3888 --- /dev/null +++ b/get_coor_base_interval.py @@ -0,0 +1,81 @@ +from util import * +import math +import argparse + + +def get_coor_base_interval(water_mask, severe_glint, output_csvpath, interval=100): + dataset_water_mask = gdal.Open(water_mask) + # geotransform_water_mask = dataset_water_mask.GetGeoTransform() # 仿射矩阵 + data_water_mask = dataset_water_mask.GetRasterBand(1).ReadAsArray() + + dataset_severe_glint = gdal.Open(severe_glint) + data_severe_glint = dataset_severe_glint.GetRasterBand(1).ReadAsArray() + im_width = dataset_severe_glint.RasterXSize # 栅格矩阵的列数 + im_height = dataset_severe_glint.RasterYSize # 栅格矩阵的行数 + geotransform_input = dataset_severe_glint.GetGeoTransform() # 仿射矩阵 + inv_geotransform_input = gdal.InvGeoTransform(geotransform_input) + + x_max = geotransform_input[0] + im_width * geotransform_input[1] + im_height * geotransform_input[2] + y_min = geotransform_input[3] + im_width * geotransform_input[4] + im_height * geotransform_input[5] + + x_range = [geotransform_input[0], x_max] + y_range = [y_min, geotransform_input[3]] + + spatial_resolution = geotransform_input[1] * interval + dx = math.ceil((x_range[1] - x_range[0]) / spatial_resolution) + dy = math.ceil((y_range[1] - y_range[0]) / spatial_resolution) + geotransform_out = (x_range[0], spatial_resolution, 0, y_range[1], 0, -spatial_resolution) + + # 确定像元的投影范围 + valid_area_ = (data_water_mask > 0) & (~ (data_severe_glint > 0)) + valid_area = np.logical_and(data_water_mask, np.logical_not(data_severe_glint)) # 逻辑矩阵，true代表是水且不是耀斑 + valid_area_path = os.path.splitext(output_csvpath)[0] + "_valid_area" + write_bands(water_mask, valid_area_path, valid_area) + + x_out = [] + y_out = [] + with open(output_csvpath, "w") as f: + for row in range(dy): + for column in range(dx): + top_left = gdal.ApplyGeoTransform(geotransform_out, column, row) # (x,y) + bottom_right = gdal.ApplyGeoTransform(geotransform_out, column + 1, row + 1) + + top_left2 = gdal.ApplyGeoTransform(inv_geotransform_input, top_left[0], top_left[1]) # (x,y) + bottom_right2 = gdal.ApplyGeoTransform(inv_geotransform_input, bottom_right[0], bottom_right[1]) + + top_left2 = [int(num) for num in top_left2] + bottom_right2 = [int(num) for num in bottom_right2] + valid_area_local = valid_area[top_left2[1]:bottom_right2[1], top_left2[0]:bottom_right2[0]] # numpy索引时(y,x) + + re = np.where(valid_area_local == True) + coor_y_local, coor_x_local = getnearest(valid_area_local) + if re[0].shape[0] > 0: + x = top_left2[0] + int(re[1][0]) + y = top_left2[1] + int(re[0][0]) + + x = top_left2[0] + coor_x_local + y = top_left2[1] + coor_y_local + + # 将坐标和对应像元的光谱写入txt文件 + dest_coordinate = gdal.ApplyGeoTransform(geotransform_input, x, y) # (x,y) + f.write("%s\t %s\t %s\t %s\n" % (str(dest_coordinate[0]), str(dest_coordinate[1]), str(x), str(y))) + x_out.append(dest_coordinate[0]) + y_out.append(dest_coordinate[1]) + + return x_out, y_out + + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + parser = argparse.ArgumentParser(description="此程序通过给定间隔输出坐标。") + + parser.add_argument('-i1', '--water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + parser.add_argument('-i2', '--severe_glint', type=str, required=True, help='输入耀斑严重文件的路径') + parser.add_argument('-i3', '--interval', type=float, required=True, help='输入间隔，单位为米') + parser.add_argument('-o', '--output_csvpath', type=str, required=True, help='输出csv文件的路径') + + parser.add_argument('-v', '--verbose', action='store_true', help='启用详细模式') + + args = parser.parse_args() + + tmp = get_coor_base_interval(args.water_mask, args.severe_glint, args.output_csvpath, args.interval) diff --git a/get_spectral_in_coor.py b/get_spectral_in_coor.py new file mode 100644 index 0000000..cbe0721 --- /dev/null +++ b/get_spectral_in_coor.py @@ -0,0 +1,54 @@ +from util import * +from osgeo import gdal +import argparse +import numpy as np + + +def get_spectral_in_coor(imgpath, coorpath, outpath): + coor_spectral = np.loadtxt(coorpath, delimiter="\t") + + dataset = gdal.Open(imgpath) + im_width = dataset.RasterXSize # 栅格矩阵的列数 + im_height = dataset.RasterYSize # 栅格矩阵的行数 + num_bands = dataset.RasterCount # 栅格矩阵的波段数 + geotransform = dataset.GetGeoTransform() # 仿射矩阵 + im_proj = dataset.GetProjection() # 地图投影信息 + + new_columns = np.zeros((coor_spectral.shape[0], num_bands)) + coor_spectral = np.hstack((coor_spectral, new_columns)) + + for i in range(coor_spectral.shape[0]): + pixel_coor_x = int(coor_spectral[i, 2]) + pixel_coor_y = int(coor_spectral[i, 3]) + coor_spectral[i, 4:coor_spectral.shape[1]] = \ + dataset.ReadAsArray(int(pixel_coor_x), int(pixel_coor_y), 1, 1).reshape(num_bands) + + del dataset + + header = "image_coor_x\timage_coor_y\tpixel_coor_x\tpixel_coor_y\trest of colums are spectral" + + in_hdr_dict = spectral.envi.read_envi_header(get_hdr_file_path(imgpath)) + wavelengths = np.array(in_hdr_dict['wavelength']).astype('float64') + new_row = np.zeros((1, coor_spectral.shape[1])) + new_row[:, -num_bands:] = wavelengths + coor_spectral = np.vstack((new_row, coor_spectral)) + + np.savetxt(outpath, coor_spectral, fmt='%.4f', delimiter="\t", header=header) + # position_content_restore = np.loadtxt(outpath) + + return coor_spectral + + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + parser = argparse.ArgumentParser(description="此程序获取给定坐标的光谱曲线。") + + parser.add_argument('-i1', '--imgpath', type=str, required=True, help='输入影像文件的路径') + parser.add_argument('-i2', '--coorpath', type=str, required=True, help='输入坐标文件的路径') + parser.add_argument('-o', '--output_path', type=str, required=True, help='输出csv文件的路径') + + parser.add_argument('-v', '--verbose', action='store_true', help='启用详细模式') + + args = parser.parse_args() + + tmp = get_spectral_in_coor(args.imgpath, args.coorpath, args.output_path) diff --git a/interpolate.py b/interpolate.py new file mode 100644 index 0000000..bc710fd --- /dev/null +++ b/interpolate.py @@ -0,0 +1,77 @@ +from util import * +import numpy as np +from osgeo import gdal +import time +import argparse + + +@timeit +def interpolate_kriging_pykrige(x, y, z, proj, spatial_resolution, output_path=None): + # https://pygis.io/docs/e_interpolation.html + from pykrige.ok import OrdinaryKriging + + x_min = x.min() + x_max = x.max() + y_min = y.min() + y_max = y.max() + + # Generate a regular grid + step_x = int(x_max - x_min / spatial_resolution) + step_y = int(y_max - y_min / spatial_resolution) + grid_x = np.linspace(x_min, x_max, step_x) + grid_y = np.linspace(y_min, y_max, step_y) + + # Create ordinary kriging object: + OK = OrdinaryKriging( + x, + y, + z, + variogram_model="spherical", + verbose=False, + enable_plotting=False, + coordinates_type="euclidean", + ) + + # Execute on grid: + print("执行插值") + start1 = time.perf_counter() + + z1, ss1 = OK.execute("grid", grid_x, grid_y, backend="loop", n_closest_points=12) + + end1 = time.perf_counter() + print("插值，运行时间：", end1 - start1, "秒") + + if output_path is not None: + format_str = "GTiff" + driver = gdal.GetDriverByName(format_str) + # driver.Create()函数中RasterXSize代表影像的sample（列数），RasterYSize代表影像的line（行数） + dataset_out = driver.Create(output_path, step_x, step_y, 1, gdal.GDT_Float64) + + dataset_out.SetProjection(proj) + geotransform_out = (x_min, spatial_resolution, 0, y_max, 0, -spatial_resolution) + dataset_out.SetGeoTransform(geotransform_out) + + dataset_out.GetRasterBand(1).WriteArray(z1) + + del dataset_out + + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + parser = argparse.ArgumentParser(description="此程序基于反演结果进行插值。") + + parser.add_argument('-i1', '--pos_content_path', type=str, required=True, help='输入含量文件的路径') + parser.add_argument('-i2', '--img_path', type=str, required=True, help='输入影像文件的路径') + parser.add_argument('-i3', '--spatial_resolution', type=float, default=1.0, help='输出控件分辨率') + parser.add_argument('-o', '--outpath', required=True, type=str, help='输出文件的路径') + + parser.add_argument('-v', '--verbose', action='store_true', help='启用详细模式') + + args = parser.parse_args() + + pos_content = np.loadtxt(args.pos_content_path) + + dataset = gdal.Open(args.img_path) + im_proj = dataset.GetProjection() # 地图投影信息 + + interpolate_kriging_pykrige(pos_content[:, 0], pos_content[:, 1], pos_content[:, 2], im_proj, args.spatial_resolution, args.outpath) diff --git a/model_correction.py b/model_correction.py new file mode 100644 index 0000000..cc6d27c --- /dev/null +++ b/model_correction.py @@ -0,0 +1,461 @@ +import numpy as np +import os +from osgeo import gdal +from util import * +from type_define import * +import math +from pyproj import CRS +from pyproj import Transformer +import argparse +import json + + +def get_coor_base_lang_lat(img_path, water_mask, severe_glint, csv_pos, coor_type, radius, point_pos_strategy=PointPosStrategy.nearest_single, img_type=ImgType.ref): + """ + 此函数用于获取影像（img_path）在经纬度坐标（csv_pos）处的值，并输出到文件（out_txt）； + :param img_path: data ignore value应该为0；hdr不应该带bil等后缀 + :param img_type: 枚举ImgType + :param csv_pos: + :param coor_type: 枚举CoorType + :param radius: 单位为米 + :param point_pos_strategy: 枚举PointPosStrategy + :param severe_glint: 单波段图像，1代表耀斑严重的区域 + :return: + """ + dataset = gdal.Open(img_path) + im_width = dataset.RasterXSize # 栅格矩阵的列数 + im_height = dataset.RasterYSize # 栅格矩阵的行数 + num_bands = dataset.RasterCount # 栅格矩阵的波段数 + geotransform = dataset.GetGeoTransform() # 仿射矩阵 + im_proj = dataset.GetProjection() # 地图投影信息 + inv_geotransform_input = gdal.InvGeoTransform(geotransform) + + band_number_used = 1 # 找绿波段干啥？ + if img_type == ImgType.ref: + green_wave = 560 # 绿光信号最强 + band_number_used = find_band_number(green_wave, img_path) + elif img_type == ImgType.content: + band_number_used = 0 + + position_content = np.loadtxt(csv_pos, delimiter="\t") + new_columns_num = 4 # + num_bands + new_columns = np.zeros((position_content.shape[0], new_columns_num)) + position_content = np.hstack((position_content, new_columns)) + + dataset_water_mask = gdal.Open(water_mask) + # geotransform_water_mask = dataset_water_mask.GetGeoTransform() # 仿射矩阵 + # data_water_mask = dataset_water_mask.GetRasterBand(1).ReadAsArray() + + dataset_severe_glint = gdal.Open(severe_glint) + geotransform_severe_glint = dataset_severe_glint.GetGeoTransform() # 仿射矩阵 + data_severe_glint = dataset_severe_glint.GetRasterBand(1).ReadAsArray() + inv_geotransform_severe_glint = gdal.InvGeoTransform(geotransform_severe_glint) + + if coor_type == CoorType.latlong: # gps + epsg_code = 32600 + math.floor(position_content[0, 0] / 6) + 31 + crs = CRS.from_epsg(epsg_code) + projTransformer = Transformer.from_crs(crs.geodetic_crs, crs) + + for i in range(position_content.shape[0]): + if coor_type == CoorType.latlong: # gps + utm_cor = projTransformer.transform(position_content[i, 1], position_content[i, 0]) + elif coor_type == CoorType.utm: # utm + utm_cor = (position_content[i, 0], position_content[i, 1]) + + img_coor = gdal.ApplyGeoTransform(inv_geotransform_input, utm_cor[0], utm_cor[1]) # (x,y) + img_coor = [int(i) for i in img_coor] + radius_pixel_data = math.ceil(radius / geotransform[1] / 2) + data_local = dataset.GetRasterBand(band_number_used + 1).ReadAsArray(img_coor[0] - radius_pixel_data, + img_coor[1] - radius_pixel_data, + 2 * radius_pixel_data + 1, + 2 * radius_pixel_data + 1) + + water_mask_local = dataset_water_mask.GetRasterBand(1).ReadAsArray(img_coor[0] - radius_pixel_data, + img_coor[1] - radius_pixel_data, + 2 * radius_pixel_data + 1, + 2 * radius_pixel_data + 1) + + # if img_type == ImgType.content: + # data_local = dataset.ReadAsArray(img_coor[0] - radius_pixel_data, img_coor[1] - radius_pixel_data, + # 2 * radius_pixel_data + 1, 2 * radius_pixel_data + 1) + # elif img_type == ImgType.ref: + # data_local = dataset.GetRasterBand(band_number_used + 1).ReadAsArray(img_coor[0] - radius_pixel_data, + # img_coor[1] - radius_pixel_data, + # 2 * radius_pixel_data + 1, + # 2 * radius_pixel_data + 1) + + img_coor_severe_glint = gdal.ApplyGeoTransform(inv_geotransform_severe_glint, utm_cor[0], utm_cor[1]) + img_coor_severe_glint = [int(i) for i in img_coor_severe_glint] + if geotransform[1] == geotransform_severe_glint[1]: + radius_severe_glint = math.ceil(radius / geotransform_severe_glint[1] / 2) # + severe_glint_local = data_severe_glint[img_coor_severe_glint[1] - radius_severe_glint:img_coor_severe_glint[1] + radius_severe_glint + 1, + img_coor_severe_glint[0] - radius_severe_glint:img_coor_severe_glint[0] + radius_severe_glint + 1] + else: # 当分辨率不一致时，还没完成???????????????????????????????????????????????????????? + severe_glint_local = np.zeros(data_local.shape) + + for x in range(2 * radius_pixel_data + 1): + for y in range(2 * radius_pixel_data + 1): + img_coor_tmp = gdal.ApplyGeoTransform(geotransform, utm_cor[0], utm_cor[1]) + img_coor_tmp2 = gdal.ApplyGeoTransform(inv_geotransform_severe_glint, img_coor_tmp[0], img_coor_tmp[1]) + + severe_glint_local[y, x] = data_severe_glint[img_coor_tmp[1], img_coor_tmp[0]] + + invalid_value = 0 + data_local[severe_glint_local == 1] = invalid_value + data_local[water_mask_local == 0] = invalid_value + + # 根据PointPosStrategy确定使用哪些位置的数据 + if point_pos_strategy == PointPosStrategy.nearest_single: + coor_y_local, coor_x_local = getnearest(data_local, invalid_value) + + if coor_x_local is None: + continue + + coor_y = coor_y_local + img_coor[1] - radius_pixel_data + coor_x = coor_x_local + img_coor[0] - radius_pixel_data + + utm_x, utm_y = gdal.ApplyGeoTransform(geotransform, coor_x, coor_y) + + position_content[i, position_content.shape[1] - new_columns_num] = utm_x + position_content[i, position_content.shape[1] - new_columns_num + 1] = utm_y + + position_content[i, position_content.shape[1] - new_columns_num + 2] = int(coor_x) + position_content[i, position_content.shape[1] - new_columns_num + 3] = int(coor_y) + + # position_content[i, position_content.shape[1] - new_columns_num + 4:position_content.shape[1]] = \ + # dataset.ReadAsArray(int(coor_x), int(coor_y), 1, 1).reshape(num_bands) + + a = 1 + + elif point_pos_strategy == PointPosStrategy.four_quadrant: # 还没完成???????????????????????????????????????????????????????? + a = 1 + a = 1 + + rows_to_keep = np.where(position_content[:, position_content.shape[1] - new_columns_num + 3] != 0) + position_content2 = position_content[rows_to_keep] + + position_content3 = position_content2[~np.isnan(position_content2).any(axis=1)] + + # np.savetxt(out_txt, position_content3) + + del dataset, dataset_severe_glint, dataset_water_mask + + return position_content3 + + +def fit(x1, x2, y): + A = np.column_stack((x1, x2, np.ones((x2.shape[0], 1)))) + coefficients, _, _, _ = np.linalg.lstsq(A, y, rcond=None) + + return coefficients + + +def accuracy_evaluation(x1, x2, y_real, coefficients): + A = np.column_stack((x1, x2, np.ones((x2.shape[0], 1)))) + y_pred = A.dot(coefficients) + + accuracy = np.absolute((y_real - y_pred) / y_real * 100) + + return accuracy + + +def accuracy_evaluation_tss(x1, x2, y_real, coefficients): + A = np.column_stack((x1, x2, np.ones((x2.shape[0], 1)))) + y = A.dot(coefficients) + + y_pred = np.exp(y) + + accuracy = np.absolute((y_real - y_pred) / y_real * 100) + + return accuracy + + +def get_x_in_coor(coor, *args): + new_columns_counter = len(args) + new_columns = np.zeros((coor.shape[0], new_columns_counter)) + coor_extend = np.hstack((coor, new_columns)) + + for i in range(coor.shape[0]): + for j in range(new_columns_counter): + coor_extend[i, coor_extend.shape[1] - (new_columns_counter - j)] = args[j][ + int(coor_extend[i, coor_extend.shape[1] - new_columns_counter - 1]), + int(coor_extend[i, coor_extend.shape[1] - new_columns_counter - 2])] + + return coor_extend + + +def write_model_info(model_type, coefficients, accuracy, long, lat, outpath): + # 将 NumPy 数组转换为列表 + np_dict = { + 'model_type': model_type, + 'model_info': coefficients.tolist(), + 'accuracy': accuracy.tolist(), + 'long': long.tolist(), + 'lat': lat.tolist() + } + # 将字典写入 JSON 文件，使用 indent 参数进行格式化（每一级缩进4个空格） + with open(outpath, 'w') as f: + json.dump(np_dict, f, indent=4) + + +def chl_a(img_path, coor, outpath_coeff, window=5): # 叶绿素 + wave1 = 651 + band_651 = average_bands(wave1 - window, wave1 + window, img_path) + + wave2 = 707 + band_707 = average_bands(wave2 - window, wave2 + window, img_path) + + wave3 = 670 + band_670 = average_bands(wave3 - window, wave3 + window, img_path) + + x = (band_651 - band_707) / (band_707 - band_670) + + coor_extend = get_x_in_coor(coor, x) + + # 修正模型参数并输出 + x = coor_extend[:, -1] + y_real = coor_extend[:, 2] + + coefficients = np.polyfit(list(x), list(y_real), 1) + + y_pred = np.polyval(coefficients, list(x)) + accuracy = np.absolute((y_real - y_pred) / y_real * 100) + + write_model_info("chl-a", coefficients, accuracy, coor_extend[:, 0], coor_extend[:, 1], outpath_coeff) + model_type, model_info, accuracy_ = load_numpy_dict_from_json(outpath_coeff) + + return coefficients + + +def nh3(img_path, coor, outpath_coeff, window=5): # 氨氮 + wave1 = 600 + band_600 = average_bands(wave1 - window, wave1 + window, img_path) + + wave2 = 500 + band_500 = average_bands(wave2 - window, wave2 + window, img_path) + + wave3 = 850 + band_850 = average_bands(wave3 - window, wave3 + window, img_path) + + x13 = np.log(band_500 / band_850) + x23 = np.exp(band_600 / band_500) + + coor_extend = get_x_in_coor(coor, x13, x23) + + # 修正模型参数并输出 + x1 = coor_extend[:, -2] + x2 = coor_extend[:, -1] + y_real = coor_extend[:, 2] + coefficients = fit(x1, x2, y_real) + + accuracy = accuracy_evaluation(x1, x2, y_real, coefficients) + + write_model_info("nh3", coefficients, accuracy, coor_extend[:, 0], coor_extend[:, 1], outpath_coeff) + model_type, model_info, accuracy_ = load_numpy_dict_from_json(outpath_coeff) + + return coefficients + + +def mno4(img_path, coor, outpath_coeff, window=5): # 高猛酸盐 + wave1 = 500 + band_500 = average_bands(wave1 - window, wave1 + window, img_path) + + wave2 = 440 + band_440 = average_bands(wave2 - window, wave2 + window, img_path) + + wave3 = 610 + band_610 = average_bands(wave3 - window, wave3 + window, img_path) + + wave4 = 800 + band_800 = average_bands(wave4 - window, wave4 + window, img_path) + + x3 = band_500 / band_440 + x6 = band_610 / band_800 + + coor_extend = get_x_in_coor(coor, x3, x6) + + # 修正模型参数并输出 + x1 = coor_extend[:, -2] + x2 = coor_extend[:, -1] + y_real = coor_extend[:, 2] + coefficients = fit(x1, x2, y_real) + + accuracy = accuracy_evaluation(x1, x2, y_real, coefficients) + + write_model_info("mno4", coefficients, accuracy, coor_extend[:, 0], coor_extend[:, 1], outpath_coeff) + model_type, model_info, accuracy_ = load_numpy_dict_from_json(outpath_coeff) + + return coefficients + + +def tn(img_path, coor, outpath_coeff, window=5): # 总氮 + wave1 = 600 + band_600 = average_bands(wave1 - window, wave1 + window, img_path) + + wave2 = 500 + band_500 = average_bands(wave2 - window, wave2 + window, img_path) + + wave3 = 850 + band_850 = average_bands(wave3 - window, wave3 + window, img_path) + + x13 = np.log(band_500 / band_850) + x23 = np.exp(band_600 / band_500) + + coor_extend = get_x_in_coor(coor, x13, x23) + + # 修正模型参数并输出 + x1 = coor_extend[:, -2] + x2 = coor_extend[:, -1] + y_real = coor_extend[:, 2] + coefficients = fit(x1, x2, y_real) + + accuracy = accuracy_evaluation(x1, x2, y_real, coefficients) + + write_model_info("tn", coefficients, accuracy, coor_extend[:, 0], coor_extend[:, 1], outpath_coeff) + model_type, model_info, accuracy_ = load_numpy_dict_from_json(outpath_coeff) + + return coefficients + + +def tp(img_path, coor, outpath_coeff, window=5): # 总磷 + wave1 = 600 + band_600 = average_bands(wave1 - window, wave1 + window, img_path) + + wave2 = 500 + band_500 = average_bands(wave2 - window, wave2 + window, img_path) + + wave3 = 850 + band_850 = average_bands(wave3 - window, wave3 + window, img_path) + + x13 = np.log(band_500 / band_850) + x23 = np.exp(band_600 / band_500) + + coor_extend = get_x_in_coor(coor, x13, x23) + + # 修正模型参数并输出 + x1 = coor_extend[:, -2] + x2 = coor_extend[:, -1] + y_real = coor_extend[:, 2] + coefficients = fit(x1, x2, y_real) + + accuracy = accuracy_evaluation(x1, x2, y_real, coefficients) + + write_model_info("tp", coefficients, accuracy, coor_extend[:, 0], coor_extend[:, 1], outpath_coeff) + model_type, model_info, accuracy_ = load_numpy_dict_from_json(outpath_coeff) + + return coefficients + + +def tss(img_path, coor, outpath_coeff, window=5): # 总悬浮物？？？？？？？？？？？？？？？？？？ + wave1 = 555 + band_555 = average_bands(wave1 - window, wave1 + window, img_path) + + wave2 = 670 + band_670 = average_bands(wave2 - window, wave2 + window, img_path) + + wave3 = 490 + band_490 = average_bands(wave3 - window, wave3 + window, img_path) + + x1 = band_555 + band_670 + x2 = band_490 / band_555 + + coor_extend = get_x_in_coor(coor, x1, x2) + + # 修正模型参数并输出 + x1 = coor_extend[:, -2] + x2 = coor_extend[:, -1] + y_real = coor_extend[:, 2] + y = np.log(y_real) + coefficients = fit(x1, x2, y) + + accuracy = accuracy_evaluation_tss(x1, x2, y_real, coefficients) + + write_model_info("tss", coefficients, accuracy, coor_extend[:, 0], coor_extend[:, 1], outpath_coeff) + model_type, model_info, accuracy_ = load_numpy_dict_from_json(outpath_coeff) + + return coefficients + + +def main(): + parser = argparse.ArgumentParser(description="此程序用于通过实测数据修正模型参数。") + + # parser.add_argument("--global_arg", type=str, help="A global argument for all modes", required=True) + + # 创建子命令解析器 + subparsers = parser.add_subparsers(dest="algorithm", required=True, help="Choose a mode") + + chl_a_ = subparsers.add_parser("chl_a", help="叶绿素") + chl_a_.add_argument('-i1', '--img', type=str, required=True, help='输入影像文件的路径') + chl_a_.add_argument('-i2', '--water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + chl_a_.add_argument('-i3', '--severe_glint', type=str, required=True, help='输入耀斑严重文件的路径') + chl_a_.add_argument('-i4', '--measured_data', type=str, required=True, help='输入实测含量数据的路径') + chl_a_.add_argument('-i5', '--radius', type=float, default=2.0, help='输入实测坐标半径') + chl_a_.add_argument('-o', '--model_info_outpath', required=True, type=str, help='输出模型信息文件的路径') + chl_a_.set_defaults(func=chl_a) + + nh3_ = subparsers.add_parser("nh3", help="氨氮") + nh3_.add_argument('-i1', '--img', type=str, required=True, help='输入影像文件的路径') + nh3_.add_argument('-i2', '--water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + nh3_.add_argument('-i3', '--severe_glint', type=str, required=True, help='输入耀斑严重文件的路径') + nh3_.add_argument('-i4', '--measured_data', type=str, required=True, help='输入实测含量数据的路径') + nh3_.add_argument('-i5', '--radius', type=float, default=2.0, help='输入实测坐标半径') + nh3_.add_argument('-o', '--model_info_outpath', required=True, type=str, help='输出模型信息文件的路径') + nh3_.set_defaults(func=nh3) + + mno4_ = subparsers.add_parser("mno4", help="高猛酸盐") + mno4_.add_argument('-i1', '--img', type=str, required=True, help='输入影像文件的路径') + mno4_.add_argument('-i2', '--water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + mno4_.add_argument('-i3', '--severe_glint', type=str, required=True, help='输入耀斑严重文件的路径') + mno4_.add_argument('-i4', '--measured_data', type=str, required=True, help='输入实测含量数据的路径') + mno4_.add_argument('-i5', '--radius', type=float, default=2.0, help='输入实测坐标半径') + mno4_.add_argument('-o', '--model_info_outpath', required=True, type=str, help='输出模型信息文件的路径') + mno4_.set_defaults(func=mno4) + + tn_ = subparsers.add_parser("tn", help="总氮") + tn_.add_argument('-i1', '--img', type=str, required=True, help='输入影像文件的路径') + tn_.add_argument('-i2', '--water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + tn_.add_argument('-i3', '--severe_glint', type=str, required=True, help='输入耀斑严重文件的路径') + tn_.add_argument('-i4', '--measured_data', type=str, required=True, help='输入实测含量数据的路径') + tn_.add_argument('-i5', '--radius', type=float, default=2.0, help='输入实测坐标半径') + tn_.add_argument('-o', '--model_info_outpath', required=True, type=str, help='输出模型信息文件的路径') + tn_.set_defaults(func=tn) + + tp_ = subparsers.add_parser("tp", help="总磷") + tp_.add_argument('-i1', '--img', type=str, required=True, help='输入影像文件的路径') + tp_.add_argument('-i2', '--water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + tp_.add_argument('-i3', '--severe_glint', type=str, required=True, help='输入耀斑严重文件的路径') + tp_.add_argument('-i4', '--measured_data', type=str, required=True, help='输入实测含量数据的路径') + tp_.add_argument('-i5', '--radius', type=float, default=2.0, help='输入实测坐标半径') + tp_.add_argument('-o', '--model_info_outpath', required=True, type=str, help='输出模型信息文件的路径') + tp_.set_defaults(func=tp) + + tss_ = subparsers.add_parser("tss", help="总悬浮物") + tss_.add_argument('-i1', '--img', type=str, required=True, help='输入影像文件的路径') + tss_.add_argument('-i2', '--water_mask', type=str, required=True, help='输入水域掩膜文件的路径') + tss_.add_argument('-i3', '--severe_glint', type=str, required=True, help='输入耀斑严重文件的路径') + tss_.add_argument('-i4', '--measured_data', type=str, required=True, help='输入实测含量数据的路径') + tss_.add_argument('-i5', '--radius', type=float, default=2.0, help='输入实测坐标半径') + tss_.add_argument('-o', '--model_info_outpath', required=True, type=str, help='输出模型信息文件的路径') + tss_.set_defaults(func=tss) + + # 解析参数 + args = parser.parse_args() + tmp = get_coor_base_lang_lat(args.img, args.water_mask, args.severe_glint, args.measured_data, CoorType.latlong, + args.radius) + if args.algorithm == "chl_a": + args.func(args.img, tmp, args.model_info_outpath) + elif args.algorithm == "nh3": + args.func(args.img, tmp, args.model_info_outpath) + elif args.algorithm == "mno4": + args.func(args.img, tmp, args.model_info_outpath) + elif args.algorithm == "tn": + args.func(args.img, tmp, args.model_info_outpath) + elif args.algorithm == "tp": + args.func(args.img, tmp, args.model_info_outpath) + elif args.algorithm == "tss": + args.func(args.img, tmp, args.model_info_outpath) + + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + main() diff --git a/readme.txt b/readme.txt new file mode 100644 index 0000000..04d18bd --- /dev/null +++ b/readme.txt @@ -0,0 +1,214 @@ +1、环境搭建注意事项 +python版本：3.13.1 +gdal whl文件下载路径：https://github.com/cgohlke/geospatial-wheels/releases +安装PyKrige时，pip换成国内源，否则会失败：https://blog.csdn.net/weixin_50679163/article/details/122392249 + +2、命令行工具使用说明 +2.1 数据要求 +（1）输入原始影像是反射率数据 +（2）无效像元对应的值必须为0，data ignore value = 0 +（3）输入影像必须为bsq +（4）实测数据的格式：以tab分割的3列：经度，纬度，含量 + +2.2 二级指令 +水体处理分为8个步骤，每步都是可调用的单独的程序，但是由于某些步骤的实现方式不止一种，所以引入了二级命令隔离参数。 +使用二级指令的程序：水域掩膜、去耀斑、基于实测值进行模型修正、反演。 + +3、命令行程序文档 +3.1 水域掩膜 +（1）此程序用于提取水域区域，输出的水域栅格和输入的影像具有相同的行列数。支持2种二级指令（算法）：rasterize_shp，ndwi。 + +（2）rasterize_shp + + -h, --help show this help message and exit + -i1, --img_path IMG_PATH + 输入影像文件的路径 + -i2, --shp_path SHP_PATH + 输入shp文件的路径 + -o, --water_mask_outpath WATER_MASK_OUTPATH + 输出水体掩膜文件的路径 + +命令示例： +rasterize_shp -i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m -i2 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_cutarea.shp -o D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_water_mask + +（3）ndwi + +options: + -h, --help show this help message and exit + -i1, --img_path IMG_PATH + 输入影像文件的路径 + -i2, --ndwi_threshold NDWI_THRESHOLD + 输入ndwi水体阈值，大于此值的为水域 + -o, --water_mask_outpath WATER_MASK_OUTPATH + 输出水体掩膜文件的路径 + +命令示例： +ndwi -i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m -i2 0.4 -o D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_water_mask_ndwi + +3.2 找耀斑严重区域 + +此程序通过大律法分割图像，提取耀斑最严重的区域，输出的栅格和输入的影像具有相同的行列数。 + +options: + -h, --help show this help message and exit + -i1, --input INPUT 输入影像文件的路径 + -i2, --input_water_mask INPUT_WATER_MASK + 输入水域掩膜文件的路径 + -gw, --glint_wave GLINT_WAVE + 用于提取耀斑严重区域的波段. + -o, --output OUTPUT 输出文件的路径 + +命令示例： +-i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_bsq -i2 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_water_mask -o D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_severe_glint + +3.3 去耀斑 +（1）此程序用于去除水域数据的耀斑，暂时支持3种算法。支持3种二级命令（算法）：subnir,regression_slope,oxygen_absorption + +（2）subnir +options: + -h, --help show this help message and exit + -i1, --input INPUT 输入影像文件的路径 + -i2, --input_water_mask INPUT_WATER_MASK + 输入水域掩膜文件的路径 + -o, --output OUTPUT 输出文件的路径 + -s, --start_wave START_WAVE + nir开始波段 + -e, --end_wave END_WAVE + nir结束波段 + +命令示例： +subnir -i1 D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510 -o D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510_nir_730-750 -s 730 -e 750 +subnir -i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_bsq -i2 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_water_mask -o D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_nir_730-750 -s 730 -e 750 + +（3）regression_slope +options: + -h, --help show this help message and exit + -i1, --input INPUT 输入影像文件的路径 + -i2, --input_water_mask INPUT_WATER_MASK + 输入水域掩膜文件的路径 + -o, --output OUTPUT 输出文件的路径 + -s, --start_wave START_WAVE + nir开始波段 + -e, --end_wave END_WAVE + nir结束波段 + -r, --roi ROI 输入roi文件的路径 + +命令示例： +regression_slope -i1 D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510 -o D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510_slope -s 730 -e 750 -r D:\PycharmProjects\sun_glint\test_data\roi2.xml +regression_slope -i1 D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510_wgs84.dat -o D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510_slope -s 730 -e 750 -r D:\PycharmProjects\sun_glint\test_data\roi2_latlong.xml +regression_slope -i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_bsq -i2 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_water_mask -o D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_slope -s 730 -e 750 -r D:\PycharmProjects\0water_rlx\test_data\roi_glint.xml + +（4）oxygen_absorption +options: + -h, --help show this help message and exit + -i1, --input INPUT 输入影像文件的路径 + -i2, --input_water_mask INPUT_WATER_MASK + 输入水域掩膜文件的路径 + -o, --output OUTPUT 输出文件的路径 + -l, --left_shoulder_wave LEFT_SHOULDER_WAVE + 氧气吸收谷左肩波长 + -v, --valley_wave VALLEY_WAVE + 氧气吸收谷底波长 + -r, --right_shoulder_wave RIGHT_SHOULDER_WAVE + 氧气吸收谷右肩波长 + -sw, --shoulder_window SHOULDER_WINDOW + 氧气吸收谷左右肩平均窗口，半径 + -vw, --valley_window VALLEY_WINDOW + 氧气吸收谷底平均窗口，半径 + +命令示例： +oxygen_absorption -i1 D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510 -o D:\PycharmProjects\sun_glint\test_data\ref_deepWater_20230914_154510_oxygen -l 752.49 -v 762.49 -r 772.49 -sw 1 -vw 1 +oxygen_absorption -i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_bsq -i2 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_water_mask -o D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_oxygen -l 752.49 -v 762.49 -r 772.49 -sw 1 -vw 1 + +3.4 基于实测值进行模型修正 +（1）此程序用于通过实测数据修正模型参数。支持6种二级命令（算法）：chl_a,nh3,mno4,tn,tp,tss。 + +（2）6种二级命令（算法）的参数都一样 + +options: + -h, --help show this help message and exit + -i1, --img IMG 输入影像文件的路径 + -i2, --water_mask WATER_MASK + 输入水域掩膜文件的路径 + -i3, --severe_glint SEVERE_GLINT + 输入耀斑严重文件的路径 + -i4, --measured_data MEASURED_DATA + 输入实测含量数据的路径 + -i5, --radius RADIUS 输入实测坐标半径 + -o, --model_info_outpath MODEL_INFO_OUTPATH + 输出模型信息文件的路径 + +命令示例： +chl_a -i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_bsq -i2 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_water_mask -i3 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_severe_glint -i4 D:\PycharmProjects\0water_rlx\test_data\content_retrieval\chl-a\Measured_point_pos_chl-a.txt -o D:\PycharmProjects\0water_rlx\test_data\content_retrieval\chl-a\chl-a.json + +3.5 获取等间隔点坐标 + +此程序通过给定间隔输出坐标。 + +options: + -h, --help show this help message and exit + -i1, --water_mask WATER_MASK + 输入水域掩膜文件的路径 + -i2, --severe_glint SEVERE_GLINT + 输入耀斑严重文件的路径 + -i3, --interval INTERVAL + 输入间隔，单位为米 + -o, --output_csvpath OUTPUT_CSVPATH + 输出csv文件的路径 + +命令示例： +-i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_water_mask -i2 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_severe_glint -i3 100 -o D:\PycharmProjects\0water_rlx\test_data\coor_interval.csv + +3.6 基于（5）所得坐标获取光谱 + +此程序获取给定坐标的光谱曲线。 + +options: + -h, --help show this help message and exit + -i1, --imgpath IMGPATH + 输入影像文件的路径 + -i2, --coorpath COORPATH + 输入坐标文件的路径 + -o, --output_path OUTPUT_PATH + 输出csv文件的路径 + +命令示例： +-i1 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_bsq -i2 D:\PycharmProjects\0water_rlx\test_data\coor_interval.csv -o D:\PycharmProjects\0water_rlx\test_data\coor-spectral.csv + +3.7 反演 + +（1）此程序使用修正后的模型进行反演。支持6种二级命令（算法）：chl_a,nh3,mno4,tn,tp,tss。 + +（2）6种二级命令（算法）的参数都一样 + +options: + -h, --help show this help message and exit + -i1, --model_info_path MODEL_INFO_PATH + 输入模型信息文件的路径 + -i2, --coor_spectral_path COOR_SPECTRAL_PATH + 输入坐标-光谱文件的路径 + -i3, --wave_radius WAVE_RADIUS + 输入波长平均半径 + -o, --outpath OUTPATH + 输出文件的路径 + +命令示例： +chl_a -i1 D:\PycharmProjects\0water_rlx\test_data\content_retrieval\chl-a\chl-a.json -i2 D:\PycharmProjects\0water_rlx\test_data\coor-spectral.csv -o D:\PycharmProjects\0water_rlx\test_data\coor-spectral-retrieval.csv + +3.8 插值 + +此程序基于反演结果进行插值。 + +options: + -h, --help show this help message and exit + -i1, --pos_content_path POS_CONTENT_PATH + 输入含量文件的路径 + -i2, --img_path IMG_PATH + 输入影像文件的路径 + -i3, --spatial_resolution SPATIAL_RESOLUTION + 输出控件分辨率 + -o, --outpath OUTPATH + 输出文件的路径 + +命令示例： +-i1 D:\PycharmProjects\0water_rlx\test_data\coor-spectral-retrieval.csv -i2 D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m -o D:\PycharmProjects\0water_rlx\test_data\content.tiff diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000000000000000000000000000000000000..c9190eb04e9f9bad2fb5bf809d60daad6ef5ca11 GIT binary patch literal 302 zcmYL@Q3}E^5Jcx&@G2p#(SCUXf4qW(YAe--wjy3$o!LYr><0GD&g|zM(PF}a9%t5` zuZCv@I8^-A1J+pKjt6c?comf^RP4y$jzU# zHkGL^@0pGZvlIKsH*lha23ylO-kWtDPE@s2{BQKs>uE|+XbkxW} K38KL+$q8TRo+{S> literal 0 HcmV?d00001 diff --git a/retrieval.py b/retrieval.py new file mode 100644 index 0000000..0d1ea88 --- /dev/null +++ b/retrieval.py @@ -0,0 +1,169 @@ +import numpy as np +import os +from osgeo import gdal +from util import * +from type_define import * +import math +from pyproj import CRS +from pyproj import Transformer +import argparse + + +def find_index(wavelength, array): + differences = np.abs(array - wavelength) + min_position = np.argmin(differences) + return min_position + + +def get_mean_value(index, array, window): + window = int(window) + result = array[1:, index - window:index + window + 1].mean(axis=1) + + return result + + +def calculate(x1, x2, coefficients): + A = np.column_stack((x1, x2, np.ones((x2.shape[0], 1)))) + y_pred = A.dot(coefficients) + + return y_pred + + +def retrieval_chl_a(model_info_path, coor_spectral_path, output_path, window=5): + model_type, model_info, accuracy_ = load_numpy_dict_from_json(model_info_path) + coor_spectral = np.loadtxt(coor_spectral_path) + + wave1 = 651 + index_wave1 = find_index(wave1, coor_spectral[0, :]) + band_651 = get_mean_value(index_wave1, coor_spectral, window) + + wave2 = 707 + index_wave2 = find_index(wave2, coor_spectral[0, :]) + band_707 = get_mean_value(index_wave2, coor_spectral, window) + + wave3 = 670 + index_wave3 = find_index(wave3, coor_spectral[0, :]) + band_670 = get_mean_value(index_wave3, coor_spectral, window) + + x = (band_651 - band_707) / (band_707 - band_670) + retrieval_result = np.polyval(model_info, list(x)) + + position_content = np.hstack((coor_spectral[1:, 0:2], retrieval_result.reshape((retrieval_result.shape[0], 1)))) + np.savetxt(output_path, position_content, fmt='%.4f', delimiter="\t") + + return position_content + + +def retrieval_nh3(model_info_path, coor_spectral_path, output_path=None, window=5): + model_type, model_info, accuracy_ = load_numpy_dict_from_json(model_info_path) + coor_spectral = np.loadtxt(coor_spectral_path) + + wave1 = 600 + index_wave1 = find_index(wave1, coor_spectral[0, :]) + band_600 = get_mean_value(index_wave1, coor_spectral, window) + + wave2 = 500 + index_wave2 = find_index(wave2, coor_spectral[0, :]) + band_500 = get_mean_value(index_wave2, coor_spectral, window) + + wave3 = 850 + index_wave3 = find_index(wave3, coor_spectral[0, :]) + band_850 = get_mean_value(index_wave3, coor_spectral, window) + + x13 = np.log(band_500 / band_850) + x23 = np.exp(band_600 / band_500) + retrieval_result = calculate(x13, x23, model_info) + + position_content = np.hstack((coor_spectral[1:, 0:2], retrieval_result.reshape((retrieval_result.shape[0], 1)))) + if output_path is not None: + np.savetxt(output_path, position_content, fmt='%.4f', delimiter="\t") + + return position_content + + +def retrieval_tss(model_info_path, coor_spectral_path, output_path, window=5): + position_content = retrieval_nh3(model_info_path, coor_spectral_path, window=window) + + tmp = np.exp(position_content[:, -1]) + position_content[:, -1] = tmp + + np.savetxt(output_path, position_content, fmt='%.4f', delimiter="\t") + + return position_content + + +def main(): + parser = argparse.ArgumentParser(description="此程序使用修正后的模型进行反演。") + + # parser.add_argument("--global_arg", type=str, help="A global argument for all modes", required=True) + + # 创建子命令解析器 + subparsers = parser.add_subparsers(dest="algorithm", required=True, help="Choose a mode") + + chl_a_ = subparsers.add_parser("chl_a", help="叶绿素") + chl_a_.add_argument('-i1', '--model_info_path', type=str, required=True, help='输入模型信息文件的路径') + chl_a_.add_argument('-i2', '--coor_spectral_path', type=str, required=True, + help='输入坐标-光谱文件的路径') + chl_a_.add_argument('-i3', '--wave_radius', type=float, default=5.0, help='输入波长平均半径') + chl_a_.add_argument('-o', '--outpath', required=True, type=str, help='输出文件的路径') + chl_a_.set_defaults(func=retrieval_chl_a) + + nh3_ = subparsers.add_parser("nh3", help="氨氮") + nh3_.add_argument('-i1', '--model_info_path', type=str, required=True, help='输入模型信息文件的路径') + nh3_.add_argument('-i2', '--coor_spectral_path', type=str, required=True, + help='输入坐标-光谱文件的路径') + nh3_.add_argument('-i3', '--wave_radius', type=float, default=5.0, help='输入波长平均半径') + nh3_.add_argument('-o', '--outpath', required=True, type=str, help='输出文件的路径') + nh3_.set_defaults(func=retrieval_nh3) + + mno4_ = subparsers.add_parser("mno4", help="高猛酸盐") + mno4_.add_argument('-i1', '--model_info_path', type=str, required=True, help='输入模型信息文件的路径') + mno4_.add_argument('-i2', '--coor_spectral_path', type=str, required=True, + help='输入坐标-光谱文件的路径') + mno4_.add_argument('-i3', '--wave_radius', type=float, default=5.0, help='输入波长平均半径') + mno4_.add_argument('-o', '--outpath', required=True, type=str, help='输出文件的路径') + mno4_.set_defaults(func=retrieval_nh3) + + tn_ = subparsers.add_parser("tn", help="总氮") + tn_.add_argument('-i1', '--model_info_path', type=str, required=True, help='输入模型信息文件的路径') + tn_.add_argument('-i2', '--coor_spectral_path', type=str, required=True, + help='输入坐标-光谱文件的路径') + tn_.add_argument('-i3', '--wave_radius', type=float, default=5.0, help='输入波长平均半径') + tn_.add_argument('-o', '--outpath', required=True, type=str, help='输出文件的路径') + tn_.set_defaults(func=retrieval_nh3) + + tp_ = subparsers.add_parser("tp", help="总磷") + tp_.add_argument('-i1', '--model_info_path', type=str, required=True, help='输入模型信息文件的路径') + tp_.add_argument('-i2', '--coor_spectral_path', type=str, required=True, + help='输入坐标-光谱文件的路径') + tp_.add_argument('-i3', '--wave_radius', type=float, default=5.0, help='输入波长平均半径') + tp_.add_argument('-o', '--outpath', required=True, type=str, help='输出文件的路径') + tp_.set_defaults(func=retrieval_nh3) + + tss_ = subparsers.add_parser("tss", help="总悬浮物") + tss_.add_argument('-i1', '--model_info_path', type=str, required=True, help='输入模型信息文件的路径') + tss_.add_argument('-i2', '--coor_spectral_path', type=str, required=True, + help='输入坐标-光谱文件的路径') + tss_.add_argument('-i3', '--wave_radius', type=float, default=5.0, help='输入波长平均半径') + tss_.add_argument('-o', '--outpath', required=True, type=str, help='输出文件的路径') + tss_.set_defaults(func=retrieval_tss) + + # 解析参数 + args = parser.parse_args() + if args.algorithm == "chl_a": + args.func(args.model_info_path, args.coor_spectral_path, args.outpath, args.wave_radius) + elif args.algorithm == "nh3": + args.func(args.model_info_path, args.coor_spectral_path, args.outpath, args.wave_radius) + elif args.algorithm == "mno4": + args.func(args.model_info_path, args.coor_spectral_path, args.outpath, args.wave_radius) + elif args.algorithm == "tn": + args.func(args.model_info_path, args.coor_spectral_path, args.outpath, args.wave_radius) + elif args.algorithm == "tp": + args.func(args.model_info_path, args.coor_spectral_path, args.outpath, args.wave_radius) + elif args.algorithm == "tss": + args.func(args.model_info_path, args.coor_spectral_path, args.outpath, args.wave_radius) + + +# Press the green button in the gutter to run the script. +if __name__ == '__main__': + main() diff --git a/type_define.py b/type_define.py new file mode 100644 index 0000000..b06e6a1 --- /dev/null +++ b/type_define.py @@ -0,0 +1,22 @@ +from enum import Enum, unique + +class FlareModel(Enum): + otsu = 0 + threshold = 1 + img = 2 + + +class ImgType(Enum): + ref = 0 + content = 1 + + +# @unique +class CoorType(Enum): + latlong = 0 + utm = 1 + + +class PointPosStrategy(Enum): + nearest_single = 0 + four_quadrant = 1 diff --git a/util.py b/util.py new file mode 100644 index 0000000..8b92205 --- /dev/null +++ b/util.py @@ -0,0 +1,174 @@ +import os, spectral +import time +import numpy as np +from osgeo import gdal +from enum import Enum, unique +import math +import json + + +class CoorType(Enum): + depend_on_image = 0 # 影像是啥类型坐标就是啥坐标 + pixel = 1 + + +class Timer: # Context Manager + def __enter__(self): + self.start = time.time() + return self + + def __exit__(self, exc_type, exc_value, traceback): + print(exc_type, exc_value, traceback) + print(f"Run Time: {time.time() - self.start}") + + +def timeit(f): # decorator + def wraper(*args, **kwargs): + start = time.time() + ret = f(*args, **kwargs) + print(f.__name__ + " run time: " + str(round(time.time() - start, 2)) + " s.") + return ret + + return wraper + + +def get_hdr_file_path(file_path): + return os.path.splitext(file_path)[0] + ".hdr" + + +def find_band_number(wav1, imgpath): + in_hdr_dict = spectral.envi.read_envi_header(get_hdr_file_path(imgpath)) + + wavelengths = np.array(in_hdr_dict['wavelength']).astype('float64') + + differences = np.abs(wavelengths - wav1) + min_position = np.argmin(differences) + + return int(min_position) + + +@timeit +def average_bands(start_wave, end_wave, imgpath): + start_bandnumber = find_band_number(start_wave, imgpath) + end_bandnumber = find_band_number(end_wave, imgpath) + + dataset = gdal.Open(imgpath) + + averaged_band = 1 + for i in range(start_bandnumber, end_bandnumber + 1): + if i == start_bandnumber: + averaged_band = dataset.GetRasterBand(i + 1).ReadAsArray() + else: + tmp = dataset.GetRasterBand(i + 1).ReadAsArray() + averaged_band = (averaged_band + tmp) / 2 + + del dataset + + return averaged_band + + +def exclude_by_mask(band, water_mask_path, ignore_value=0): + dataset = gdal.Open(water_mask_path) + data_tmp = dataset.GetRasterBand(1).ReadAsArray() + + del dataset + + band[np.where(data_tmp == ignore_value)] = 0 + + return band + + +@timeit +def average_bands_in_mask(start_wave, end_wave, imgpath, water_mask_path): + tmp = average_bands(start_wave, end_wave, imgpath) + + tmp = exclude_by_mask(tmp, water_mask_path) + + # raster_fn_out_tmp = append2filename(imgpath, "glint_delete") + # write_bands(imgpath, raster_fn_out_tmp, tmp) + + return tmp + + +def get_average_value(dataset, x, y, band_number, window): + spectral_tmp = dataset.ReadAsArray(x, y, 1, 1) + average_value = spectral_tmp[band_number - window:band_number + window, :, :].mean() + + return average_value + + +def get_valid_extent(dataset, data_ignore_value=0): + pass + + +def write_bands(imgpath_in, imgpath_out, *args): + # 将输入的波段（可变）写入文件 + dataset = gdal.Open(imgpath_in) + im_width = dataset.RasterXSize + im_height = dataset.RasterYSize + num_bands = dataset.RasterCount + geotransform = dataset.GetGeoTransform() + im_proj = dataset.GetProjection() + + format = "ENVI" + driver = gdal.GetDriverByName(format) + dst_ds = driver.Create(imgpath_out, im_width, im_height, len(args), gdal.GDT_Float32, + options=["INTERLEAVE=BSQ"]) + dst_ds.SetGeoTransform(geotransform) + dst_ds.SetProjection(im_proj) + + for i in range(len(args)): + dst_ds.GetRasterBand(i + 1).WriteArray(args[i]) + + del dataset, dst_ds + + +def append2filename(file_path, txt2add): + imgfile_out_tmp = os.path.splitext(file_path) + new_file_path = imgfile_out_tmp[0] + "_" + txt2add + imgfile_out_tmp[1] + + return new_file_path + + +def write_fields_to_hdrfile(source_hdr_file, dest_hdr_file): + source_fields = spectral.envi.read_envi_header(source_hdr_file) + dest_fields = spectral.envi.read_envi_header(dest_hdr_file) + + with open(dest_hdr_file, "a", encoding='utf-8') as f: + for key in source_fields.keys(): + if key in dest_fields or key == "description": + continue + + if key == "data ignore value" or key == "wavelength" or key == "wavelength units": + if type(source_fields[key]) == list: + f.write(key + " = {" + ", ".join(source_fields[key]) + "}\n") + else: + f.write(key + " = " + source_fields[key] + "\n") + + +def getnearest(m, invalid_value=0): + layer_number = math.floor(m.shape[0] / 2) + center = layer_number + + for i in range(layer_number + 1): + orig = (center - i, center - i) + + tmp = m[center - i:center + i + 1, center - i:center + i + 1] + valid_indices = np.where((tmp != invalid_value) & np.isfinite(tmp)) + + if valid_indices[0].shape[0] != 0: + return int(valid_indices[0][0] + orig[0]), int(valid_indices[1][0] + orig[0]) # (y ,x) + + return None, None + + +def load_numpy_dict_from_json(filename): + with open(filename, 'r') as f: + np_dict = json.load(f) + + # 将字典中的列表转换回 NumPy 数组 + model_type = np_dict['model_type'] + model_info = np.array(np_dict['model_info']) + precision = np.array(np_dict['accuracy']) + + return model_type, model_info, precision diff --git a/算法步骤.txt b/算法步骤.txt new file mode 100644 index 0000000..d6380d7 --- /dev/null +++ b/算法步骤.txt @@ -0,0 +1,24 @@ +20241212 +要求： +（1）输入原始影像是反射率数据，每步都是可调用的单独的程序 +（2）无效像元对应的值必须为0，data ignore value = 0 +（3）输入影像必须为bsq +（4）实测数据的格式为3列：经度，纬度，含量 + + +1 提取准确水体范围（1）自动实现ndwi阈值（会出现提取错误，例如孤岛现象）（2）矢量shp文件转栅格 +输出的水域掩膜和影像有相同行列数和分辨率 + +2 otsu找过曝/耀斑严重区域（是否通过dn值自动识别） + +3 去耀斑的方法：添加水域掩膜 + +4 修正模型参数：（1）通过实测数据（避过过曝/耀斑严重的位置）修正模型参数，【任总决定不反演】，现在还是通过envi做的（2）模型是否合适（3）csv规定格式 json + （4）输出修正模型信息（怎么输出？）（5）随机选择几个点做精度评价，需要将实测数据分为反演集合和测试集合 +输出的模型信息为json格式：模型类型、模型参数、模型精度（带坐标） + +5 根据输入（原始影像）间隔（先像素后距离）和取点策略（仅实现一种） 输出一系列是水体且【避过过曝/耀斑严重的位置】的像素 的csv坐标（修改） + +6 根据上步坐标，获取光谱和其真实坐标位置 +7 对第4、6步数据进行计算，输出第5步关键点的反演值 +8 插值：和arcgis的结果有差异