"""
@Project ：microproduct
@File ：ImageHandle.py
@Function ：实现对待处理SAR数据的读取、格式标准化和处理完后保存文件功能
@Author ：LMM
@Date ：2021/10/19 14:39
@Version ：1.0.0
"""
import os
from PIL import Image
import time
from osgeo import gdal
from osgeo import osr
import numpy as np
from PIL import Image
from xml.etree.ElementTree import ElementTree, Element
import math
import shutil
import tarfile



class ImageHandler:

    @staticmethod
    def extract_tarfile(self, tar_file, target_dir):
        with tarfile.open(tar_file, 'r') as tar_ref:
            for member in tar_ref.getmembers():
                if member.isdir():
                    continue
                filename = os.path.basename(member.name)
                source = tar_ref.extractfile(member)
                target = open(os.path.join(target_dir, filename), "wb")
                with source, target:
                    shutil.copyfileobj(source, target)

    """
    影像读取、编辑、保存
    """

    def __init__(self):
        pass

    @staticmethod
    def get_dataset(filename):
        """
       :param filename: tif路径
       :return: 图像句柄
       """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        return dataset

    def get_scope(self, filename):
        """
        :param filename: tif路径
        :return: 图像范围
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        im_scope = self.cal_img_scope(dataset)
        del dataset
        dataset=None
        return im_scope

    @staticmethod
    def get_projection(filename):
        """
         :param filename: tif路径
         :return: 地图投影信息
         """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        im_proj = dataset.GetProjection()
        del dataset
        dataset=None
        return im_proj

    @staticmethod
    def get_geotransform(filename):
        """
        :param filename: tif路径
        :return: 从图像坐标空间（行、列），也称为（像素、线）到地理参考坐标空间（投影或地理坐标）的仿射变换
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        geotransform = dataset.GetGeoTransform()
        del dataset
        dataset=None
        return geotransform

    def get_invgeotransform(filename):
        """
        :param filename: tif路径
        :return: 从地理参考坐标空间（投影或地理坐标）的到图像坐标空间（行、列
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        geotransform = dataset.GetGeoTransform()
        geotransform = gdal.InvGeoTransform(geotransform)
        del dataset
        dataset=None
        return geotransform

    @staticmethod
    def get_bands(filename):
        """
        :param filename: tif路径
        :return: 影像的波段数
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        bands = dataset.RasterCount
        del dataset
        dataset=None
        return bands

    @staticmethod
    def geo2lonlat(dataset, x, y):
        """
        将投影坐标转为经纬度坐标（具体的投影坐标系由给定数据确定）
        :param dataset: GDAL地理数据
        :param x: 投影坐标x
        :param y: 投影坐标y
        :return: 投影坐标(x, y)对应的经纬度坐标(lon, lat)
        """
        prosrs = osr.SpatialReference()
        prosrs.ImportFromWkt(dataset.GetProjection())
        geosrs = prosrs.CloneGeogCS()
        ct = osr.CoordinateTransformation(prosrs, geosrs)
        coords = ct.TransformPoint(x, y)
        return coords[:2]

    @staticmethod
    def get_band_array(filename, num=1):
        """
        :param filename: tif路径
        :param num: 波段序号
        :return: 对应波段的矩阵数据
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        bands = dataset.GetRasterBand(num)
        array = bands.ReadAsArray(0, 0, bands.XSize, bands.YSize)

        # if 'int' in str(array.dtype):
        #     array[np.where(array == -9999)] = np.inf
        # else:
        #     array[np.where(array < -9000.0)] = np.nan

        del dataset
        dataset=None
        return array

    @staticmethod
    def write_imgArray(filename, im_data, no_data='null'):
        """
        影像保存
        :param filename: 保存的路径
        :param im_proj:
        :param im_geotrans:
        :param im_data:
        :param no_data: 把无效值设置为 nodata
        :return:
        """

        gdal_dtypes = {
            'int8': gdal.GDT_Byte,
            'unit16': gdal.GDT_UInt16,
            'int16': gdal.GDT_Int16,
            'unit32': gdal.GDT_UInt32,
            'int32': gdal.GDT_Int32,
            'float32': gdal.GDT_Float32,
            'float64': gdal.GDT_Float64,
        }
        if not gdal_dtypes.get(im_data.dtype.name, None) is None:
            datatype = gdal_dtypes[im_data.dtype.name]
        else:
            datatype = gdal.GDT_Float32

        # 判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape

        # 创建文件
        if os.path.exists(os.path.split(filename)[0]) is False:
            os.makedirs(os.path.split(filename)[0])

        driver = gdal.GetDriverByName("ENVI")  # 数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

        if im_bands == 1:
            # outRaster.GetRasterBand(1).WriteArray(array)  # 写入数组数据
            outband = dataset.GetRasterBand(1)
            outband.WriteArray(im_data)
            if no_data != 'null':
                outband.SetNoDataValue(no_data)
            outband.FlushCache()
        else:
            for i in range(im_bands):
                outband = dataset.GetRasterBand(1 + i)
                outband.WriteArray(im_data[i])
                outband.FlushCache()
                # outRaster.GetRasterBand(i + 1).WriteArray(array[i])
        del dataset
        dataset=None

    @staticmethod
    def get_data(filename):
        """
        :param filename: tif路径
        :return: 获取所有波段的数据
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        im_width = dataset.RasterXSize
        im_height = dataset.RasterYSize
        im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
        del dataset
        dataset=None
        return im_data

    @staticmethod
    def get_scopes(ori_sim):
        ori_sim_data = ImageHandler.get_data(ori_sim)
        lon = ori_sim_data[0, :, :]
        lat = ori_sim_data[1, :, :]

        min_lon = np.nanmin(lon)
        max_lon = np.nanmax(lon)
        min_lat = np.nanmin(lat)
        max_lat = np.nanmax(lat)

        scopes = [[min_lon, max_lat], [max_lon, max_lat], [min_lon, min_lat], [max_lon, min_lat]]
        return scopes

    @staticmethod
    def get_dataset(filename):
        """
        :param filename: tif路径
        :return: 获取所有波段的数据
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        return dataset

    @staticmethod
    def get_all_band_array(filename):
        """
        （大气延迟算法）
        将ERA-5影像所有波段存为一个数组, 波段数在第三维度           get_data（）->（37，8，8）
        :param filename： 影像路径                       get_all_band_array ->（8，8，37）
        :return: 影像数组
        """
        dataset = gdal.Open(filename)
        x_size = dataset.RasterXSize
        y_size = dataset.RasterYSize
        nums = dataset.RasterCount
        array = np.zeros((y_size, x_size, nums), dtype=float)
        if nums == 1:
            bands_0 = dataset.GetRasterBand(1)
            array = bands_0.ReadAsArray(0, 0, x_size, y_size)
        else:
            for i in range(0, nums):
                bands = dataset.GetRasterBand(i + 1)
                arr = bands.ReadAsArray(0, 0, x_size, y_size)
                array[:, :, i] = arr
        return array

    @staticmethod
    def get_img_width(filename):
        """
        :param filename: tif路径
        :return: 影像宽度
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        width = dataset.RasterXSize

        del dataset
        dataset=None
        return width

    @staticmethod
    def get_img_height(filename):
        """
        :param filename: tif路径
        :return: 影像高度
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None
        height = dataset.RasterYSize
        del dataset
        dataset=None
        return height

    @staticmethod
    def read_img(filename):
        """
        影像读取
        :param filename:
        :return:
        """
        gdal.AllRegister()
        img_dataset = gdal.Open(filename)  # 打开文件

        if img_dataset is None:
            msg = 'Could not open ' + filename
            print("ERROR",msg)
            return None, None, None

        im_proj = img_dataset.GetProjection()  # 地图投影信息
        if im_proj is None:
            return None, None, None
        im_geotrans = img_dataset.GetGeoTransform()  # 仿射矩阵

        im_width = img_dataset.RasterXSize  # 栅格矩阵的行数
        im_height = img_dataset.RasterYSize  # 栅格矩阵的行数
        im_arr = img_dataset.ReadAsArray(0, 0, im_width, im_height)
        del img_dataset
        return im_proj, im_geotrans, im_arr

    def cal_img_scope(self, dataset):
        """
        计算影像的地理坐标范围
        根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
        :param dataset :GDAL地理数据
        :return: list[point_upleft, point_upright, point_downleft, point_downright]
        """
        if dataset is None:
            return None

        img_geotrans = dataset.GetGeoTransform()
        if img_geotrans is None:
            return None

        width = dataset.RasterXSize  # 栅格矩阵的列数
        height = dataset.RasterYSize  # 栅格矩阵的行数

        point_upleft = self.trans_rowcol2geo(img_geotrans, 0, 0)
        point_upright = self.trans_rowcol2geo(img_geotrans, width, 0)
        point_downleft = self.trans_rowcol2geo(img_geotrans, 0, height)
        point_downright = self.trans_rowcol2geo(img_geotrans, width, height)

        return [point_upleft, point_upright, point_downleft, point_downright]

    @staticmethod
    def get_scope_ori_sim(filename):
        """
        计算影像的地理坐标范围
        根据GDAL的六参数模型将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
        :param dataset :GDAL地理数据
        :return: list[point_upleft, point_upright, point_downleft, point_downright]
        """
        gdal.AllRegister()
        dataset = gdal.Open(filename)
        if dataset is None:
            return None

        width = dataset.RasterXSize  # 栅格矩阵的列数
        height = dataset.RasterYSize  # 栅格矩阵的行数

        band1 = dataset.GetRasterBand(1)
        array1 = band1.ReadAsArray(0, 0, band1.XSize, band1.YSize)

        band2 = dataset.GetRasterBand(2)
        array2 = band2.ReadAsArray(0, 0, band2.XSize, band2.YSize)

        if array1[0, 0] < array1[0, width - 1]:
            point_upleft = [array1[0, 0], array2[0, 0]]
            point_upright = [array1[0, width - 1], array2[0, width - 1]]
        else:
            point_upright = [array1[0, 0], array2[0, 0]]
            point_upleft = [array1[0, width - 1], array2[0, width - 1]]

        if array1[height - 1, 0] < array1[height - 1, width - 1]:
            point_downleft = [array1[height - 1, 0], array2[height - 1, 0]]
            point_downright = [array1[height - 1, width - 1], array2[height - 1, width - 1]]
        else:
            point_downright = [array1[height - 1, 0], array2[height - 1, 0]]
            point_downleft = [array1[height - 1, width - 1], array2[height - 1, width - 1]]

        if (array2[0, 0] < array2[height - 1, 0]):
            # 上下调换顺序
            tmp1 = point_upleft
            point_upleft = point_downleft
            point_downleft = tmp1

            tmp2 = point_upright
            point_upright = point_downright
            point_downright = tmp2

        return [point_upleft, point_upright, point_downleft, point_downright]

    @staticmethod
    def trans_rowcol2geo(img_geotrans, img_col, img_row):
        """
        据GDAL的六参数模型仿射矩阵将影像图上坐标（行列号）转为投影坐标或地理坐标（根据具体数据的坐标系统转换）
        :param img_geotrans: 仿射矩阵
        :param img_col:图像纵坐标
        :param img_row:图像横坐标
        :return: [geo_x,geo_y]
        """
        geo_x = img_geotrans[0] + img_geotrans[1] * img_col + img_geotrans[2] * img_row
        geo_y = img_geotrans[3] + img_geotrans[4] * img_col + img_geotrans[5] * img_row
        return [geo_x, geo_y]

    @staticmethod
    def write_era_into_img(filename, im_proj, im_geotrans, im_data):
        """
        影像保存
        :param filename:
        :param im_proj:
        :param im_geotrans:
        :param im_data:
        :return:
        """
        gdal_dtypes = {
            'int8': gdal.GDT_Byte,
            'unit16': gdal.GDT_UInt16,
            'int16': gdal.GDT_Int16,
            'unit32': gdal.GDT_UInt32,
            'int32': gdal.GDT_Int32,
            'float32': gdal.GDT_Float32,
            'float64': gdal.GDT_Float64,
        }
        if not gdal_dtypes.get(im_data.dtype.name, None) is None:
            datatype = gdal_dtypes[im_data.dtype.name]
        else:
            datatype = gdal.GDT_Float32

        # 判读数组维数
        if len(im_data.shape) == 3:
            im_height, im_width, im_bands = im_data.shape  # shape[0] 行数
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape

        # 创建文件
        if os.path.exists(os.path.split(filename)[0]) is False:
            os.makedirs(os.path.split(filename)[0])

        driver = gdal.GetDriverByName("ENVI")  # 数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影

        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i + 1).WriteArray(im_data[:, :, i])
                # dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
        del dataset
        dataset=None

        # 写GeoTiff文件

    @staticmethod
    def lat_lon_to_pixel(raster_dataset_path, location):
        """From zacharybears.com/using-python-to-translate-latlon-locations-to-pixels-on-a-geotiff/."""
        gdal.AllRegister()
        raster_dataset = gdal.Open(raster_dataset_path)
        if raster_dataset is None:
            return None
        ds = raster_dataset
        gt = ds.GetGeoTransform()
        srs = osr.SpatialReference()
        srs.ImportFromWkt(ds.GetProjection())
        srs_lat_lon = srs.CloneGeogCS()
        ct = osr.CoordinateTransformation(srs_lat_lon, srs)
        new_location = [None, None]
        # Change the point locations into the GeoTransform space
        (new_location[1], new_location[0], holder) = ct.TransformPoint(location[1], location[0])
        # Translate the x and y coordinates into pixel values
        Xp = new_location[0]
        Yp = new_location[1]
        dGeoTrans = gt
        dTemp = dGeoTrans[1] * dGeoTrans[5] - dGeoTrans[2] * dGeoTrans[4]
        Xpixel = (dGeoTrans[5] * (Xp - dGeoTrans[0]) - dGeoTrans[2] * (Yp - dGeoTrans[3])) / dTemp
        Yline = (dGeoTrans[1] * (Yp - dGeoTrans[3]) - dGeoTrans[4] * (Xp - dGeoTrans[0])) / dTemp
        del raster_dataset
        return (Xpixel, Yline)

    @staticmethod
    def write_img(filename, im_proj, im_geotrans, im_data, no_data='0'):
        """
        影像保存
        :param filename: 保存的路径
        :param im_proj:
        :param im_geotrans:
        :param im_data:
        :param no_data: 把无效值设置为 nodata
        :return:
        """

        gdal_dtypes = {
            'int8': gdal.GDT_Byte,
            'unit16': gdal.GDT_UInt16,
            'int16': gdal.GDT_Int16,
            'unit32': gdal.GDT_UInt32,
            'int32': gdal.GDT_Int32,
            'float32': gdal.GDT_Float32,
            'float64': gdal.GDT_Float64,
        }
        if not gdal_dtypes.get(im_data.dtype.name, None) is None:
            datatype = gdal_dtypes[im_data.dtype.name]
        else:
            datatype = gdal.GDT_Float32
        flag = False
        # 判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
            flag = True
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape

        # 创建文件
        if os.path.exists(os.path.split(filename)[0]) is False:
            os.makedirs(os.path.split(filename)[0])

        driver = gdal.GetDriverByName("ENVI")  # 数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数

        dataset.SetProjection(im_proj)  # 写入投影

        if im_bands == 1:
            # outRaster.GetRasterBand(1).WriteArray(array)  # 写入数组数据
            if flag:
                outband = dataset.GetRasterBand(1)
                outband.WriteArray(im_data[0])
                if no_data != 'null':
                    outband.SetNoDataValue(np.double(no_data))
                outband.FlushCache()
            else:
                outband = dataset.GetRasterBand(1)
                outband.WriteArray(im_data)
                if no_data != 'null':
                    outband.SetNoDataValue(np.double(no_data))
                outband.FlushCache()
        else:
            for i in range(im_bands):
                outband = dataset.GetRasterBand(1 + i)
                outband.WriteArray(im_data[i])
                if no_data != 'null':
                    outband.SetNoDataValue(np.double(no_data))
                outband.FlushCache()
                # outRaster.GetRasterBand(i + 1).WriteArray(array[i])
        del dataset
        dataset=None

        # 写GeoTiff文件

    @staticmethod
    def write_img_envi(filename, im_proj, im_geotrans, im_data, no_data='null'):
        """
        影像保存
        :param filename: 保存的路径
        :param im_proj:
        :param im_geotrans:
        :param im_data:
        :param no_data: 把无效值设置为 nodata
        :return:
        """

        gdal_dtypes = {
            'int8': gdal.GDT_Byte,
            'unit16': gdal.GDT_UInt16,
            'int16': gdal.GDT_Int16,
            'unit32': gdal.GDT_UInt32,
            'int32': gdal.GDT_Int32,
            'float32': gdal.GDT_Float32,
            'float64': gdal.GDT_Float64,
        }
        if not gdal_dtypes.get(im_data.dtype.name, None) is None:
            datatype = gdal_dtypes[im_data.dtype.name]
        else:
            datatype = gdal.GDT_Float32

        # 判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape

        # 创建文件
        if os.path.exists(os.path.split(filename)[0]) is False:
            os.makedirs(os.path.split(filename)[0])

        driver = gdal.GetDriverByName("ENVI")  # 数据类型必须有，因为要计算需要多大内存空间
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数

        dataset.SetProjection(im_proj)  # 写入投影

        if im_bands == 1:
            # outRaster.GetRasterBand(1).WriteArray(array)  # 写入数组数据
            outband = dataset.GetRasterBand(1)
            outband.WriteArray(im_data)
            if no_data != 'null':
                outband.SetNoDataValue(no_data)
            outband.FlushCache()
        else:
            for i in range(im_bands):
                outband = dataset.GetRasterBand(1 + i)
                outband.WriteArray(im_data[i])
                outband.FlushCache()
                # outRaster.GetRasterBand(i + 1).WriteArray(array[i])
        del dataset
        dataset=None

    @staticmethod
    def write_img_rpc(filename, im_proj, im_geotrans, im_data, rpc_dict):
        """
        图像中写入rpc信息
        """
        # 判断栅格数据的数据类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_Int16
        else:
            datatype = gdal.GDT_Float32

        # 判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            im_bands, (im_height, im_width) = 1, im_data.shape

        # 创建文件
        driver = gdal.GetDriverByName("ENVI")
        dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)

        dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数
        dataset.SetProjection(im_proj)  # 写入投影

        # 写入RPC参数
        for k in rpc_dict.keys():
            dataset.SetMetadataItem(k, rpc_dict[k], 'RPC')

        if im_bands == 1:
            dataset.GetRasterBand(1).WriteArray(im_data)  # 写入数组数据
        else:
            for i in range(im_bands):
                dataset.GetRasterBand(i + 1).WriteArray(im_data[i])

        del dataset
        dataset=None

    def transtif2mask(self, out_tif_path, in_tif_path, threshold):
        """
        :param out_tif_path:输出路径
        :param in_tif_path:输入的路径
        :param threshold:阈值
        """
        im_proj, im_geotrans, im_arr, im_scope = self.read_img(in_tif_path)
        im_arr_mask = (im_arr < threshold).astype(int)
        self.write_img(out_tif_path, im_proj, im_geotrans, im_arr_mask)


    def get_center_scopes(self, dataset):
        if dataset is None:
            return None

        img_geotrans = dataset.GetGeoTransform()
        if img_geotrans is None:
            return None

        width = dataset.RasterXSize  # 栅格矩阵的列数
        height = dataset.RasterYSize  # 栅格矩阵的行数

        x_split = int(width / 5)
        y_split = int(height / 5)
        img_col_start = x_split * 1
        img_col_end = x_split * 3
        img_row_start = y_split * 1
        img_row_end = y_split * 3
        cols = img_col_end - img_col_start
        rows = img_row_end - img_row_start
        if cols > 10000 or rows > 10000:
            img_col_end = img_col_start + 10000
            img_row_end = img_row_start + 10000

        point_upleft = self.trans_rowcol2geo(img_geotrans, img_col_start, img_row_start)
        point_upright = self.trans_rowcol2geo(img_geotrans, img_col_end, img_row_start)
        point_downleft = self.trans_rowcol2geo(img_geotrans, img_col_start, img_row_end)
        point_downright = self.trans_rowcol2geo(img_geotrans, img_col_end, img_row_end)

        return [point_upleft, point_upright, point_downleft, point_downright]

    def limit_field(self, out_path, in_path, min_value, max_value):
        """
        :param out_path:输出路径
        :param in_path:主mask路径，输出影像采用主mask的地理信息
        :param min_value
        :param max_value
        """
        proj = self.get_projection(in_path)
        geotrans = self.get_geotransform(in_path)
        array = self.get_band_array(in_path, 1)
        array[array < min_value] = min_value
        array[array > max_value] = max_value
        self.write_img(out_path, proj, geotrans, array)
        return True

    def band_merge(self, lon, lat, ori_sim):
        lon_arr = self.get_data(lon)
        lat_arr = self.get_data(lat)
        temp = np.zeros((2, lon_arr.shape[0], lon_arr.shape[1]), dtype=float)
        temp[0, :, :] = lon_arr[:, :]
        temp[1, :, :] = lat_arr[:, :]
        self.write_img(ori_sim, '', [0.0, 1.0, 0.0, 0.0, 0.0, 1.0], temp, '0')

    def get_scopes(self, ori_sim):
        ori_sim_data = self.get_data(ori_sim)
        lon = ori_sim_data[0, :, :]
        lat = ori_sim_data[1, :, :]

        min_lon = np.nanmin(np.where((lon != 0) & ~np.isnan(lon), lon, np.inf))
        max_lon = np.nanmax(np.where((lon != 0) & ~np.isnan(lon), lon, -np.inf))
        min_lat = np.nanmin(np.where((lat != 0) & ~np.isnan(lat), lat, np.inf))
        max_lat = np.nanmax(np.where((lat != 0) & ~np.isnan(lat), lat, -np.inf))

        scopes = [[min_lon, max_lat], [max_lon, max_lat], [min_lon, min_lat], [max_lon, min_lat]]
        return scopes

    @staticmethod
    def dem_merged(in_dem_path, out_dem_path):
        '''
        DEM重采样函数，默认坐标系为WGS84
        agrs:
            in_dem_path: 输入的DEM文件夹路径
            meta_file_path: 输入的xml元文件路径
            out_dem_path: 输出的DEM文件夹路径
        '''
        # 读取文件夹中所有的DEM
        dem_file_paths = [os.path.join(in_dem_path, dem_name) for dem_name in os.listdir(in_dem_path) if
                          dem_name.find(".tif") >= 0 and dem_name.find(".tif.") == -1]
        spatialreference = osr.SpatialReference()
        spatialreference.SetWellKnownGeogCS("WGS84")  # 设置地理坐标，单位为度 degree # 设置投影坐标，单位为度 degree
        spatialproj = spatialreference.ExportToWkt()  # 导出投影结果
        # 将DEM拼接成一张大图
        mergeFile = gdal.BuildVRT(os.path.join(out_dem_path, "mergedDEM_VRT.tif"), dem_file_paths)
        out_DEM = os.path.join(out_dem_path, "mergedDEM.tif")
        gdal.Warp(out_DEM,
                  mergeFile,
                  format="ENVI",
                  dstSRS=spatialproj,
                  dstNodata=-9999,
                  outputType=gdal.GDT_Float32)
        time.sleep(3)
        # gdal.CloseDir(out_DEM)
        return out_DEM

    @staticmethod
    def landcover_merged(in_dem_path, out_dem_path):
        '''
        DEM重采样函数，默认坐标系为WGS84
        agrs:
            in_dem_path: 输入的DEM文件夹路径
            meta_file_path: 输入的xml元文件路径
            out_dem_path: 输出的DEM文件夹路径
        '''
        # 读取文件夹中所有的DEM
        dem_file_paths = [os.path.join(in_dem_path, dem_name) for dem_name in os.listdir(in_dem_path) if
                          dem_name.find(".tif") >= 0 and dem_name.find(".tif.") == -1]
        spatialreference = osr.SpatialReference()
        spatialreference.SetWellKnownGeogCS("WGS84")  # 设置地理坐标，单位为度 degree # 设置投影坐标，单位为度 degree
        spatialproj = spatialreference.ExportToWkt()  # 导出投影结果
        # 将DEM拼接成一张大图
        mergeFile = gdal.BuildVRT(os.path.join(out_dem_path, "mergedDEM_VRT.tif"), dem_file_paths)
        out_DEM = os.path.join(out_dem_path, "mergedDEM.tif")
        gdal.Warp(out_DEM,
                  mergeFile,
                  format="ENVI",
                  dstSRS=spatialproj,
                  dstNodata=-9999,
                  outputType=gdal.GDT_Byte)
        time.sleep(3)
        # gdal.CloseDir(out_DEM)
        return out_DEM

    @staticmethod
    def get_inc_angle(inc_xml, rows, cols, out_path):
        tree = ElementTree()
        tree.parse(inc_xml)  # 影像头文件
        root = tree.getroot()
        values = root.findall('incidenceValue')
        angle_value = [value.text for value in values]
        angle_value = np.array(angle_value)
        inc_angle = np.tile(angle_value, (rows, 1))
        ImageHandler.write_img(out_path, '', [0.0, 1.0, 0.0, 0.0, 0.0, 1.0], inc_angle)

        pass


if __name__ == '__main__':
    cols = 7086
    rows = 8064
    inc_xml = r"D:\micro\WorkSpace\Dem\Temporary\processing\product\GF3_SAY_FSI_001614_E113.2_N34.5_20161129_L1A_HHHV_L10002015686-DEM.tiff"
    # ImageHandler().write_quick_view(inc_xml)
    # fn = r'E:\202306hb\result\GF3B_SYC_QPSI_008316_E116.1_N43.3_20230622_L1A_AHV_L10000202892-cal-SMC.tif'
    # out = r'E:\202306hb\result\soil.tif'
    # #
    # im_proj, im_geotrans, im_arr = ImageHandler.read_img(fn)
    # im_arr[np.where(np.isnan(im_arr))] = 0
    # # h,w = im_arr.shape
    # # arr = np.random.rand(h,w)*0.4
    # ImageHandler.write_img(out, im_proj, im_geotrans, im_arr, '-1')

    #     path = r'D:\BaiduNetdiskDownload\GZ\lon.rdr'
    #     path2 = r'D:\BaiduNetdiskDownload\GZ\lat.rdr'
    #     path3 = r'D:\BaiduNetdiskDownload\GZ\lon_lat.tif'
    #     s = ImageHandler().band_merge(path, path2, path3)
    #     print(s)
    #     pass
    fileDir = r'D:\micro\WorkSpace\ortho\Temporary\VH_db'
    outDir = r'D:\micro\WorkSpace\ortho\Temporary\VH_db\test'
    files = os.listdir(fileDir)
    for file in files:
        tifFile = os.path.join(fileDir, file)
        outFile = os.path.join(outDir, file)
        im_proj, im_geotrans, im_arr = ImageHandler.read_img(tifFile)
        im_arr[np.isnan(im_arr)] = 0
        im_arr[np.isinf(im_arr)] = 0
        ImageHandler.write_img(outFile, im_proj, im_geotrans, im_arr, '0')
