# -*- coding: UTF-8 -*- """ @Project ：microproduct @File ：CalculateIncident.py @Function ：计算、局部入射角计算 @Author ：LMM @Date ：2021/8/25 14:17 @Version ：1.0.0 """ import os import numpy as np from osgeo import gdal from osgeo import gdalconst import gc import math from xml.dom import minidom # 不需要安装，默认环境里就有 class CalculateIncident: def __init__(self): pass @staticmethod def add_round(npgrid): """ 边缘填充一圈,然后输出填充得到的矩阵 param：npgrid dem数组 """ ny, nx = npgrid.shape # ny:行数，nx:列数 zbc = np.zeros((ny + 2, nx + 2)) zbc[1:-1, 1:-1] = npgrid # 四边 zbc[0, 1:-1] = npgrid[0, :] zbc[-1, 1:-1] = npgrid[-1, :] zbc[1:-1, 0] = npgrid[:, 0] zbc[1:-1, -1] = npgrid[:, -1] # 角点 zbc[0, 0] = npgrid[0, 0] zbc[0, -1] = npgrid[0, -1] zbc[-1, 0] = npgrid[-1, 0] zbc[-1, -1] = npgrid[-1, -1] print("输出填充后的数组的形状", zbc.shape) return zbc @staticmethod def cal_dxdy(zbc, dx): """ 计算dx,dy param：zbc填充后的数组 param：dx dem数据像元大小 """ we_x = ((zbc[1:-1, :-2]) - (zbc[1:-1, 2:])) / dx / 2 # WE方向 ns_y = ((zbc[2:, 1:-1]) - (zbc[:-2, 1:-1])) / dx / 2 # NS方向 print("输出Sx的数组的形状", we_x.shape, "输出Sy的数组的形状", ns_y.shape) sx = we_x[1:-1, 1:-1] sy = ns_y[1:-1, 1:-1] # np.savetxt("dxdy.csv",dx,delimiter=",") print("输出Sx2的数组的形状", sx.shape, "输出Sy2的数组的形状", sy.shape) return sx, sy @staticmethod def cal_slopasp(dx, dy): # 计算坡度\坡向 # 坡度计算 slope slope = (np.arctan(np.sqrt(dx * dx + dy * dy))) * 57.29578 # 转换成°,57.29578=180/math.pi slope = slope[1:-1, 1:-1] # 坡向计算 aspect aspect = np.ones([dx.shape[0], dx.shape[1]]).astype(np.float32) # 生成一个全是0的数组 # dx = dx.astype(np.float32) # dy = dy.astype(np.float32) # a1=(np.where(dx==0) and np.where(dy ==0)) # print(a1) # aspect[a1]=-1 # a2 = (np.where(dx == 0) and np.where(dy > 0)) # aspect[a2] =0.0 # a3 = (np.where(dx == 0) and np.where(dy <0)) # aspect[a3] =180.0 # a4 = (np.where(dx > 0) and np.where(dy ==0)) # aspect[a4] =90.0 # a5 = (np.where(dx < 0) and np.where(dy ==0)) # aspect[a5] =270.0 # a6 = (np.where(dx != 0) or np.where(dy !=0)) # b=dy[a6] # print(":", 1) # aspect[a6] =float(math.atan2(dy[i, j], dx[i, j])) * 57.29578 # a7=np.where(aspect[a6]< 0.0) # aspect[a7] = 90.0 - aspect[a7] # a8=np.where(aspect[a6]> 90.0) # aspect[a8] = 450.0- aspect[a8] # a9 =np.where(aspect[a6] >= 0 or aspect[a6] <= 90) # aspect[a9] =90.0 - aspect[a9] for i in range(dx.shape[0]): for j in range(dx.shape[1]): x = float(dx[i, j]) y = float(dy[i, j]) if (x == 0.0) & (y == 0.0): aspect[i, j] = -1 elif x == 0.0: if y > 0.0: aspect[i, j] = 0.0 else: aspect[i, j] = 180.0 elif y == 0.0: if x > 0.0: aspect[i, j] = 90.0 else: aspect[i, j] = 270.0 else: aspect[i, j] = float(math.atan2(y, x)) * 57.29578 # 范围（-Π/2，Π/2） if aspect[i, j] < 0.0: aspect[i, j] = 90.0 - aspect[i, j] elif aspect[i, j] > 90.0: aspect[i, j] = 450.0 - aspect[i, j] else: aspect[i, j] = 90.0 - aspect[i, j] print("输出aspect形状:", aspect.shape) # 3599, 3599 print("输出aspect:", aspect) return slope, aspect def creat_twofile(self, dem_file_path, slope_out_path, aspect_out_path): """ 生成坡度图、坡向图 param: path_file1 为输入文件tif数据的文件路径 """ if os.path.isfile(dem_file_path): print("高程数据文件存在") else: print("高程数据文件不存在") dataset_caijian = gdal.Open(dem_file_path) x_size = dataset_caijian.RasterXSize y_size = dataset_caijian.RasterYSize geo = dataset_caijian.GetGeoTransform() pro = dataset_caijian.GetProjection() array0 = dataset_caijian.ReadAsArray(0, 0, x_size, y_size) print("输出dem数据的数组", array0) zbc = self.add_round(array0) sx, sy = self.cal_dxdy(zbc, 30) slope, aspect = self.cal_slopasp(sx, sy) driver = gdal.GetDriverByName("GTiff") # 创建一个数据格式 driver.Register() newfile = driver.Create(slope_out_path, x_size, y_size, 1, gdal.GDT_Float32) # 存放路径文件名，长，宽，波段，数据类型 newfile.SetProjection(pro) geo = [geo[0], geo[1], 0, geo[3], 0, -geo[1]] newfile.SetGeoTransform(geo) newfile.GetRasterBand(1).WriteArray(slope) driver2 = gdal.GetDriverByName("GTiff") # 创建一个数据格式 driver2.Register() newfile2 = driver2.Create(aspect_out_path, x_size, y_size, 1, gdal.GDT_Float32) # 存放路径文件名，长，宽，波段，数据类型 geo = [geo[0], geo[1], 0, geo[3], 0, -geo[1]] newfile2.SetGeoTransform(geo) newfile2.GetRasterBand(1).WriteArray(aspect) @staticmethod def resampling(input_file1, input_file2, ref_file, output_file, output_file2): """ 采用gdal.Warp()方法进行重采样，差值法为双线性插值 :param input_file1 slope path :param input_file2 aspect path :param ref_file: 参考图像路径 :param output_file: slope path :param output_file2 aspect path :return: """ gdal.AllRegister() in_ds1 = gdal.Open(input_file1) in_ds2 = gdal.Open(input_file2) ref_ds = gdal.Open(ref_file, gdal.GA_ReadOnly) # 获取输入影像信息 input_file_proj = in_ds1.GetProjection() # inputefileTrans = in_ds1.GetGeoTransform() reference_file_proj = ref_ds.GetProjection() reference_file_trans = ref_ds.GetGeoTransform() nbands = in_ds1.RasterCount bandinputfile1 = in_ds1.GetRasterBand(1) bandinputfile2 = in_ds2.GetRasterBand(1) x = ref_ds.RasterXSize y = ref_ds.RasterYSize # 创建重采样输出文件（设置投影及六参数） driver1 = gdal.GetDriverByName('GTiff') output1 = driver1.Create(output_file, x, y, nbands, bandinputfile1.DataType) output1.SetGeoTransform(reference_file_trans) output1.SetProjection(reference_file_proj) # options = gdal.WarpOptions(srcSRS=inputProj, dstSRS=referencefileProj, resampleAlg=gdalconst.GRA_Bilinear) # resampleAlg = gdalconst.GRA_NearestNeighbour gdal.ReprojectImage(in_ds1, output1, input_file_proj, reference_file_proj, gdalconst.GRA_Bilinear) driver2 = gdal.GetDriverByName('GTiff') output2 = driver2.Create(output_file2, x, y, nbands, bandinputfile2.DataType) output2.SetGeoTransform(reference_file_trans) output2.SetProjection(reference_file_proj) # options = gdal.WarpOptions(srcSRS=inputProj, dstSRS=referencefileProj, resampleAlg=gdalconst.GRA_Bilinear) # resampleAlg = gdalconst.GRA_NearestNeighbour gdal.ReprojectImage(in_ds2, output2, input_file_proj, reference_file_proj, gdalconst.GRA_Bilinear) @staticmethod def getorbitparameter(xml_path): """ 从轨道参数文件xml中获取升降轨信息、影像四个角的经纬度坐标 """ # 打开xml文档,根据路径初始化DOM doc = minidom.parse(xml_path) # 得到xml文档元素对象,初始化root对象 root = doc.documentElement # 输出升降轨信息，DEC降轨，ASC升轨 direction = root.getElementsByTagName("Direction")[0] # print("输出Direction的子节点列表",Direction.firstChild.data) pd = direction.firstChild.data imageinfo = root.getElementsByTagName("imageinfo")[0] # 输出topLeft的纬度和经度 top_left = imageinfo.getElementsByTagName("topLeft")[0] latitude = top_left.getElementsByTagName("latitude")[0] longitude = top_left.getElementsByTagName("longitude")[0] # print("输出topLeft的纬度lat和经度lon：", latitude.firstChild.data,longitude.firstChild.data) tl_lat, tl_lon = latitude.firstChild.data, longitude.firstChild.data # 输出topRight的纬度和经度 top_right = imageinfo.getElementsByTagName("topRight")[0] latitude = top_right.getElementsByTagName("latitude")[0] longitude = top_right.getElementsByTagName("longitude")[0] # print("输出topRight的纬度lat和经度lon：", latitude.firstChild.data,longitude.firstChild.data) tr_lat, tr_lon = latitude.firstChild.data, longitude.firstChild.data # 输出 bottomLeft的纬度和经度 bottom_left = imageinfo.getElementsByTagName("bottomLeft")[0] latitude = bottom_left.getElementsByTagName("latitude")[0] longitude = bottom_left.getElementsByTagName("longitude")[0] # print("输出bottomLeft的纬度lat和经度lon：", latitude.firstChild.data,longitude.firstChild.data) bl_lat, bl_lon = latitude.firstChild.data, longitude.firstChild.data # 输出 bottomRight的纬度和经度 bottom_right = imageinfo.getElementsByTagName("bottomRight")[0] latitude = bottom_right.getElementsByTagName("latitude")[0] longitude = bottom_right.getElementsByTagName("longitude")[0] # print("输出bottomRight的纬度lat和经度lon：", latitude.firstChild.data,longitude.firstChild.data) br_lat, br_lon = latitude.firstChild.data, longitude.firstChild.data print("pd,tl_lat,tl_lon,tr_lat,tr_lon,bl_lat,bl_lon,br_lat,br_lon", pd, tl_lat, tl_lon, tr_lat, tr_lon, bl_lat, bl_lon, br_lat, br_lon) return pd, tl_lat, tl_lon, tr_lat, tr_lon, bl_lat, bl_lon, br_lat, br_lon def get_rparademeter(self, xml_path): """ 计算雷达视线向方向角R """ pd, tl_lat, tl_lon, tr_lat, tr_lon, bl_lat, bl_lon, br_lat, br_lon = self.getorbitparameter(xml_path) tl_lat = float(tl_lat) # 原来的数是带有小数点的字符串，int会报错，使用float tl_lon = float(tl_lon) # tr_lat = float(tr_lat) # tr_lon = float(tr_lon) bl_lat = float(bl_lat) bl_lon = float(bl_lon) # br_lat = float(br_lat) # br_lon = float(br_lon) if pd == "DEC": # 降轨 b = np.arctan((tl_lat - bl_lat) / (tl_lon - bl_lon)) * 57.29578 r = 270 + b return r # tl_lat, tl_lon = lonlat2geo(tl_lat, tl_lon) # tr_lat, tr_lon = lonlat2geo(tr_lat, tr_lon) # bl_lat, bl_lon = lonlat2geo(bl_lat, bl_lon) # br_lat, br_lon = lonlat2geo(br_lat, br_lon) # B2 = np.arctan((tl_lat - bl_lat) / (tl_lon - bl_lon)) * 57.29578 # R2 = 270 + B2 # print(("输出R2：", R2)) if pd == "ASC": # 升轨 b = np.arctan((tl_lat - bl_lat) / (tl_lon - bl_lon)) * 57.29578 return b def clau(self, pathfile1, pathfile2, pathfile3, xml_path, save_localangle_path): """ 计算局部入射角 param: pathfile1是slope的坡度图路径 param: pathfile2是aspect的坡向图路径 param: pathfile3是入射角文件的路径 param: xml_path是轨道参数文件 r是雷达视线向方位角 """ r = self.get_rparademeter(xml_path) pd, tl_lat, tl_lon, tr_lat, tr_lon, bl_lat, bl_lon, br_lat, br_lon = self.getorbitparameter(xml_path) print("输出升降轨：", pd) dataset = gdal.Open(pathfile1) x = dataset.RasterXSize y = dataset.RasterYSize print("输出slope的行、列：", x, y) slope_array = dataset.ReadAsArray(0, 0, x, y) dataset2 = gdal.Open(pathfile2) x2 = dataset2.RasterXSize y2 = dataset2.RasterYSize print("输出aspect的行、列：", x2, y2) aspect_array = dataset2.ReadAsArray(0, 0, x2, y2) dataset3 = gdal.Open(pathfile3) x3 = dataset3.RasterXSize y3 = dataset3.RasterYSize geo3 = dataset3.GetGeoTransform() pro3 = dataset3.GetProjection() print("输出入射角文件的行、列：", x3, y3) rushe_array = dataset3.ReadAsArray(0, 0, x3, y3) # b0 = np.where(rushe_array > 0.00001, 0, 1) radina_value = 0 if pd == "DEC": # 降轨数据 # 雷达视线角-坡度角在90度到270度之间 where_0 = np.where(rushe_array == 0) bb1 = (r-aspect_array).all() and (r-aspect_array).all() bb2 = np.where(90 < bb1 < 270, 1, 0) b1 = (bb1 and bb2) # b1 = np.where(90 < ((r-aspect_array).all()) and ((r-aspect_array).all()) < 270, 1, 0) c1 = np.cos(rushe_array*(math.pi/180)) * np.cos(slope_array*(math.pi/180)) - np.sin(slope_array*(math.pi/180)) * np.sin( rushe_array*(math.pi/180)) * np.cos((r - aspect_array)*(math.pi/180)) d1 = b1 * c1 # 雷达视线角-坡度角=90度或=270度时 b2 = np.where((r-aspect_array == 90) | (r-aspect_array == 270), 1, 0) d2 = b2*c1 # 雷达视线角-坡度角在90度到270度之间 b3 = 1-b1-b2 c3 = np.cos(rushe_array*(math.pi/180)) * np.cos(slope_array*(math.pi/180)) + np.sin( slope_array*(math.pi/180)) * np.sin(rushe_array*(math.pi/180)) * np.cos((r - aspect_array)*(math.pi/180)) d3 = b3 * c3 del b1, b2, b3, c3, c1 gc.collect() radina_value = d1 + d2 + d3 radina_value[where_0] = 0 del d1, d2, d3 gc.collect() if pd == "ASC": # 升轨数据 # 坡度-雷达视线角在90度到270度之间 where_0 = np.where(rushe_array == 0) bb1 = (r-aspect_array).all() and (r-aspect_array).all() bb2 = np.where(90 < bb1 < 270, 1, 0) b1 = (bb1 and bb2) # b1 = np.where(90 < ((r-aspect_array).all()) and ((r-aspect_array).all()) < 270, 1, 0) c1 = np.cos(rushe_array*(math.pi/180)) * np.cos(slope_array*(math.pi/180)) + np.sin( slope_array*(math.pi/180)) * np.sin(rushe_array*(math.pi/180)) * np.cos((r - aspect_array)*(math.pi/180)) d1 = b1 * c1 # 坡度-雷达视线角=90或=270时 b2 = np.where((aspect_array-r == 90) | (aspect_array-r == 270), 1, 0) d2 = b2 * c1 # 坡度-雷达视线角在0-90度或270-360度之间 b3 = 1 - b1-b2 c3 = np.cos(rushe_array*(math.pi/180)) * np.cos(slope_array*(math.pi/180)) - np.sin(slope_array*(math.pi/180)) *\ np.sin(rushe_array*(math.pi/180)) * np.cos((r - aspect_array)*(math.pi/180)) d3 = b3 * c3 radina_value = d1 + d2 + d3 radina_value[where_0] = 0 del b1, b2, b3, c3, c1, d1, d2, d3 gc.collect() jubu_o = 57.29578 * np.arccos(radina_value) print("输出局部入射角", jubu_o) driver = gdal.GetDriverByName("GTiff") # 创建一个数据格式 driver.Register() newfile = driver.Create(save_localangle_path, x3, y3, 1, gdal.GDT_Float32) # 存放路径文件名，长，宽，波段，数据类型 newfile.SetProjection(pro3) newfile.SetGeoTransform(geo3) newfile.GetRasterBand(1).WriteArray(jubu_o) def localangle(self, dem_path, incidence_angle_path, orbital_parameters_path): """ 获取输入文件的路径 计算坡度图、坡向图 计算局部入射角 """ para_names = ["Dem", "IncidenceAngle", "OrbitalParameters", "经验A"] if len(para_names) == 0: return False # 获取三个文件的路径 # print("输出三个文件路径",Dem_path,IncidenceAngle_path,OrbitalParameters_path) # 确定坡度、坡向的输出路径，输出坡度、坡向图 slope_out_path = r"D:\MicroWorkspace\LeafAreaIndex\Temporary\UnClipslope.tif" aspect_out_path = r"D:\MicroWorkspace\LeafAreaIndex\Temporary\UnClipaspect.tif" print("slope_out_path的路径是", slope_out_path) print("aspect_out_path的路径是", aspect_out_path) self.creat_twofile(dem_path, slope_out_path, aspect_out_path) # 根据入射角文件对坡度坡向图进行裁剪与重采样 slope_out_path2 = r"D:\MicroWorkspace\LocaLangle\Temporary\Clipslope.tif" aspect_out_path2 = r"D:\MicroWorkspace\LocaLangle\Temporary\Clipaspect.tif" self.resampling(slope_out_path, aspect_out_path, incidence_angle_path, slope_out_path2, aspect_out_path2) # 输出局部入射角文件 save_localangle_path = r"D:\\MicroWorkspace\\LocaLangle\\Temporary\\\localangle.tif" self.clau(slope_out_path2, aspect_out_path2, incidence_angle_path, orbital_parameters_path, save_localangle_path) # if __name__ == '__main__': # calu_incident = CalculateIncident() # Dem_path = "D:\\MicroWorkspace\\LocaLangle\\Input\\dem.tif" # IncidenceAngle_path = "D:\\MicroWorkspace\\LocaLangle\\Input\\RSJ.tif" # OrbitalParameters_path = "D:\\MicroWorkspace\\LocaLangle\\Input\\" \ # "GF3_KAS_FSII_020008_E113.2_N23.1_20200528_L1A_HHHV_L10004829485.meta.xml" # calu_incident.localangle(Dem_path, IncidenceAngle_path, OrbitalParameters_path) # print('done')