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)