KJ135
/
microproduct
mirror of http://172.16.0.12:5000/WBJY/microproduct.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
microproduct/vegetationHeight-L-SAR/Orthorectification.py

1045 lines
48 KiB
Python
Raw Permalink Normal View History

入所测试定稿上传
2023-08-28 10:17:29 +00:00
from tool.algorithm.image.ImageHandle import ImageHandler
from tool.algorithm.algtools.ROIAlg import ROIAlg
#from Orthorectification.Indirect_pstn import PSTN
import os
from types import resolve_bases
from typing import OrderedDict
# import gdal
from numpy.ma.core import power, reshape
import xmltodict
from xml.dom.minidom import parse
import pandas
import numpy as np
import json
#from matplotlib import pyplot as plt
import datetime
import yaml
import math
from decimal import Decimal # 为了更好的精度
# from multiprocessing import Pool # 进程池
import multiprocessing # 多进程类
# import numba
# from numba import jit
from SatelliteOrbit import ReconstructionSatelliteOrbit
from util import GetVectorNorm,FindInfomationFromJson,appendInfomationForCsv
#from util import XYZ_to_LLA,LLA_to_XYZ
#import indirect_pstn
#import scipy as sp
#import pylab as pl
from scipy.optimize import leastsq
def Interpolation_sinc():
''' 基于sinc的插值算法，专门用于影像校正
'''
pass
class Orthorectification(object):
'''
正射校正模块类。
'''
def __init__(self,configPath="./config.ymal") -> None:
super().__init__()
# 常量声明
# 加载配置文件
with open(configPath,'r',encoding='utf-8') as f:
cont = f.read()
self.config = yaml.load(cont,Loader=yaml.FullLoader)
self.config['SatelliteOribit']['StartTime']=datetime.datetime.strptime(self.config['SatelliteOribit']['StartTime']['Value'],self.config['SatelliteOribit']['StartTime']['format']).timestamp()
self.R_e=self.config['SatelliteOribit']['ReferenceSpheroid']['Re'] # m
self.R_p=self.config['SatelliteOribit']['ReferenceSpheroid']['Rp'] # m
self.w_e=self.config['SatelliteOribit']['ReferenceSpheroid']['we'] # rad/s
def ParseHearderFile(self,HeaderFilePath_str,GPSNodeNames_list=['TimeStamp', 'xPosition', 'yPosition', 'zPosition', 'xVelocity', 'yVelocity', 'zVelocity']):
'''
解析头文件，返回计算所需要的必要信息:
args:
HeaderFilePath_str: 头文件地址，
return:
HeaderFileInfomation_json:头文件信息包
'''
with open(HeaderFilePath_str,'r',encoding='utf-8') as fp:
HeaderFile_dom_str=fp.read()
HeaderFile_dom=xmltodict.parse(HeaderFile_dom_str)
HeaderFile_dom_json=json.loads(json.dumps(HeaderFile_dom))
# 获取头文件
HeaderInformation_json={}
# 解析轨道节点，假定是GPS节点。
HeaderInformation_json['GPS']=[]
# GPS 解析信息
GPSNode_Path=self.config["GPS"]['GPSNode_Path']
GPSNodeNames_list=self.config["GPS"]['NodeInfomation_Name']
GPSTime_Format=self.config['GPS']['Time_format']
GPSPoints=FindInfomationFromJson(HeaderFile_dom_json,GPSNode_Path)
for GPSPoint in GPSPoints:
GPSPoint=[
datetime.datetime.strptime(GPSPoint[GPSNodeNames_list[0]],GPSTime_Format).timestamp(), # TimeStamp
float(GPSPoint[GPSNodeNames_list[1]]), # Xp
float(GPSPoint[GPSNodeNames_list[2]]), # Yp
float(GPSPoint[GPSNodeNames_list[3]]), # Zp
float(GPSPoint[GPSNodeNames_list[4]]), # Vx
float(GPSPoint[GPSNodeNames_list[5]]), # Vy
float(GPSPoint[GPSNodeNames_list[6]])] # VZ
HeaderInformation_json["GPS"].append(GPSPoint)
# 提取成像时间信息
HeaderInformation_json['ImageInformation']={}
# 1、开始成像时间
HeaderInformation_json['ImageInformation']['StartTime']=datetime.datetime.strptime(
FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['StartImageTime']['NodePath']),
self.config['imageinfo']['StartImageTime']['Format']
).timestamp()
# 2、影像结束成像时间
HeaderInformation_json['ImageInformation']['EndTime']=datetime.datetime.strptime(
FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['EndImageTime']['NodePath']),
self.config['imageinfo']['StartImageTime']['Format']
).timestamp()
# 3、影像的宽高
HeaderInformation_json['ImageInformation']['height']=int(FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['ImageHeight']['NodePath']))
HeaderInformation_json['ImageInformation']['width']=int(FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['ImageWidth']['NodePath']))
# 4、影像的近斜距
HeaderInformation_json['ImageInformation']['NearRange']=float(FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['NearRange']['NodePath']))
# 5、入射角
HeaderInformation_json['ImageInformation']['incidenceAngle']={
'NearRange':float(FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['incidenceAngle']['NearRange']['NodePath'])),
'FarRange':float(FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['incidenceAngle']['FarRange']['NodePath']))
}
# 6、成像时刻影像带宽
HeaderInformation_json['ImageInformation']['bandWidth']=float(FindInfomationFromJson(HeaderFile_dom_json,self.config['sensor']['bandWidth']['NodePath']))
# 7、多普勒质心常数
DopplerCentroidCoefficients_list=FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['DopplerCentroidCoefficients']['NodePath'])
HeaderInformation_json['ImageInformation']['DopplerCentroidCoefficients']=[
float(DopplerCentroidCoefficients_list[self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][0]]),
float(DopplerCentroidCoefficients_list[self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][1]]),
float(DopplerCentroidCoefficients_list[self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][2]]),
float(DopplerCentroidCoefficients_list[self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][3]]),
float(DopplerCentroidCoefficients_list[self.config['imageinfo']['DopplerCentroidCoefficients']['DopplerCentroidCoefficients_Name'][4]]),
]
# 8、波长
HeaderInformation_json['ImageInformation']['lambda']=float(FindInfomationFromJson(HeaderFile_dom_json,self.config['sensor']['lambda']['NodePath']))
# 9、中心像元对应得中心经纬度
CenterPoint=FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['CenterImagePositon']['NodePath'])
HeaderInformation_json['ImageInformation']['centerPosition']=[
float(CenterPoint[self.config['imageinfo']['CenterImagePositon']['Value'][0]]), # 纬度
float(CenterPoint[self.config['imageinfo']['CenterImagePositon']['Value'][1]]), # 经度
]
# 10、多普勒参数参考时间
HeaderInformation_json['ImageInformation']['DopplerParametersReferenceTime']=float(FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['DopplerParametersReferenceTime']['NodePath']))
# 11、参考斜距
HeaderInformation_json['ImageInformation']['refRange']=float(FindInfomationFromJson(HeaderFile_dom_json,self.config['imageinfo']['ReferenceRange']['NodePath']))
HeaderInformation_json['PRF']=float(FindInfomationFromJson(HeaderFile_dom_json,self.config['sensor']['PRF']['NodePath']))
return HeaderInformation_json
pass
def getRByColumncode(self,c):
'''
根据列号计算斜距
args:
c：列号
'''
return self.R0+c*self.delta_R
def getTimeByLineCode(self,r):
'''
根据行号计算成像时间
args:
r：行号
'''
return np.matmul(np.concatenate([ np.array([r**i]).reshape(1,1) for i in range(self.r2t_A_arr.shape[0])],axis=1),self.r2t_A_arr)
def ReconstuctionTimesOflyDirectionPositionModel(self,timePoints_list):
''' 构建飞行向坐标与时间的转换模型,注意这里没有调整时间的起算点。
args:
timePoints_list:时间点坐标 [r,t]
'''
# 根据点数确定模型，最高为3次项模型
Y=np.zeros((len(timePoints_list),1)) # 确定Y
X=np.zeros((len(timePoints_list),1))
for i in range(len(timePoints_list)):
Y[i,0]=timePoints_list[i][1]
X[i,0]=timePoints_list[i][0]
if len(timePoints_list)==2:
# 一次项
X=np.concatenate([np.ones((len(timePoints_list),1)),X],axis=1)
A=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)
elif len(timePoints_list)>=3: # 二次项
X=np.concatenate([np.ones(len(timePoints_list),1),X,np.power(X,2)],axis=1)
A=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)
pass
elif len(timePoints_list)>=8: # 三次项
X=np.concatenate([np.ones(len(timePoints_list),1),X,np.power(X,2),np.power(X,3)],axis=1)
A=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)
pass
self.r2t_A_arr=A
def PrepareConverteSystem(self):
''' 计算常量
在数据计算之前，提前处理部分常量信息
1、构建r,c 到坐标斜距的计算公式
2、计算行号对应的时间频率
'''
R0=self.header_info['ImageInformation']['NearRange'] # 起始斜距，近斜距
delta_R=self.config['LightSpeed']/(2*self.header_info['ImageInformation']['bandWidth']*10**6)# 计算距离向分辨率 ,参考《InSAR原理与应用》P36，距离向分辨率 带宽单位为MHZ
self.ReconstuctionTimesOflyDirectionPositionModel([
[0,self.header_info['ImageInformation']['StartTime']],
[self.header_info['ImageInformation']['height']-1,self.header_info['ImageInformation']['EndTime']]
]) # 构建坐标到时间的变换
self.R0=R0
self.delta_R=delta_R
# 计算行号分辨率
self.delta_t=1/self.header_info['PRF']
def Orthorectification(self,FilePath_str,ori_sim):
'''
正射校正组件。
正射校正整体可以分为两个步骤：
Step 1：计算出栅格坐标对应的真实的大地坐标系坐标
Step 2：根据大地坐标系将影像重采样到大地坐标系空间中。
args:
FilePath_str:影像所在文件夹
'''
# 分离需要校正对象
file_name_list=os.listdir(FilePath_str)
header_name=[file_name for file_name in file_name_list if file_name.rfind(".meta.xml")==len(file_name)-len('.meta.xml')][0] # 头文件
tiff_name=[file_name for file_name in file_name_list if file_name.rfind(".tiff")==len(file_name)-len('.tiff')] # 影像文件
# 解析头文件
self.header_info=self.ParseHearderFile(os.path.join(FilePath_str,header_name))
# 构建轨道模型
self.SatelliteOrbitModel=ReconstructionSatelliteOrbit(self.header_info['GPS'],starttime=self.header_info['ImageInformation']['StartTime']) # 构建卫星轨道模型,取第0个节点的时间
# 求解模型前计算常数变换模型
self.PrepareConverteSystem()
# # 测试第0，0个点的空间坐标
# self.ConverteCoordinaryPoint_single(0,0)
# 批量计算空间坐标
lamda = self.header_info['ImageInformation']['lambda']
height_array,width_array,width,height,geo,pro=self.analy_image(ori_sim)
p0_array=np.zeros((3,height, width),dtype=float) # 创建3波段空数组存放p0坐标
sat_array=np.zeros((3,height, width),dtype=float) # 创建3波段空数组存放卫星位置
im_data, im_data2=self.ConverteCoordinary(height_array, width_array,width,height, geo, pro, p0_array, sat_array)
return im_data, im_data2, lamda
def analy_image(self, ori_sim):
"""
计算影像的行列号的范围、投影系、坐标系
"""
height_array = ImageHandler.get_band_array(ori_sim, 1)
width_array = ImageHandler.get_band_array(ori_sim, 2)
geo = ImageHandler.get_geotransform(ori_sim)
pro = ImageHandler.get_projection(ori_sim)
width=ImageHandler.get_img_width(ori_sim)
height=ImageHandler.get_img_height(ori_sim)
return height_array,width_array,width,height,geo,pro
def ConverteCoordinary(self, height_array, width_array,width,height, geo, pro, p0_array, sat_array):
'''
批量求解点坐标
'''
return p0_array, sat_array
pass
def ConverteCoordinaryPoint_single(self,r,c,h=0,breakCondition=1):
'''
求解单个点坐标
args:
r: 行号 height
c: 列号 width
h：平均大地高
breakCondition: 迭代终止条件，默认值为1
'''
pass
class Orthorectification_RD(Orthorectification):
'''
正射校正模块类。
'''
def __init__(Orthorectification_RD,configPath="./config.ymal") -> None:
super().__init__(configPath)
def ConverteCoordinary(self, height_array, width_array,width,height, geo, pro, p0_array, sat_array):
'''
param： height_min:
批量求解点坐标
'''
# 构建映射表 ----------------------------------需要重写代码
result={}
print("开始时间：",datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S.%f"))
alpha0=0
beta0=0
# for r in range(self.header_info['ImageInformation']['height']):
# if r%10==0:
# print("r={} 时间：".format(str(r)),datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S.%f"))
# for c in range(self.header_info['ImageInformation']['width']):
# result[(r,c)]=self.ConverteCoordinaryPoint_single(r,c,alpha0=alpha0,beta0=beta0)
# alpha0=result[(r,c)][1]
# beta0=result[(r,c)][2]
# print("终止时间：", datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S.%f"))
# 计算ori_smi对应的卫星空间位置及地面p0坐标，并存入三维数组中。
for r in range(0, height):
if r%10==0:
print("r={} 时间：".format(str(r)),datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S.%f"))
for c in range(0, width):
num_r=height_array[r,c]
num_c=width_array[r,c]
result=self.ConverteCoordinaryPoint_single(num_r,num_c,alpha0=alpha0,beta0=beta0)
xyz=result[0]
alpha0=result[1]
beta0=result[2]
Sate_x, Sate_y, Sate_z=result[3],result[4],result[5]
p0_array[0,r, c] = xyz[0, 0] # 存放p0坐标(地固参心坐标)
p0_array[1,r, c] = xyz[0, 1]
p0_array[2,r, c] = xyz[0, 2]
sat_array[0,r, c] = Sate_x # 存放卫星位置(地固参心坐标)
sat_array[1,r, c] = Sate_y
sat_array[2,r, c] = Sate_z
print("终止时间：", datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S.%f"))
return p0_array, sat_array
pass
# def ConverteCoordinary(self,OutputPath=None,DEM=None,BlocksPath="./OrthorecticationGrid",BlockSize=1):
# '''批量转换坐标
# 采用矩阵结构，批量转换点，并将坐标点保存到指定文件分块文件中，并最终汇总成一个大型文件。
# 分块文件夹内容格式：
#
# Args:
# OutputPath: 最终文件地址，如果为None，则不把分块文件进行汇总成一个整体文件。
# DEM：与坐标对应的高程，default:None
# BlocksPath:分块存储地址与最终集合地址
# BlockSize:块大小，default：1 GB
# return:
# Result_path: 最终文件地址
# '''
# # 数据参数预读取
# # 多普勒参数
# Dopper_d=np.array(self.header_info['ImageInformation']['DopplerCentroidCoefficients']).reshape(-1,1)
# LightSpeed=self.config['LightSpeed']
# # T0
# DopplerParametersReferenceTime=self.header_info['ImageInformation']['refRange']*2/LightSpeed
# # 波长
# lamda=self.header_info['ImageInformation']['lambda']
# # Re,Rp,we
# Re=self.config['SatelliteOribit']['ReferenceSpheroid']['Re']
# Rp=self.config['SatelliteOribit']['ReferenceSpheroid']['Rp']
# we=self.config['SatelliteOribit']['ReferenceSpheroid']['we'] # we 地球自转速度矢量 --- 废弃
# # 影像宽高
# height=self.header_info['ImageInformation']['height']
# width=self.header_info['ImageInformation']['width']
# # 粗估分块文件的大小
# size_max=BlockSize*1024*1024*1024 # 粗估最大字节数
# col_num_block=int(size_max/(height*3*8))-1 # 预留一行作为其他信息的存储空间
# alpha0=0
# beta0=0
# for r in range(self.header_info['ImageInformation']['height']):
# if r%10==0:
# print("r={} 时间：".format(str(r)),datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S.%f"))
# for c in range(self.header_info['ImageInformation']['width']):
# result=self.ConverteCoordinaryPoint_single(r,c,alpha0=alpha0,beta0=beta0)
# alpha0=result[1]
# beta0=result[2]
# print("终止时间：",datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S.%f"))
# # 构建分块矩阵计算代码
# pass
def ConverteCoordinaryPoint_single(self,r,c,h=0,breakCondition=1,alpha0=0,beta0=0):
'''
求解单个点坐标
args:
r: 行号 height
c: 列号 width
h：平均大地高
breakCondition: 迭代终止条件，默认值为1
return:
TargetPosition_:点坐标
版本改动：
2021.8.19 尝试将椭球方程进行变形
2021.8.22 方法重写
'''
# 计算对应成像时间和成像斜距
R=self.getRByColumncode(c) # 斜距
t=self.getTimeByLineCode(r) # 成像时间
# 解析成像时刻，卫星空间位置
SatelliteState=self.SatelliteOrbitModel.SatelliteSpaceState(t)
[Sate_time,Sate_x,Sate_y,Sate_z,Sate_vx,Sate_vy,Sate_vz]=SatelliteState.reshape(1,7).tolist()[0]
# 多普勒
Dopper_d=np.array(self.header_info['ImageInformation']['DopplerCentroidCoefficients']).reshape(-1,1)
LightSpeed=self.config['LightSpeed']
# DopplerParametersReferenceTime=self.header_info['ImageInformation']['DopplerParametersReferenceTime']
DopplerParametersReferenceTime=self.header_info['ImageInformation']['refRange']*2/LightSpeed
# 波长
lamda=self.header_info['ImageInformation']['lambda']
# Re,Rp,we
Re=self.config['SatelliteOribit']['ReferenceSpheroid']['Re']
Rp=self.config['SatelliteOribit']['ReferenceSpheroid']['Rp']
we=self.config['SatelliteOribit']['ReferenceSpheroid']['we'] # we 地球自转速度矢量
# 使用牛顿迭代法，求解Xp,Yp,Zp
# 椭球体方程常数 a+h ==>a b==>b
a=(Re+h)**2
b=Rp**2
# 初始迭代点为中心像元经纬度对应的大地坐标系
#XYZ_=self.LLA_to_XYZ(self.header_info['ImageInformation']['centerPosition'][0],
#self.header_info['ImageInformation']['centerPosition'][1],
#h)
XYZ=np.array([Sate_x,Sate_y,Sate_z]).reshape(1,3)
# 确定初始值，作为求解起始值
TargetPosition=math.sqrt( 1/((Sate_x**2+Sate_y**2)/a+(Sate_z**2)/b))
TargetPosition=TargetPosition*np.array([Sate_x,Sate_y,Sate_z]).reshape(1,3)
TargetPosition=TargetPosition.reshape(1,3)
R_threshold= np.sqrt( np.sum( np.power(TargetPosition-XYZ,2) ))
SatellitePosition=np.array([Sate_x,Sate_y,Sate_z]).reshape(3,1)
SatelliteVelocity=np.array([Sate_vx,Sate_vy,Sate_vz]).reshape(3,1)
TargetPosition_,alpha,beta=self.ASF(R,SatellitePosition,SatelliteVelocity,TargetPosition,Dopper_d,LightSpeed,DopplerParametersReferenceTime,lamda,R_threshold,H=h,alpha0=alpha0,beta0=beta0)
# 中间检查有无错误
# TargetPositionLLA=XYZ_to_LLA(TargetPosition[0,0],TargetPosition[0,1],TargetPosition[0,2]).reshape(1,-1)
# TargetPosition_LLA=XYZ_to_LLA(TargetPosition_[0,0],TargetPosition_[0,1],TargetPosition_[0,2]).reshape(1,-1)
np.set_printoptions(suppress=True)
# print(TargetPosition.reshape(1,-1))
# print(TargetPosition_.reshape(1,-1))
# print(TargetPositionLLA.reshape(1,-1))
# print(TargetPosition_LLA.reshape(1,-1))
# print(TargetPosition_-TargetPosition)
# print(TargetPosition_LLA-TargetPositionLLA)
np.set_printoptions(suppress=False)
return TargetPosition_.reshape(1,-1),alpha,beta, Sate_x,Sate_y,Sate_z
'''
RD模型解算模块,这里主要使用ASF方法
reference: 陈尔学. 星载合成孔径雷达影像正射校正方法研究[D].中国林业科学研究院,2004. 第三章内容
且假定地表运动速度为0
'''
def ASF(self,R,SatellitePosition,SatelliteVelocity,TargetPosition,Dopper_d,LightSpeed,DopplerParametersReferenceTime,lamda,R_threshold,H=0,alpha0=0,beta0=0):
# 部分参数初始化
alpha=alpha0
beta=beta0
delta_alpha=0
Dopper_d=Dopper_d.reshape(-1,1)
# 这里需要重构
TSR=2*R/LightSpeed
T0=DopplerParametersReferenceTime
TSR=TSR-T0
# 这个公式是根据《RADARSAT-2 PRODUCT FORMAT DEFINITION》(RadarSat2的数据格式手册)
# “table 5-17 Doppler Rate Values” 中 dopplerRateValuesCoefficients 一行
# ist of up to 5 Doppler rate values coefficients as a function of slant range time: r0, r1, r2, r3, and r4 where:
# Doppler rate values = r0 + r1(tSR - t0) + r2(tSR-t0)^2 + r3(tSR-t0)^3 + r4(tSR-t0)^4
fd=np.matmul(np.array([TSR**0,TSR**1,TSR**2,TSR**3,TSR**4]).reshape(1,5),Dopper_d)
we=self.w_e
M=self.SatelliteCoordination2ECR(SatellitePosition,SatelliteVelocity)
i=0 # 统计迭代次数
while True:
TargetVelocity=np.array([0,0,0]).reshape(3,1)# np.array([-1*we*TargetPosition[1,0],we*TargetPosition[0,0],0]) # 计算地面速度,粗糙计算时，默认为0
beta=self.GetLookFromRangeYaw(R,alpha,beta,SatellitePosition,SatelliteVelocity,R_threshold,H=H)
FD_=self.FD(alpha,beta,SatelliteVelocity,TargetVelocity,R,lamda,M)
delta_fd=FD_-fd
delta_alpha=-1*delta_fd*lamda/(2* np.sum(np.sqrt(np.power(SatelliteVelocity-TargetVelocity,2))))
if np.abs(delta_alpha*R)<0.1:
XYZ=self.GetXYZByBetaAlpha(alpha,beta,SatellitePosition,R,M)
break
alpha=alpha+delta_alpha
i=i+1
return XYZ.reshape(1,-1),alpha,beta
def GetLookFromRangeYaw(self,R,alpha,beta,SatellitePosition,SatelliteVelocity,R_threshold,H=0):
'''
根据R和alpha 计算出beta来。
根据参考论文，式3-13，完成此方法。（注意这是近似法
args:
R：斜距
alpha:雷达侧视角
beta0:beta初始值
SatellitePosition:3x1 卫星空间位置
TargetPosition：3x1 目标对象
'''
M=self.SatelliteCoordination2ECR(SatellitePosition,SatelliteVelocity)
# 地球半径
Rp=self.R_p+H
Rs_norm=math.sqrt(np.matmul(SatellitePosition.T,SatellitePosition).reshape(1,1)[0,0])
# 初始化beta
if beta==0:
beta_cos=(Rs_norm**2+R**2-Rp**2)/(2*Rs_norm*R)
beta=math.acos(beta_cos)
# 迭代
i=0
while True:
# 计算斜率的增加
delta_R=R-self.FR(alpha,beta,SatellitePosition,M,R_threshold,R,H=H)
# 计算入射角
sin_eta=Rs_norm*math.sin(beta)/Rp
tan_eta=sin_eta/math.sqrt(1-sin_eta**2)
# 计算增量
delta_beta=delta_R/(R*tan_eta)
# 更新
beta=beta+delta_beta
# 判断循环终止条件
if np.abs(delta_R)<0.1:
break
i=i+1
if i>=10000: # 达到终止条件
return None
return beta
def XYZOuter(self,A,B):
'''
外积（叉乘）,日后版本换成可以任意维度的外积运算方程
args:
A:3x1
B:3x1
'''
C=np.cross(A.reshape(1,3),B.reshape(1,3)).reshape(-1,1)
return C
def SatelliteCoordination2ECR(self,SatellitePosition,SatelliteVelocity):
'''
将卫星为中心的直角坐标系Sxyz转换为ECR坐标系。
M=[x,y,z]
args:
SatellitePosition:np.array 3x1 ECR坐标系下，卫星空间位置
SatelliteVelocity:np.array 3x1 ECR坐标系下，卫星运动速度
return:
坐标系变换矩阵 M
ECR=M*Sxyz
P‘=MP
'''
z=SatellitePosition/GetVectorNorm(SatellitePosition)
y=self.XYZOuter(SatellitePosition,SatelliteVelocity)/(np.sqrt(np.matmul(SatellitePosition.T,SatellitePosition))*np.sqrt(np.matmul(SatelliteVelocity.T,SatelliteVelocity))) # 外积
x=self.XYZOuter(y,z)
M=np.concatenate([x.reshape(3,1),y.reshape(3,1),z.reshape(3,1)],axis=1) # 以列组合
return M
def UnitVectorOfSatelliteAndTarget(self,alpha,beta,M):
'''
获取从卫星指向地面目标的单位矢量
'''
# 目标T的单位向量 ST=R*M*P
P=np.array([math.sin(alpha),-math.cos(alpha)*math.sin(beta),-math.cos(alpha)*math.cos(beta),]).reshape(3,1)
P_=np.matmul(M,P).reshape(-1,1)
return P_
def FD(self,alpha,beta,SatelliteVelocity,TargetVelocity,R,lamda,M):
'''
根据beta,alpha,卫星空间位置,斜距，转换矩阵
'''
we=self.w_e
#
TargetVelocity=0 # 假设地面运动为0
UnitVector=self.UnitVectorOfSatelliteAndTarget(alpha,beta,M)
Rst=-1*R*UnitVector
Rst_Vst=np.matmul(Rst.T,SatelliteVelocity-TargetVelocity)
fd=-1*(2/lamda)*Rst_Vst/R
return fd
def GetXYZByBetaAlpha(self,alpha,beta,SatellitePosition,R,M):
'''
从beta，alpha获取获取目标的空间位置
'''
P=np.array([math.sin(alpha),-math.cos(alpha)*math.sin(beta),-math.cos(alpha)*math.cos(beta),]).reshape(3,1)
Rt=SatellitePosition+R*np.matmul(M,P)
return Rt
def FR(self,alpha,beta,SatellitePosition,M,R_threshold,R_ref,H=0):
'''
根据 雷达侧视角，雷达视角，卫星位置，坐标系变换矩阵，大地高等参数，根据参考论文中公式3-1 ==》3-10的计算过程，设计此方程
args:
alpha:雷达侧视角 卫星与地物间速度运动差导致的（即多普勒频移的结果）
beta: 雷达视角
Satellite_position:卫星空间位置
M：坐标系转换矩阵
H：大地高
return:
R：斜距
'''
UnitVector=self.UnitVectorOfSatelliteAndTarget(alpha,beta,M)
# 解方程 3-10 纠正错误，A公式有问题
a=(self.R_e+H)**2
b=self.R_p**2
[Xs,Ys,Zs]=SatellitePosition.reshape(-1,1)
[Xp,Yp,Zp]=UnitVector.reshape(-1,1)
A=(Xp**2+Yp**2)/a+(Zp**2)/b
B=2*(Xs*Xp+Ys*Yp)/a+2*Zs*Zp/b
C=(Xs**2+Ys**2)/a+Zs**2/b
C=C-1
t=B**2-4*A*C
#
R1=(-B-math.sqrt(t))/(2*A)
if R1>R_threshold:
return R1
else:
return None
#return (-B+math.sqrt(t))/(2*A) # 这里通过几何分析排除这个解
return None
class Orthorectification_PSTN(Orthorectification):
'''间接定位法，生产GTC
'''
def Orthorectification(self,FilePath_str, ori_sim):
'''
正射校正组件。
正射校正整体可以分为两个步骤：
Step 1：计算出栅格坐标对应的真实的大地坐标系坐标
Step 2：根据大地坐标系将影像重采样到大地坐标系空间中。
args:
FilePath_str:影像所在文件夹
'''
# 分离需要校正对象
file_name_list=os.listdir(FilePath_str)
header_name=[file_name for file_name in file_name_list if file_name.rfind(".meta.xml")==len(file_name)-len('.meta.xml')][0] # 头文件
tiff_name=[file_name for file_name in file_name_list if file_name.rfind(".tiff")==len(file_name)-len('.tiff')] # 影像文件
# 解析头文件
self.header_info=self.ParseHearderFile(os.path.join(FilePath_str,header_name))
# 构建轨道模型
self.SatelliteOrbitModel=ReconstructionSatelliteOrbit(self.header_info['GPS'],starttime=self.config['SatelliteOribit']['StartTime']) # 构建卫星轨道模型,取第0个节点的时间
# 求解模型前计算常数变换模型
self.PrepareConverteSystem()
# 测试第0，0个点的空间坐标
Doppler_d=np.array(self.header_info['ImageInformation']['DopplerCentroidCoefficients']).reshape(-1,1)
LightSpeed=self.config['LightSpeed']
T0=self.header_info['ImageInformation']['refRange']*2/LightSpeed
# -2312340.457 5393726.182 2489909.349
TargetPosition=np.array([-2312340.457, 5393726.182 ,2489909.349]).reshape(1,-1)
#TargetPosition=np.array([-2366713.39445398 , 5350367.51812531 , 2531732.36706542]).reshape(1,-1)
lamda=self.header_info['ImageInformation']['lambda']
StartTime=self.header_info['ImageInformation']['StartTime']
PRF=self.header_info['PRF']
#self.PSTN(TargetPosition,Doppler_d,StartTime,lamda,T0,LightSpeed,PRF)
def ConverteCoordinary(self,DEMClip):
'''
批量求解点坐标
args:
DEMClip:裁剪后的DEM
return:
rcArray:像元的对应的行列号
'''
pass
def ConverteCoordinaryPoint_single(self,r,c,h=0,breakCondition=1):
'''
求解单个点坐标
args:
r: 行号 height
c: 列号 width
h：平均大地高
breakCondition: 迭代终止条件，默认值为1
'''
pass
def PSTN(self,TargetPostion,Doppler_d,StartTime,lamda,T0,LightSpeed,PRF):
'''间接定位法
Args:
TargetPostion:地面坐标
Doppler_d:多普勒系数
StartTime:成像开始时间
lamda: 波长
T0: 多普勒参考时间
LightSpeed: 光速
return:
r,c:影像的行列号
'''
ti=StartTime
Doppler_d=Doppler_d.reshape(-1,1)
dt=0.001
delta_t=1/PRF
delta_R=self.delta_R
i=0
while True:
# 卫星轨道
statellitePosition=self.SatelliteOrbitModel.SatelliteSpaceState(ti) # 卫星轨道
statellite_XYZ=statellitePosition[0,1:4]
statellite_velocity=statellitePosition[0,4:]
# 位置向量
R_sp=statellite_XYZ-TargetPostion.reshape(1,-1)
V_sp=statellite_velocity
inc_r=GetVectorNorm(V_sp)
# 斜距
R=GetVectorNorm(R_sp)
# 计算多普勒频率
fde=-(2/lamda)*np.matmul(R_sp.reshape(1,3),V_sp.reshape(3,1))/R # 理论
TSR=R*2/LightSpeed-T0
fd=np.matmul(np.array([TSR**0,TSR**1,TSR**2,TSR**3,TSR**4]).reshape(1,5),Doppler_d) #t=ti
# 计算多普勒频率的导数
t=ti+dt
satellitionPostion_t=self.SatelliteOrbitModel.SatelliteSpaceState(t)
sate_position=satellitionPostion_t[0,1:4]
sate_velocity=satellitionPostion_t[0,4:]
R_sp_t=sate_position-TargetPostion
V_sp_t=sate_velocity
R_t=GetVectorNorm(R_sp_t)
fd_dt=-(2/(lamda*R_t))*np.matmul(R_sp_t.reshape(1,3),V_sp_t.reshape(3,1))
fd_grad=(fd_dt-fde)/dt
inc_t=(fde-fd)/fd_grad # 求得变量
if np.abs(inc_t)<delta_t*0.001:
break
ti=ti-inc_t
i=i+1
# 计算行号
r=(ti-StartTime)/delta_t
c=(R-self.header_info['ImageInformation']['NearRange'])/delta_R
# 取整
r=np.round(r)
c=np.round(c)
return r,c
class fine_registration():
def __init__(self):
pass
def choose_same_point(self, master_sim_file, auxil_sim_file):
"""
master_sim_file、auxil_sim_file是经过配准裁剪后的，具有相同的范围，映射关系为（rm,cm）->(rs,cs)
:param master_sim_file：主影像的ori_sim.tif路径（正射输出的成果）
:param auxil_sim_file：辅影像的ori_sim.tif路径（正射输出的成果）
return 数组rm、cm、rs、cs
"""
width= ImageHandler.get_img_width(master_sim_file) #8182
height=ImageHandler.get_img_height(master_sim_file)
img_geotrans = ImageHandler.get_geotransform(master_sim_file)
pro = ImageHandler.get_projection(master_sim_file)
master_height_array = ImageHandler.get_band_array(master_sim_file, 1)
master_width_array = ImageHandler.get_band_array(master_sim_file, 2)
auxil_height_array = ImageHandler.get_band_array(auxil_sim_file, 1)
auxil_width_array = ImageHandler.get_band_array(auxil_sim_file, 2)
# 判断影像的大小，1.若影像的行列小于6 无法选择同名点；2.影像行列在6-2000，则每隔两个点取一个同名点；3.影像行列在2000以上，每隔20个点取一个同名点
if (width < 6) and (height < 6):
raise Exception('overlap so small,less of 9 pixel!')
elif (width > 6 and height > 6) and (width < 2000 and height < 2000) : # 取一半的同名点
point_dict1={}
wid_param1=int(width/2)
hig_param1 = int(height/2)
for i in range(1,hig_param1+1):
for j in range(1,wid_param1+1):
hig_value = (height - 2) * (i / hig_param1)
wid_value = (width - 2) * (j / wid_param1)
value=np.array([hig_value,wid_value])
# print(value)
name="a"+str(i)+str(j)
point_dict1.update({name:value})
elif (width > 2000) and (height > 2000): # 取一半的同名点
point_dict1 = {}
wid_param1 = int(width / 20)
hig_param1 = int(height / 20)
for i in range(1, hig_param1 + 1):
for j in range(1, wid_param1 + 1):
hig_value = (height - 2) * (i / hig_param1)
wid_value = (width - 2) * (j / wid_param1)
value = np.array([hig_value, wid_value])
name = "a" + str(i) +"_"+ str(j)
point_dict1.update({name: value})
# 将同名点的行、列分别存到不同的数组中
master_rm=[]
auxil_rs=[]
master_cm=[]
auxil_cs=[]
for key,value in point_dict1.items():
x,y=int(value[0]),int(value[1])
# print((x,y))
master_rm_array = master_height_array[x, y] # 主影像行列6006,8182
master_cm_array = master_width_array[x, y]
auxil_rs_array = auxil_height_array[x, y] # 辅影像行列
auxil_cs_array = auxil_width_array[x, y]
if master_rm_array>1 and master_cm_array>1 and auxil_rs_array>1 and auxil_cs_array>1: # 剔除边界值
master_rm.append(master_rm_array)
master_cm.append(master_cm_array)
auxil_rs.append(auxil_rs_array)
auxil_cs.append(auxil_cs_array)
return np.array(master_rm), np.array(master_cm), np.array(auxil_rs), np.array(auxil_cs)
def fun(self,p1,rm,cm):
"""
最小二乘法的多项式公式
:param p1:系数的数组
:param rm：同名点的列
:param cm：同名点的行
"""
a0, a1, a2, a3, a4, a5=p1
return a0+a1*rm+a2*cm+a3*rm*rm+a4*cm*cm+a5*rm*cm
def error(self,p1,rm,cm,y):
"""
误差公式
:param p1:系数的数组
:param rm：同名点的列
:param cm：同名点的行
:param y：误差值
"""
error_value=self.fun(p1,rm,cm)-y
return error_value
def fun2(self,p2,rm,cm):
b0, b1, b2, b3, b4, b5=p2
return b0+b1*rm+b2*cm+b3*rm*rm+b4*cm*cm+b5*rm*cm
def error2(self,p2,rm,cm,y):
error_value2=self.fun2(p2,rm,cm)-y
return error_value2
def least_sqrt(self, master_sim_file, auxil_sim_file):
"""
最小二乘法拟合多项式
"""
x1,x2,y1,y2=self.choose_same_point(master_sim_file, auxil_sim_file) # 获取同名点的矩阵
p0=np.array([0,0,0,0,0,0]) # 指定多项式系数的初始值
r_plsq=leastsq(self.error,p0,args=(x1,x2,y1)) # 最小二乘法拟合含有多变量的多项式
# print(r_plsq[0])
c_plsq2=leastsq(self.error2,p0,args=(x1,x2,y2))
# print(c_plsq2[0])
return r_plsq[0], c_plsq2[0]
def get_new_master_sim(self, r_plsq, c_plsq2, master_ori_sim, out_new_master_ori_sim):
"""
主影像的ori_sim数据，通过变换模型生成新的ori_sim数据
:param r_plsq:列变换模型系数矩阵
:param c_plsq2：行变换模型系数矩阵
:param master_ori_sim：主影像的ori_sim数据（正射输出的结果）
:param out_new_master_ori_sim：输出由主影像生成的ori_sim数据
"""
height=ImageHandler.get_img_height(master_ori_sim)
width= ImageHandler.get_img_width(master_ori_sim)
img_geotrans = ImageHandler.get_geotransform(master_ori_sim)
img_pro = ImageHandler.get_projection(master_ori_sim)
new_array=np.zeros((2,height,width),dtype=float)
master_height_rm = ImageHandler.get_band_array(master_ori_sim, 1)
master_width_cm = ImageHandler.get_band_array(master_ori_sim, 2)
rm=master_height_rm
cm=master_width_cm
a0, a1, a2, a3, a4, a5=r_plsq[0],r_plsq[1],r_plsq[2],r_plsq[3],r_plsq[4],r_plsq[5]
b0, b1, b2, b3, b4, b5=c_plsq2[0],c_plsq2[1],c_plsq2[2],c_plsq2[3],c_plsq2[4],c_plsq2[5],
rs_array=a0+a1*rm+a2*cm+a3*rm*rm+a4*cm*cm+a5*rm*cm
cs_array=b0+b1*rm+b2*cm+b3*rm*rm+b4*cm*cm+b5*rm*cm
new_array[0,:,:]=rs_array
new_array[1,:,:]=cs_array
ImageHandler.write_img(out_new_master_ori_sim, img_pro, img_geotrans, new_array, no_data='null')
def calu_vector_angle(self, S1, S2, P):
"""
计算垂直基线
S1:s1坐标
S2:s2坐标
R1：P坐标
输出基线B
"""
# 计算垂直基线绝对值
S1_P_x=P[0,:,:]-S1[0,:,:] # 矢量S1-P在x方向
S1_P_y=P[1,:,:]-S1[1,:,:] # 矢量S1-P在y方向
S1_P_z=P[2,:,:]-S1[2,:,:] # 矢量S1-P在z方向
S1_S2_x=S2[0,:,:]-S1[0,:,:] # 矢量S1-S2在x方向
S1_S2_y=S2[1,:,:]-S1[1,:,:] # 矢量S1-S2在y方向
S1_S2_z=S2[2,:,:]-S1[2,:,:] # 矢量S1-S2在z方向
S2_P_x=P[0,:,:]-S2[0,:,:]
S2_P_y=P[1,:,:]-S2[1,:,:]
S2_P_z=P[2,:,:]-S2[2,:,:]
P_0_x=-P[0,:,:] # 矢量地心点-P在x方向
P_0_y=-P[1,:,:] # 矢量地心点-P在y方向
P_0_z=-P[2,:,:] # 矢量地心点-P在z方向
S1P_S1S2=S1_P_x*S1_S2_x+S1_P_y*S1_S2_y+S1_P_z*S1_S2_z # 矢量的乘积S1P*S2P
len_S1_P=np.sqrt(S1_P_x*S1_P_x+S1_P_y*S1_P_y+S1_P_z*S1_P_z) # 矢量的模S1_P0
len_S1_S2=np.sqrt(S1_S2_x*S1_S2_x+S1_S2_y*S1_S2_y+S1_S2_z*S1_S2_z) # 矢量的模S1_S2
cos_angle_s2s1p=S1P_S1S2/(len_S1_P*len_S1_S2) # 夹角S1 余弦值
sin_angle_s2s1p = np.sqrt(1 - cos_angle_s2s1p * cos_angle_s2s1p) # 夹角S1 正弦值
B_vertical=len_S1_S2*sin_angle_s2s1p # 垂直基线的绝对值
# 判断垂直基线的符号（地球观测与导航技术丛书InSAR原理与应用.pdf）68页
len_P_0=np.sqrt(P_0_x*P_0_x+ P_0_y* P_0_y+P_0_z*P_0_z) # 矢量的模P0_P0
S1P_P0=S1_P_x*P_0_x+S1_P_y*P_0_y+S1_P_z*P_0_z
cos_angle1 = S1P_P0 / (len_S1_P * len_P_0)
S2P_P0=S2_P_x*P_0_x+S2_P_y*P_0_y+S2_P_z*P_0_z
len_S2_P=np.sqrt(S2_P_x*S2_P_x+S2_P_y*S2_P_y+S2_P_z*S2_P_z) # 矢量的模S1_P0
cos_angle2 =S2P_P0/(len_S2_P*len_P_0)
# 垂直基线绝对值x垂直基线的符号
cos_angle=cos_angle1-cos_angle2 # 若cos_angle1<cos_angle2 垂直基线的符号f = 1 ;反之等于-1
pos_neg=np.where(cos_angle<0,1,-1)
B_vertical = B_vertical * pos_neg
return B_vertical, len_S1_P, len_S2_P
def calu_phi(self, r1_file, r2_file, lamda_01):
"""
计算phi
输出数组
"""
r1 = ImageHandler().get_data(r1_file)
r2 = ImageHandler().get_data(r2_file)
phi_land = 4 * math.pi * (r1 - r2) * lamda_01
return phi_land
def calu_kz(self, master_incident_file, B_vertical_path,r1_file, lamda_01):
"""
计算kz
输出数组
"""
v1 = ImageHandler().get_data(master_incident_file) # v1是cos值
# sin_v1=np.sqrt(1 - v1 * v1) # 夹角S1 正弦值
sin_v1=np.sin(np.arccos(v1))
v_b= ImageHandler().get_data(B_vertical_path) # v1是cos值
r1 = ImageHandler().get_data(r1_file)
aaaa=4 * math.pi * v_b
bbbb=lamda_01 * r1 * sin_v1
kz_array=aaaa/bbbb
# kz_array = (4 * math.pi * v_b) / (lamda_01 * r1 * sin_v1)
return kz_array
from xml.etree.ElementTree import ElementTree, Element
class header_file_read_angle():
"""从_incidence.xml文件中读取角度，计算雷达入射角。"""
def __init__(self):
pass
def read_all_angle(self,incidence_file,hight,width):
"""从_incidence.xml文件中读取角度，存到一个数组中"""
tree = ElementTree()
tree.parse(incidence_file) # 影像头文件
root = tree.getroot()
element_trees = list(root)
all_angle_array=np.zeros((hight,width),dtype=float)
a=0
for i in range(0, len(element_trees)):
if element_trees[i].tag=='incidenceValue':
a=i
break
for j in range(0, width):
all_angle_array[:,j]=element_trees[j+a].text
return all_angle_array
def get_orisim_angle(self, ori_sim_file, all_angle_array, header_height, header_width):
"""根据ori_sim.tif文件记录的行列号去取对应的雷达入射角
ori_sim_file:正射产品
all_angle_array：从incidence.xml读取的所有雷达入射角数组
header_height：轨道参数文件记录的高
header_width：轨道参数文件记录的宽
"""
height_array = ImageHandler.get_band_array(ori_sim_file, 1)
width_array = ImageHandler.get_band_array(ori_sim_file, 2)
img_height=ImageHandler.get_img_height(ori_sim_file) # 读取影像的高度
img_width= ImageHandler.get_img_width(ori_sim_file) # 读取影像的宽度
angle_array=np.zeros((img_height, img_width),dtype=float) # 存放角度的空矩阵
for i in range(0,img_height):
for j in range(0, img_width):
x=height_array[i,j] # 取出行数值
y=width_array[i,j] # 取出列数值
if x>0 and y>0:
if x<header_height and y<header_width:
angle_array[i, j]=all_angle_array[int(x), int(y)] # 根据行列号去取对应的雷达入射角
else:
angle_array[i, j] =0
else:
angle_array[i, j] = 0
return angle_array
# 数据测试模块
# if __name__=="__main__":
#
# # 0、 准备数据
# Mas_tiff_path = r"C:\Users\Administrator\Desktop\新建文件夹6\MasterSar" # 主影像头文件存放的文件夹路径 已知
# Aux_tiff_path = r"C:\Users\Administrator\Desktop\新建文件夹6\AuxiliarySar" # 辅影像头文件存放的文件夹路径 已知
# Mas_ori_sim = r"C:\Users\Administrator\Desktop\新建文件夹6\MasterSar\ori_sim_MasterSar_preprocessed.tif" # 主影像地面p点影像 已知
# Aux_ori_sim = r"C:\Users\Administrator\Desktop\新建文件夹6\AuxiliarySar\ori_sim_AuxiliarySar_preprocessed.tif" # 辅影像地面p点影像 已知
#
# #########
# # 1.构建主辅影像行列变换模型，生成新的影像(p点坐标):out_new_master_ori_sim
# out_new_master_ori_sim= r"C:\Users\Administrator\Desktop\新建文件夹6\new_master_ori_sim.tif" # 写入新的地面p点路径
# r_plsq, c_plsq2=fine_registration().least_sqrt(Mas_ori_sim, Aux_ori_sim)
# fine_registration().get_new_master_sim(r_plsq, c_plsq2, Mas_ori_sim, out_new_master_ori_sim)
# # 2.计算卫星的空间位置和p0点的位置
#
# Orthorectification_RD_object=Orthorectification_RD() # 使用生成的out_new_master_ori_sim来计算空间位置
# Mas_p0_array, Mas_sat_array, Mas_lamda= Orthorectification_RD_object.Orthorectification(Mas_tiff_path, Mas_ori_sim)
# Aux_p0_array, Aux_sat_array, Aux_lamda= Orthorectification_RD_object.Orthorectification(Aux_tiff_path, out_new_master_ori_sim)
#
# S1=Mas_sat_array # 主影像卫星位置
# S2=Aux_sat_array # 辅影像卫星位置
# P=Mas_p0_array # 主影像计算出的p0
# # 3.计算垂直基线数组B_vertical、斜距len_S1_P、斜距len_S2_P
# B_vertical, len_S1_P, len_S2_P=fine_registration().calu_vector_angle(S1, S2, P)
# print("垂直基线计算完成")
# ########
#
# # 开始计算kz
# geo = ImageHandler.get_geotransform(Mas_ori_sim)
# pro = ImageHandler.get_projection(Mas_ori_sim)
# B_vertical_path=r"C:\Users\Administrator\Desktop\新建文件夹6\B_vertical_path.tif"
# ImageHandler.write_img(B_vertical_path, pro, geo, B_vertical, no_data='null')
#
# len_S1_P_path=r"C:\Users\Administrator\Desktop\新建文件夹6\len_S1_P_path.tif"
# ImageHandler.write_img(len_S1_P_path, pro, geo, len_S1_P, no_data='null')
# len_S2_P_path=r"C:\Users\Administrator\Desktop\新建文件夹6\len_S2_P_path.tif"
# ImageHandler.write_img(len_S2_P_path, pro, geo, len_S2_P, no_data='null')
#
#
# master_incident_file=r"C:\Users\Administrator\Desktop\新建文件夹6\MasterSar\IncidenceAngle.tif"
# Mas_lamda=0.055517
# kz_array=fine_registration().calu_kz(master_incident_file, B_vertical_path,len_S1_P_path, Mas_lamda)
# kz_array_path=r"C:\Users\Administrator\Desktop\新建文件夹6\kz_array_path.tif"
# ImageHandler.write_img(len_S2_P_path, pro, geo, kz_array, no_data='null')
# pass
# m_incidence_file=r"I:\西藏未正射\西藏\GF3_MYN_QPSI_011437_E98.7_N31.1_20181012_L1A_AHV_L10003514923\GF3_MYN_QPSI_011437_E98.7_N31.1_20181012_L1A_AHV_L10003514923.incidence.xml"
# m_ori_sim_file =r"D:\MicroWorkspace\L-SAR\VegetationHeight\Input\MasterSar\ori_sim.tif"
# m_all_angle_array=header_file_read_angle().read_all_angle(m_incidence_file,6000,8062)
# m_angle_array=header_file_read_angle().get_orisim_angle(m_ori_sim_file, m_all_angle_array)
#
#
# a_incidence_file = r"I:\西藏未正射\西藏\GF3_MYN_QPSI_011437_E98.7_N31.3_20181012_L1A_AHV_L10003514924\GF3_MYN_QPSI_011437_E98.7_N31.3_20181012_L1A_AHV_L10003514924.incidence.xml"
# a_ori_sim_file =r"D:\MicroWorkspace\L-SAR\VegetationHeight\Input\AuxiliarySar\ori_sim2.tif"
# a_all_angle_array = header_file_read_angle().read_all_angle(a_incidence_file, 6000, 8062)
# a_angle_array = header_file_read_angle().get_orisim_angle(a_ori_sim_file, a_all_angle_array)
#
# v1=m_angle_array[0:561,0:1262]
# v2=a_angle_array[0:561,0:1262]
# v1=v1/180*math.pi
# v2=v2/180 * math.pi
#
# mean = (v1 + v2) / 2
# chazhi = v1 - v2
# kz_array = 4 * math.pi * (abs(v2 - v1)) / (np.sin(mean) * 0.055517)
# pass