415 lines
18 KiB
Python
415 lines
18 KiB
Python
# -*- 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')
|