实现水体含量反演程序。

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()

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))

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))

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))

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()

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)

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)

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)

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)

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()

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

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()

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

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

算法步骤.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的结果有差异