KJ135
/
ISCE_INSAR
mirror of http://172.16.0.12:5000/WBJY/ISCE_INSAR.git
4
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
ISCE_INSAR/components/isceobj/Sensor/LT1ABLT1ABREPEAT.py

1864 lines
71 KiB
Python
Raw Permalink Normal View History

LT1AB重复轨道适配
2023-12-06 09:44:22 +00:00
#!/usr/bin/env python3
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Copyright 2010 California Institute of Technology. ALL RIGHTS RESERVED.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# LT1AB 重复轨道模式
#
# Author: Chenzenghui
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
from xml.etree.ElementTree import ElementTree
import datetime
import isceobj
from isceobj.Util import Poly1D
from isceobj.Scene.Frame import Frame
from isceobj.Planet.Planet import Planet
from isceobj.Orbit.Orbit import StateVector, Orbit
from isceobj.Orbit.OrbitExtender import OrbitExtender
from isceobj.Planet.AstronomicalHandbook import Const
from iscesys.Component.Component import Component
from iscesys.DateTimeUtil.DateTimeUtil import DateTimeUtil as DTUtil
import os
import numpy as np
import math
import datetime
import time
from math import sin,cos
from scipy.optimize import leastsq
import numpy as np
from isceobj.Util.Poly2D import Poly2D
sep = "\n"
tab = " "
lookMap = { 'RIGHT' : -1,
'LEFT' : 1}
TIFF = Component.Parameter(
'tiff',
public_name='TIFF',
default='',
type=str,
mandatory=True,
doc='RadarSAT2 tiff imagery file'
)
XML = Component.Parameter(
'xml',
public_name='XML',
default='',
type=str,
mandatory=True,
doc='RadarSAT2 xml metadata file'
)
ORBIT_DIRECTORY = Component.Parameter(
'orbitDirectory',
public_name = 'orbit directory',
default=None,
type=str,
mandatory=False,
doc='Directory with Radarsat2 precise orbits')
ORBIT_FILE = Component.Parameter(
'orbitFile',
public_name = 'orbit file',
default = None,
type = str,
mandatory = False,
doc = 'Precise orbit file to use')
from .Sensor import Sensor
#######################################################################################################
# 增加轨道求解模型，为了方便轨道模型的求解，这里对轨道模型进行修改
# 这里采用《Insar原理和应用 》（刘国祥著）中基于空间定位的轨道参数计算方法一章中的内容（Gabriel and Goldstein）
# 注意为了方便计算，采用基准时间的方法
######################################################
def FindInfomationFromJson(HeaderFile_dom_json, node_path_list):
"""
在Json文件中，按照指定路径解析出制定节点
"""
result_node = HeaderFile_dom_json
for nodename in node_path_list:
result_node = result_node[nodename]
return result_node
def GetVectorNorm(Vecter):
"""
得到向量的模
"""
Vecter = Vecter.reshape(-1,1)
Vecter_Norm_pow = np.matmul(Vecter.T,Vecter)
return np.sqrt(Vecter_Norm_pow)
def XYZOuterM2(A, B):
"""
外积（叉乘）,日后版本换成可以任意维度的外积运算方程
args:
A:nx3
B:nx3
"""
cnt = A.shape[0]
C = np.zeros((cnt, 3))
C[:, 0] = A[:, 1] * B[:, 2] - A[:, 2] * B[:, 1]
C[:, 1] = A[:, 2] * B[:, 0] - A[:, 0] * B[:, 2]
C[:, 2] = A[:, 0] * B[:, 1] - A[:, 1] * B[:, 0]
return C
class SatelliteOrbit(object):
def __init__(self) -> None:
super().__init__()
self.starttime = 1262275200.0
self.modelName=""
def get_starttime(self):
'''
返回卫星轨道时间起算点
'''
return self.starttime
def ReconstructionSatelliteOrbit(self, GPSPoints_list):
'''
重建卫星轨道，使用多项式拟合法
args:
GPSPoints_list:GPS 卫星轨道点
return：
SatelliteOrbitModel 卫星轨道模型
'''
self.SatelliteOrbitModel = None
def SatelliteSpaceState(self, time_float):
'''
根据时间戳，返回对应时间的卫星的轨迹状态
args:
time_float:时间戳
return：
State_list:[time,Xp,Yp,Zp,Vx,Vy,Vz]
'''
return None
class SatelliteOrbitFitPoly(SatelliteOrbit):
'''
继承于SatelliteOribit类，为拟合多项式实现方法
'''
def __init__(self) -> None:
super().__init__()
self.modelName="多项式"
self.polynum=4
def ReconstructionSatelliteOrbit(self, GPSPoints_list, starttime):
if len(GPSPoints_list)==2:
self.polynum=1
self.starttime = starttime
record_count = len(GPSPoints_list)
time_arr = np.zeros((record_count, 1), dtype=np.float64) # 使用np.float64只是为了精度高些；如果32位也能满足需求，请用32位
state_arr = np.zeros((record_count, 6), dtype=np.float64)
A_arr = np.zeros((self.polynum+1, 6), dtype=np.float64) # 四次项
X=np.ones((record_count,self.polynum+1),dtype=np.float64) # 记录时间坐标
# 将点记录转换为自变量矩阵、因变量矩阵
for i in range(record_count):
GPSPoint = GPSPoints_list[i]
time_ = GPSPoint[0] - self.starttime # 为了保证精度，对时间进行缩放
X[i,:]=np.array([1,time_])
state_arr[i, :] = np.array(GPSPoint[1:],dtype=np.float64).reshape(1,6) # 空间坐标
self.model_f=[]
for i in range(6):
Y = state_arr[:, i].reshape(-1,1)
A_arr[:,i]=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)[:,0]
self.A_arr=copy.deepcopy(A_arr.copy())
return self.A_arr
elif len(GPSPoints_list) > 6:
# 多项式的节点数，理论上是超过5个可以起算，这里为了精度选择10个点起算。
# 多项式 XA=Y ==> A=(X'X)^X'Y，其中 A 为待求系数，X为变量，Y为因变量
# 这里使用三次项多项式，共有6组参数。
# 声明自变量，因变量，系数矩阵
self.starttime = starttime
record_count = len(GPSPoints_list)
time_arr = np.zeros((record_count, 1), dtype=np.float64) # 使用np.float64只是为了精度高些；如果32位也能满足需求，请用32位
state_arr = np.zeros((record_count, 6), dtype=np.float64)
A_arr = np.zeros((self.polynum, 6), dtype=np.float64) # 四次项
X=np.ones((record_count,self.polynum),dtype=np.float64) # 记录时间坐标
# 将点记录转换为自变量矩阵、因变量矩阵
for i in range(record_count):
GPSPoint = GPSPoints_list[i]
time_ = GPSPoint[0] - self.starttime # 为了保证精度，对时间进行缩放
X[i,:]=np.array(list(map(lambda ii:time_**ii,range(self.polynum))))
state_arr[i, :] = np.array(GPSPoint[1:],dtype=np.float64).reshape(1,6) # 空间坐标
self.model_f=[]
for i in range(6):
Y = state_arr[:, i].reshape(-1,1)
A_arr[:,i]=np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),Y)[:,0]
self.A_arr=A_arr.copy()
return self.A_arr
else:
self.A_arr = None
return None
def SatelliteSpaceState(self, time_float):
'''
逐像素求解
根据时间戳，返回对应时间的卫星的轨迹状态，会自动计算与起算时间之差
args:
time_float:时间戳
return：
State_list:[time,Xp,Yp,Zp,Vx,Vy,Vz]
'''
if self.model_f is None:
return None
result_arr=np.zeros((1,7))
time_float = time_float - self.starttime
#
px=0
py=0
pz=0
vx=0
vy=0
vz=0
for ii in range(self.polynum):
px+=self.A_arr[ii,0]*time_float**ii
py+=self.A_arr[ii,1]*time_float**ii
pz+=self.A_arr[ii,2]*time_float**ii
vx+=self.A_arr[ii,3]*time_float**ii
vy+=self.A_arr[ii,4]*time_float**ii
vz+=self.A_arr[ii,5]*time_float**ii
return [time_float,[px,py,pz,vx,vy,vz]]
def getTimeOrbitStamp(self,UTCStartTime_float):
sv=_StateVector()
temp_sv=self.SatelliteSpaceState(UTCStartTime_float)
sv.timeStamp=datetime.datetime.fromtimestamp(UTCStartTime_float)
sv.xPosition = temp_sv[1][0,0]
sv.yPosition = temp_sv[1][0,1]
sv.zPosition = temp_sv[1][0,2]
sv.xVelocity = temp_sv[1][0,3]
sv.yVelocity = temp_sv[1][0,4]
sv.zVelocity = temp_sv[1][0,5]
return sv
def getTimeOrbits(self,UTCStartTime,UTCEndTime,orbitnum=1000):
#
startTime_stamp=datetime.datetime.strptime(UTCStartTime,"%Y-%m-%dT%H:%M:%S.%f").timestamp()-1
endTime_stamp=datetime.datetime.strptime(UTCEndTime,"%Y-%m-%dT%H:%M:%S.%f").timestamp()+1
if startTime_stamp>endTime_stamp:
raise
delta_t=(endTime_stamp-startTime_stamp)/orbitnum
extractOrbits=[]
#
temptime=startTime_stamp
while temptime<endTime_stamp:
temptime=temptime+delta_t
newOrbit=self.getTimeOrbitStamp(temptime)
extractOrbits.append(newOrbit)
newOrbit=self.getTimeOrbitStamp(endTime_stamp)
extractOrbits.append(newOrbit)
return extractOrbits # 扩展的轨道节点
def ReconstructionSatelliteOrbit(stateVectors, starttime):
'''
构建卫星轨道
args:
GPSPoints_list:卫星轨道点
starttime:起算时间
'''
GPSPoint_list=[]
for sv in stateVectors:
GPSPoint=[sv.timeStamp.timestamp(),
sv.xPosition,
sv.yPosition,
sv.zPosition,
sv.xVelocity,
sv.yVelocity,
sv.zVelocity
]
GPSPoint_list.append(GPSPoint)
SatelliteOrbitModel = SatelliteOrbitFitPoly()
if SatelliteOrbitModel.ReconstructionSatelliteOrbit(GPSPoint_list, starttime=starttime) is None:
return None
print("orbit test")
distance=[]
for gpsPoint in GPSPoint_list:
temp_sv=SatelliteOrbitModel.SatelliteSpaceState(gpsPoint[0])
sv=np.array(temp_sv[1])
temp_distance=sv-gpsPoint[1:]
distance.append(temp_distance)
distance=np.array(distance)
print("orbit max:",np.max(distance))
print("orbit min:",np.min(distance))
return SatelliteOrbitModel
########
# 函数列表
########
def poly1dfunc(p,x): # 一次函数
k,b=p
return k*x+b
def poly2dfunc(p,x): # 二次函数
k1,k2,b=p
return b+k1*x+k2*x*x
def poly3dfunc(p,x):
k1,k2,k3,b=p
return b+k1*x+k2*x*x+k3*x*x*x
def poly4dfunc(p,x):
k1,k2,k3,k4,b=p
return b+k1*x+k2*x**2+k3*x**3+k4*x**4
def poly5dfunc(p,x):
k1,k2,k3,k4,k5,b=p
return b+k1*x+k2*x**2+k3*x**3+k4*x**4+k5*x**5
def poly1derror(p,x,y):
return poly1dfunc(p,x)-y
def poly2derror(p,x,y):
return poly2dfunc(p,x)-y
def poly3derror(p,x,y):
return poly3dfunc(p,x)-y
def poly4derror(p,x,y):
return poly4dfunc(p,x)-y
def poly5derror(p,x,y):
return poly5dfunc(p,x)-y
class orbitVector:
def __init__(self,UTCTimes,vx,vy,vz,px,py,pz,dateformat="%Y-%m-%dT%H:%M:%S.%fZ"):
self.UTCTime=datetime.datetime.strptime(UTCTimes,dateformat) # 字符串转UTC时间
self.time_stamp=self.UTCTime.timestamp() # 时间戳
self.vx=vx # Vx
self.vy=vy
self.vz=vz
self.px=px
self.py=py # Ys
self.pz=pz
self.orbitVector2GG() # 将坐标进行变换
self.Check()
pass
def orbitVector2GG(self):
Cx=self.py*self.vz-self.pz*self.vy
Cy=self.pz*self.vx-self.px*self.vz
Cz=self.px*self.vy-self.py*self.vx
C=(Cx**2+Cy**2+Cz**2)**0.5
self.rho=(self.px**2+self.py**2+self.pz**2)**0.5
self.I=math.acos(Cz/C)
self.omega=math.asin(self.pz/(self.rho*math.sin(self.I)))
self.Omega=math.atan(-Cx/Cy)
return [self.rho,self.I,self.omega,self.Omega]
def Check(self):
rho=self.rho
I=self.I
omega=self.omega
Omega=self.Omega
Xp=rho*(cos(omega)*cos(Omega)-sin(omega)*sin(Omega)*cos(I))
Yp=rho*(cos(omega)*sin(Omega)+sin(omega)*cos(Omega)*cos(I))
Zp=rho*sin(Omega)*sin(I)
print("计算插值:",self.UTCTime,self.px,self.py,self.pz,self.vx,self.vy,self.vz,"|",Xp-self.px,Yp-self.py,Zp-self.pz)
class SARorbit(object):
# 作为自定义轨道计算基类
# 定义几个基本方法：addvector(),createOrbit(),getTime(),getTimes()
# 注意均为UTC
def __init__(self):
self.vectors=[]
self.baseTime=None
def addvector(self,UTCTime,vx,vy,vz,px,py,pz):
neworbit_point=orbitVector(UTCTime,vx,vy,vz,px,py,pz)
self.vectors.append(neworbit_point)
if len(self.vectors)==1:
self.baseTime=neworbit_point.time_stamp-1 # 基线
else:
if self.baseTime>neworbit_point.time_stamp:
self.baseTime=neworbit_point.time_stamp-1
def createOrbit(self):
pass
def getTimeOrbit(self,UTCTime):
return None
def getTimeOrbitStamp(self,StampTime):
utcStr=datetime.datetime.fromtimestamp(StampTime).strftime("%Y-%m-%dT%H:%M:%S.%fZ")
return self.getTimeOrbit(utcStr)
def getTimeOrbits(self,UTCStartTime,UTCEndTime,orbitnum=100):
#
startTime_stamp=datetime.datetime.strptime(UTCStartTime,"%Y-%m-%dT%H:%M:%S.%fZ").timestamp()-10
endTime_stamp=datetime.datetime.strptime(UTCEndTime,"%Y-%m-%dT%H:%M:%S.%fZ").timestamp()+10
if startTime_stamp>endTime_stamp:
raise
delta_t=(endTime_stamp-startTime_stamp)*1000/orbitnum
delta_t=int(delta_t)/1000 # 获取
extractOrbits=[]
#
temptime=startTime_stamp
while temptime<endTime_stamp:
temptime=temptime+delta_t
newOrbit=self.getTimeOrbitStamp(temptime)
extractOrbits.append(newOrbit)
newOrbit=self.getTimeOrbitStamp(endTime_stamp)
extractOrbits.append(newOrbit)
return extractOrbits # 扩展的轨道节点
pass
class polyorbit(SARorbit):
def __init__(self):
super(polyorbit,self).__init__()
def createOrbit(self):
if len(self.vectors)<8:
raise
ts,xs,ys,zs,vx,vy,vz=[],[],[],[],[],[],[]
for orbit_point in self.vectors:
ts.append(orbit_point.time_stamp-self.baseTime)
xs.append(orbit_point.px)
ys.append(orbit_point.py)
zs.append(orbit_point.pz)
vx.append(orbit_point.vx)
vy.append(orbit_point.vy)
vz.append(orbit_point.vz)
ts=np.array(ts)
xs=np.array(xs)
ys=np.array(ys)
zs=np.array(zs)
vx=np.array(vx)
vy=np.array(vy)
vz=np.array(vz)
self.p_xs=leastsq(poly5derror,[100,100,100,100,100,100],args=(ts,xs))
self.p_ys=leastsq(poly5derror,[100,100,100,100,100,100],args=(ts,ys))
self.p_zs=leastsq(poly5derror,[100,100,100,100,100,100],args=(ts,zs))
self.p_vxs=leastsq(poly5derror,[100,100,100,100,100,100],args=(ts,vx))
self.p_vy=leastsq(poly5derror,[100,100,100,100,100,100],args=(ts,vy))
self.p_vz=leastsq(poly5derror,[100,100,100,100,100,100],args=(ts,vz))
return
def getTimeOrbit(self,UTCTime):
staptime=datetime.datetime.strptime(UTCTime,"%Y-%m-%dT%H:%M:%S.%fZ").timestamp()-self.baseTime
xs=poly5dfunc(self.p_xs,staptime)
ys=poly5dfunc(self.p_ys,staptime)
zs=poly5dfunc(self.p_zs,staptime)
vx=poly5dfunc(self.p_vx,staptime)
vy=poly5dfunc(self.p_vy,staptime)
vz=poly5dfunc(self.p_vz,staptime)
utmstr=datetime.datetime.strptime(UTCTime,"%Y-%m-%dT%H:%M:%S.%fZ").strftime("%Y-%m-%d %H:%M:%S.%f")
neworbit_point=orbitVector(utmstr,vx,vy,vz,xs,ys,zs,dateformat="%Y-%m-%d %H:%M:%S.%f")
return neworbit_point
#######################################################################################################
class LT1ABLT1ABREPEAT(Sensor):
"""
A Class representing RADARSAT 2 data
"""
family='LT1ABLT1ABREPEAT'
parameter_list = (XML, TIFF, ORBIT_DIRECTORY, ORBIT_FILE ) + Sensor.parameter_list
def __init__(self, family='', name=''):
super().__init__(family if family else self.__class__.family, name=name)
self.product = _Product()
self.frame = Frame()
self.frame.configure()
def getFrame(self):
return self.frame
def parse(self):
try:
fp = open(self.xml,'r')
except IOError as strerr:
print("IOError: %s" % strerr)
return
self._xml_root = ElementTree(file=fp).getroot()
self.product.set_from_etnode(self._xml_root) # 解析meta xml
self.populateMetadata()
self.extractDoppler()
fp.close()
def populateMetadata(self):
"""
Create metadata objects from the metadata files
"""
mission = self.product.mission
swath = self.product.swath
frequency = self.product.frequency
orig_prf =self.product.orig_prf
rangePixelSize = self.product.rangePixelSize
rangeSamplingRate = Const.c/(2*rangePixelSize)
pulseLength =self.product.pulseLength
pulseBandwidth = self.product.pulseBandwidth
polarization = self.product.polarization
lookSide = self.product.lookSide
facility = self.product.facility
version = self.product.version
lines = self.product.lines
samples = self.product.samples
startingRange = self.product.startingRange
incidenceAngle = self.product.incidenceAngle
# some RS2 scenes have oversampled SLC images because processed azimuth bandwidth larger than PRF EJF 2015/08/15
azimuthPixelSize = self.product.azimuthPixelSize # ground spacing in meters
totalProcessedAzimuthBandwidth = self.product.totalProcessedAzimuthBandwidth
prf = orig_prf * np.ceil(totalProcessedAzimuthBandwidth / orig_prf) # effective PRF can be double original, suggested by Piyush
print('-------------------------------------')
print("mission",mission)
print("swath",swath)
print("frequency",frequency)
print("orig_prf",orig_prf)
print("prf",prf)
print("rangePixelSize",rangePixelSize)
print("rangeSamplingRate",rangeSamplingRate)
print("pulseLength",pulseLength)
print("pulseBandwidth",pulseBandwidth)
print("polarization",polarization)
print("lookSide",lookSide)
print("lines",lines)
print("samples",samples)
print("startingRange",startingRange)
print("incidenceAngle",incidenceAngle)
print(self.product.imageGenerationParameters.dopplerCentroid.dopplerCentroidCoefficients)
print(self.product.imageGenerationParameters.dopplerRateValues.dopplerRateValuesCoefficients)
print("effective PRF %f, original PRF %f" % (prf, orig_prf) )
print('-------------------------------------------')
#lineFlip = (self.product.imageAttributes.rasterAttributes.lineTimeOrdering.upper() == 'DECREASING')
dataStartTime = self.product.dataStartTime
dataStopTime = self.product.dataStopTime
passDirection = self.product.passDirection
####Populate platform
platform = self.frame.getInstrument().getPlatform()
platform.setPlanet(Planet(pname="Earth"))
platform.setMission(mission)
platform.setPointingDirection(lookSide)
platform.setAntennaLength(15.0)
####Populate instrument
instrument = self.frame.getInstrument()
instrument.setRadarFrequency(frequency)
instrument.setPulseRepetitionFrequency(prf)
instrument.setPulseLength(pulseLength)
instrument.setChirpSlope(pulseBandwidth/pulseLength)
instrument.setIncidenceAngle(incidenceAngle)
#self.frame.getInstrument().setRangeBias(0)
instrument.setRangePixelSize(rangePixelSize)
instrument.setRangeSamplingRate(rangeSamplingRate)
instrument.setBeamNumber(swath)
instrument.setPulseLength(pulseLength)
#Populate Frame
#self.frame.setSatelliteHeight(height)
self.frame.setSensingStart(dataStartTime)
self.frame.setSensingStop(dataStopTime)
diffTime = DTUtil.timeDeltaToSeconds(dataStopTime - dataStartTime)/2.0
sensingMid = dataStartTime + datetime.timedelta(microseconds=int(diffTime*1e6))
self.frame.setSensingMid(sensingMid)
self.frame.setPassDirection(passDirection)
self.frame.setPolarization(polarization)
self.frame.setStartingRange(startingRange)
self.frame.setFarRange(startingRange + (samples-1)*rangePixelSize)
self.frame.setNumberOfLines(lines)
self.frame.setNumberOfSamples(samples)
self.frame.setProcessingFacility("LAMP") # 处理软件
self.frame.setProcessingSoftwareVersion(version)
self.frame.setPassDirection(passDirection)
self.frame.setProcessingSystem("LT1AB") # 陆探数据
self.frame.getOrbit().setOrbitSource("LT1AB.META.XML")
print("======= 修改轨道代码部分： =================")
if (self.orbitFile is None) and (self.orbitDirectory is None):
self.extractOrbit()
elif (self.orbitDirectory is not None):
self.orbitFile = findPreciseOrbit(self.orbitDirectory, self.frame.getOrbit().getOrbitSource(), self.frame.sensingStart.year)
if self.orbitFile is not None:
self.extractPreciseOrbit(self.orbitFile, self.frame.sensingStart, self.frame.sensingStop)
# save the Doppler centroid coefficients, converting units from product.xml file
# units in the file are quadratic coefficients in Hz, Hz/sec, and Hz/(sec^2)
# ISCE expects Hz, Hz/(range sample), Hz((range sample)^2
# note that RS2 Doppler values are estimated at time dc.dopplerCentroidReferenceTime,
# so the values might need to be adjusted for ISCE usage
# added EJF 2015/08/17
dc = self.product.dopplerCentroid
poly = dc.dopplerCentroidCoefficients
# need to convert units
poly[1] = poly[1]/rangeSamplingRate
poly[2] = poly[2]/rangeSamplingRate**2
self.doppler_coeff = poly
# similarly save Doppler azimuth fm rate values, converting units
# units in the file are quadratic coefficients in Hz, Hz/sec, and Hz/(sec^2)
# Guessing that ISCE expects Hz, Hz/(range sample), Hz((range sample)^2
# note that RS2 Doppler values are estimated at time dc.dopplerRateReferenceTime,
# so the values might need to be adjusted for ISCE usage
# added EJF 2015/08/17
dr = self.product.dopplerRateValues
fmpoly = dr.dopplerRateValuesCoefficients
# need to convert units
fmpoly[1] = fmpoly[1]/rangeSamplingRate
fmpoly[2] = fmpoly[2]/rangeSamplingRate**2
self.azfmrate_coeff = fmpoly
# now calculate effective PRF from the azimuth line spacing after we have the orbit info EJF 2015/08/15
# this does not work because azimuth spacing is on ground. Instead use bandwidth ratio calculated above EJF
# SCHvelocity = self.frame.getSchVelocity()
# SCHvelocity = 7550.75 # hard code orbit velocity for now m/s
# prf = SCHvelocity/azimuthPixelSize
# instrument.setPulseRepetitionFrequency(prf)
def extractOrbit(self):
'''
Extract the orbit state vectors from the XML file.
'''
# Initialize orbit objects
增加HJ2星的代码适配
2025-02-05 01:37:37 +00:00
# Read into temp orbit first.frame
LT1AB重复轨道适配
2023-12-06 09:44:22 +00:00
# Radarsat 2 needs orbit extensions.
print("构建新的轨道程序代码：")
self.frame.getOrbit().setOrbitSource('Header: ' + self.frame.getOrbit().getOrbitSource())
stateVectors = self.product.OrbitModelInfo.stateVectors
tempOrbit= ReconstructionSatelliteOrbit(stateVectors,self.product.orbitstarttime.timestamp())
As_arr=tempOrbit.A_arr.tolist()
""" 构建轨道 """
self.frame.getOrbit().setsessionMode(self.product.mission)
增加HJ2星的代码适配
2025-02-05 01:37:37 +00:00
print("卫星型号")
print(self.frame.orbit.sessionMode)
LT1AB重复轨道适配
2023-12-06 09:44:22 +00:00
self.frame.getOrbit().setPolyParams(tempOrbit.polynum,tempOrbit.starttime,As_arr)
self.frame.getOrbit().setDoppler(self.product.dopplerRateValues.dopplerRateReferenceTime,
self.product.dopplerRateValues.dopplerRateValuesCoefficients)
""" """
newOrb=self.frame.getOrbit().getTimeOrbits(self.frame.sensingStart.strftime("%Y-%m-%dT%H:%M:%S.%f"),
self.frame.sensingStop.strftime("%Y-%m-%dT%H:%M:%S.%f"),
orbitnum=1000)
for svect in newOrb:
self.frame.getOrbit().addStateVector(svect)
print('Successfully read state vectors from product XML')
def extractPreciseOrbit(self, orbitfile, tstart, tend):
'''
Extract precise orbits for given time-period from orbit file.
'''
self.frame.getOrbit().setOrbitSource('File: ' + orbitfile)
tmin = tstart - datetime.timedelta(seconds=30.)
tmax = tstart + datetime.timedelta(seconds=30.)
fid = open(orbitfile, 'r')
for line in fid:
if not line.startswith('; Position'):
continue
else:
break
for line in fid:
if not line.startswith(';###END'):
tstamp = convertRSTimeToDateTime(line)
if (tstamp >= tmin) and (tstamp <= tmax):
sv = StateVector()
sv.configure()
sv.setTime( tstamp)
sv.setPosition( [float(x) for x in fid.readline().split()])
sv.setVelocity( [float(x) for x in fid.readline().split()])
self.frame.getOrbit().addStateVector(sv)
else:
fid.readline()
fid.readline()
dummy = fid.readline()
if not dummy.startswith(';'):
raise Exception('Expected line to start with ";". Got {0}'.format(dummy))
fid.close()
print('Successfully read {0} state vectors from {1}'.format( len(self.frame.getOrbit()._stateVectors), orbitfile))
def extractImage(self, verbose=True):
'''
Use gdal to extract the slc.
'''
try:
from osgeo import gdal
except ImportError:
raise Exception('GDAL python bindings not found. Need this for RSAT2 / TandemX / Sentinel1A.')
self.parse()
width = self.frame.getNumberOfSamples()
lgth = self.frame.getNumberOfLines()
lineFlip = False #(self.product.imageAttributes.rasterAttributes.lineTimeOrdering.upper() == 'DECREASING')
pixFlip = False #(self.product.imageAttributes.rasterAttributes.pixelTimeOrdering.upper() == 'DECREASING')
src = gdal.Open(self.tiff.strip(), gdal.GA_ReadOnly)
cJ = np.complex64(1.0j)
####Images are small enough that we can do it all in one go - Piyush
real = src.GetRasterBand(1).ReadAsArray(0,0,width,lgth)
imag = src.GetRasterBand(2).ReadAsArray(0,0,width,lgth)
if (real is None) or (imag is None):
raise Exception('Input Radarsat2 SLC seems to not be a 2 band Int16 image.')
data = real+cJ*imag
real = None
imag = None
src = None
if lineFlip:
if verbose:
print('Vertically Flipping data')
data = np.flipud(data)
if pixFlip:
if verbose:
print('Horizontally Flipping data')
data = np.fliplr(data)
data.tofile(self.output)
####
slcImage = isceobj.createSlcImage()
slcImage.setByteOrder('l')
slcImage.setFilename(self.output)
slcImage.setAccessMode('read')
slcImage.setWidth(width)
slcImage.setLength(lgth)
slcImage.setXmin(0)
slcImage.setXmax(width)
#slcImage.renderHdr()
self.frame.setImage(slcImage)
def extractDoppler(self):
'''
self.parse()
Extract doppler information as needed by mocomp
'''
ins = self.frame.getInstrument()
dc = self.product.dopplerCentroid
quadratic = {}
r0 = self.frame.startingRange
fs = ins.getRangeSamplingRate()
tNear = 2*r0/Const.c
tMid = tNear + 0.5*self.frame.getNumberOfSamples()/fs
t0 = dc.dopplerCentroidReferenceTime
poly = dc.dopplerCentroidCoefficients
fd_mid = 0.0
for kk in range(len(poly)):
fd_mid += poly[kk] * (tMid - t0)**kk
####For insarApp
quadratic['a'] = fd_mid / ins.getPulseRepetitionFrequency()
quadratic['b'] = 0.
quadratic['c'] = 0.
####For roiApp
####More accurate
coeffs = poly
dr = self.frame.getInstrument().getRangePixelSize()
rref = 0.5 * Const.c * t0
r0 = self.frame.getStartingRange()
norm = 0.5*Const.c/dr
dcoeffs = []
for ind, val in enumerate(coeffs):
dcoeffs.append( val / (norm**ind))
poly = Poly1D.Poly1D()
poly.initPoly(order=len(coeffs)-1)
poly.setMean( (rref - r0)/dr - 1.0)
poly.setCoeffs(dcoeffs)
pix = np.linspace(0, self.frame.getNumberOfSamples(), num=len(coeffs)+1)
evals = poly(pix)
fit = np.polyfit(pix,evals, len(coeffs)-1)
self.frame._dopplerVsPixel = list(fit[::-1])
self.frame._dopplerVsPixel_LT1AB = list([self.product.dopplerRateValues.dopplerRateReferenceTime]
+self.product.dopplerRateValues.dopplerRateValuesCoefficients)
print('Doppler Fit: ', fit[::-1])
# print("---- radar ------------------------------------------print")
# print("****" * 10)
# print("t0",t0)
# print('lightspeed',Const.c)
# print('rref',rref)
# print('dr',dr)
# print('dcoeff',dcoeffs)
# print('r0',r0)
# print("****" * 10)
# print('pix',pix)
# print('evals',evals)
# print('fit',fit)
# print('Doppler Fit: ', fit[::-1])
# print('---------------------------------------------------')
# """ 测试 多普勒"""
# print("======= Test LTInSAR =========")
# """ 测试 轨道情况"""
# coeff = self.frame._dopplerVsPixel_LT1AB
# doppler = Poly2D()
# doppler._meanRange = self.frame.startingRange
# doppler._normRange = self.frame.instrument.rangePixelSize
# doppler.initPoly(azimuthOrder=0, rangeOrder=len(coeff)-1, coeffs=[coeff])
# lookSide = self.frame.instrument.platform.pointingDirection
# planet = Planet(pname='Earth')
# wvl = self.frame.instrument.getRadarWavelength()
# test_lla=[47.75,130.82,159.50]
# taz, rgm = self.frame.orbit.geo2rdr(test_lla, side=lookSide, doppler=doppler, wvl=wvl)
# print("taz",taz)
# print('n11',datetime.datetime.strptime("2023-03-27T21:28:27.551796","%Y-%m-%dT%H:%M:%S.%f"))
# print("rgm",rgm)
# print('r11',4.96786423292669768E-03*Const.c/2)
# print("lon,lat,ati",test_lla)
# print("==============================================")
return quadratic
class LT1ABNamespace(object):
def __init__(self):
self.uri = ""
def elementName(self,element):
return element
def convertToDateTime(self,string):
dt = datetime.datetime.strptime(string,"%Y-%m-%dT%H:%M:%S.%f")
return dt
class _Product(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.productId = None
self.documentId = None
self.mission=None
self.swath=None
self.frequency=None
self.orig_prf=None # 原始PRF --- 参考PRF
self.rangePixelSize=None # 距离向分辨率
self.pulseLength=None
self.pulseBandwidth=None
self.polarization=None
self.lookSide=None
self.facility=None
self.version=None
self.lines=None # 行数
self.samples=None # 列数
self.startingRange=None
self.incidenceAngle=None
self.azimuthPixelSize=None
self.totalProcessedAzimuthBandwidth=None
self.prf=None
"""
入射角
"""
self.topRight_IncidenceAngle=None
self.topLeft_IncidenceAngle=None
self.bottomRight_IncidenceAngle=None
self.bottomLeft_IncidenceAngle=None
self.Near_incidenceAngle=None
self.Far_incidenceAngle=None
""" 斜距 """
self.slantRangeTimeToFirstRangeSample=None
self.startingRange=None
""" 轨道"""
self.passDirection=None
""" 成像 """
self.dataStartTime=None
self.dataStopTime=None
self.orbitstarttime=None
""" 多普勒系数 """
self.dopplerCentroid=_DopplerCentroid()
self.dopplerRateValues = _DopplerRateValues()
self.OrbitModelInfo=_OrbitInformation()
self.sourceAttributes = _SourceAttributes()
self.imageGenerationParameters = _ImageGenerationParameters()
self.imageAttributes = _ImageAttributes()
def set_from_etnode(self,node):
self.productId=node.find("productInfo").find("generationInfo").find("logicalProductID").text
self.documentId=node.find("generalHeader").find("referenceDocument").text
self.mission=node.find("generalHeader").find("mission").text
self.swath=node.find("productInfo").find("acquisitionInfo").find("imagingMode").text
self.frequency=float(node.find("processing").find("processingParameter").find("rangeCompression").find("chirps").find("referenceChirp").find("centerFrequency").text)
self.rangePixelSize=float(node.find("productInfo").find("imageDataInfo").find("imageRaster").find("columnSpacing").text)
self.pulseLength=float(node.find("processing").find("processingParameter").find("rangeCompression").find("chirps").find("referenceChirp").find("pulseLength").text)
self.pulseBandwidth=float(node.find("processing").find("processingParameter").find("rangeCompression").find("chirps").find("referenceChirp").find("pulseBandwidth").text)
self.polarization=node.find("productInfo").find("acquisitionInfo").find("polarisationMode").text
self.lookSide=lookMap[node.find("productInfo").find("acquisitionInfo").find("lookDirection").text.upper()]
self.facility=node.find("productInfo").find("generationInfo").find("level1ProcessingFacility").text
self.version="V1.1"
self.lines=int(node.find("productInfo").find("imageDataInfo").find("imageRaster").find("numberOfRows").text)
self.samples=int(node.find("productInfo").find("imageDataInfo").find("imageRaster").find("numberOfColumns").text)
sceneCorners=node.find("productInfo").find("sceneInfo").findall("sceneCornerCoord")
for sceneCorner in sceneCorners:
incidenceAngle_temp=float(sceneCorner.find("incidenceAngle").text)
if sceneCorner.attrib["name"]=="topRight":
self.topRight_IncidenceAngle=incidenceAngle_temp
elif sceneCorner.attrib["name"]=="topLeft":
self.topLeft_IncidenceAngle=incidenceAngle_temp
elif sceneCorner.attrib["name"]=="bottomRight":
self.bottomRight_IncidenceAngle=incidenceAngle_temp
elif sceneCorner.attrib["name"]=="bottomLeft":
self.bottomLeft_IncidenceAngle=incidenceAngle_temp
self.Near_incidenceAngle=self.topRight_IncidenceAngle
self.Near_incidenceAngle=self.Near_incidenceAngle if self.Near_incidenceAngle < self.topLeft_IncidenceAngle else self.topLeft_IncidenceAngle
self.Near_incidenceAngle=self.Near_incidenceAngle if self.Near_incidenceAngle < self.bottomRight_IncidenceAngle else self.bottomRight_IncidenceAngle
self.Near_incidenceAngle=self.Near_incidenceAngle if self.Near_incidenceAngle < self.bottomLeft_IncidenceAngle else self.bottomLeft_IncidenceAngle
self.Far_incidenceAngle=self.topRight_IncidenceAngle
self.Far_incidenceAngle=self.Far_incidenceAngle if self.Far_incidenceAngle > self.topLeft_IncidenceAngle else self.topLeft_IncidenceAngle
self.Far_incidenceAngle=self.Far_incidenceAngle if self.Far_incidenceAngle > self.bottomRight_IncidenceAngle else self.bottomRight_IncidenceAngle
self.Far_incidenceAngle=self.Far_incidenceAngle if self.Far_incidenceAngle > self.bottomLeft_IncidenceAngle else self.bottomLeft_IncidenceAngle
self.incidenceAngle=float((node.find("productInfo").find("sceneInfo").find("sceneCenterCoord").find("incidenceAngle").text))
""" 斜距 """
self.slantRangeTimeToFirstRangeSample=float(node.find("productInfo").find("sceneInfo").find("rangeTime").find("firstPixel").text)
self.startingRange=self.slantRangeTimeToFirstRangeSample * (Const.c/2)
self.azimuthPixelSize=float(node.find("productInfo").find("imageDataInfo").find("imageRaster").find("rowSpacing").text)
self.totalProcessedAzimuthBandwidth=float(node.find("processing").find("processingParameter").find("totalProcessedAzimuthBandwidth").text) # 方位向
self.prf=float(node.find("instrument").find("settings").find("settingRecord").find("PRF").text)
""" 成像 """
self.dataStartTime=self.convertToDateTime(node.find("productInfo").find("sceneInfo").find("start").find("timeUTC").text)
self.dataStopTime=self.convertToDateTime(node.find("productInfo").find("sceneInfo").find("stop").find("timeUTC").text)
self.orig_prf=(self.lines-1)/((self.dataStopTime-self.dataStartTime).total_seconds())
for z in node.getchildren():
if z.tag == self.elementName('platform'):
self.OrbitModelInfo.set_from_etnode(z.find("orbit"))
elif z.tag == self.elementName('imageGenerationParameters'):
self.imageGenerationParameters.set_from_etnode(z)
elif z.tag == self.elementName('imageAttributes'):
self.imageAttributes.set_from_etnode(z)
""" 轨道 """
self.passDirection=node.find("productInfo").find("missionInfo").find("orbitDirection").text
self.orbitstarttime=(self.dataStopTime-self.dataStartTime)/2+self.dataStartTime
""" 多普勒系数 """
self.dopplerCentroid.set_from_etnode(node.find("processing").find("doppler").find("dopplerCentroid").find("dopplerEstimate").find("basebandDoppler"))
self.dopplerRateValues.set_from_etnode(node.find("processing").find("geometry").find("dopplerRate").find("dopplerRatePolynomial"))
def __str__(self):
retstr = "Product:"+sep+tab
retlst = ()
retstr += "productID=%s"+sep+tab
retlst += (self.productId,)
retstr += "documentIdentifier=%s"+sep
retlst += (self.documentId,)
retstr += "%s"+sep
retlst += (str(self.sourceAttributes),)
retstr += "%s"+sep
retlst += (str(self.imageGenerationParameters),)
retstr += ":Product"
return retstr % retlst
class _SourceAttributes(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.satellite = None
self.sensor = None
self.inputDatasetId = None
self.imageId = None
self.inputDatasetFacilityId = None
self.beamModeId = None
self.beamModeMnemonic = None
self.rawDataStartTime = None
self.radarParameters = _RadarParameters()
self.rawDataAttributes = _RawDataAttributes()
self.orbitAndAttitude = _OrbitAndAttitude()
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('satellite'):
self.satellite = z.text
elif z.tag == self.elementName('sensor'):
self.sensor = z.text
elif z.tag == self.elementName('inputDatasetId'):
self.inputDatasetId = z.text
elif z.tag == self.elementName('imageID'):
self.imageId = z.text
elif z.tag == self.elementName('inputDatasetFacilityId'):
self.inputDatasetFacilityId = z.text
elif z.tag == self.elementName('beamModeID'):
self.beamModeId = z.text
elif z.tag == self.elementName('beamModeMnemonic'):
self.beamModeMnemonic = z.text
elif z.tag == self.elementName('rawDataStartTime'):
self.rawDataStartTime = self.convertToDateTime(z.text)
elif z.tag == self.elementName('radarParameters'):
self.radarParameters.set_from_etnode(z)
elif z.tag == self.elementName('rawDataAttributes'):
self.rawDataAttributes.set_from_etnode(z)
elif z.tag == self.elementName('orbitAndAttitude'):
self.orbitAndAttitude.set_from_etnode(z)
def __str__(self):
retstr = "SourceAttributes:"+sep+tab
retlst = ()
retstr += "satellite=%s"+sep+tab
retlst += (self.satellite,)
retstr += "sensor=%s"+sep+tab
retlst += (self.sensor,)
retstr += "inputDatasetID=%s"+sep
retlst += (self.inputDatasetId,)
retstr += "%s"
retlst += (str(self.radarParameters),)
retstr += "%s"
retlst += (str(self.rawDataAttributes),)
retstr += "%s"
retlst += (str(self.orbitAndAttitude),)
retstr += ":SourceAttributes"
return retstr % retlst
class _RadarParameters(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.acquisitionType = None
self.beams = None
self.polarizations = None
self.pulses = None
self.rank = None
self.settableGains = []
self.radarCenterFrequency = None
self.prf = None
self.pulseLengths = []
self.pulseBandwidths = []
self.antennaPointing = None
self.adcSamplingRate = []
self.yawSteeringFlag = None
self.geodeticFlag = None
self.rawBitsPerSample = None
self.samplesPerEchoLine = None
self.referenceNoiseLevels = [_ReferenceNoiseLevel()]*3
def set_from_etnode(self,node):
i = 0
for z in node.getchildren():
if z.tag == self.elementName('acquisitionType'):
self.acquisitionType = z.text
elif z.tag == self.elementName('beams'):
self.beams = z.text
elif z.tag == self.elementName('polarizations'):
self.polarizations = z.text
elif z.tag == self.elementName('pulses'):
self.pulses = z.text
elif z.tag == self.elementName('rank'):
self.rank = z.text
elif z.tag == self.elementName('settableGain'):
self.settableGains.append(z.text)
elif z.tag == self.elementName('radarCenterFrequency'):
self.radarCenterFrequency = float(z.text)
elif z.tag == self.elementName('pulseRepetitionFrequency'):
self.prf = float(z.text)
elif z.tag == self.elementName('pulseLength'):
self.pulseLengths.append(float(z.text))
elif z.tag == self.elementName('pulseBandwidth'):
self.pulseBandwidths.append(float(z.text))
elif z.tag == self.elementName('antennaPointing'):
self.antennaPointing = z.text
elif z.tag == self.elementName('adcSamplingRate'):
self.adcSamplingRate.append(float(z.text))
elif z.tag == self.elementName('yawSteeringFlag'):
self.yawSteeringFlag = z.text
elif z.tag == self.elementName('rawBitsPerSample'):
self.rawBitsPerSample = int(z.text)
elif z.tag == self.elementName('samplesPerEchoLine'):
self.samplesPerEchoLine = int(z.text)
elif z.tag == self.elementName('referenceNoiseLevels'):
self.referenceNoiseLevels[i].set_from_etnode(z)
i += 1
def __str__(self):
retstr = "RadarParameters:"+sep+tab
retlst = ()
retstr += "acquisitionType=%s"+sep+tab
retlst += (self.acquisitionType,)
retstr += "beams=%s"+sep+tab
retlst += (self.beams,)
retstr += "polarizations=%s"+sep+tab
retlst += (self.polarizations,)
retstr += "pulses=%s"+sep+tab
retlst += (self.pulses,)
retstr += "rank=%s"+sep
retlst += (self.rank,)
retstr += ":RadarParameters"+sep
return retstr % retlst
class _ReferenceNoiseLevel(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.pixelFirstNoiseValue = None
self.stepSize = None
self.numberOfNoiseLevelValues = None
self.noiseLevelValues = []
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('pixelFirstNoiseValue'):
self.pixelFirstNoiseValue = int(z.text)
elif z.tag == self.elementName('stepSize'):
self.stepSize = int(z.text)
elif z.tag == self.elementName('numberOfNoiseLevelValues'):
self.numberOfNoiseLevelValues = int(z.text)
elif z.tag == self.elementName('noiseLevelValues'):
self.noiseLevelValues = list(map(float,z.text.split()))
def __str__(self):
retstr = "ReferenceNoiseLevel:"+sep+tab
retlst = ()
retstr += "pixelFirstNoiseValue=%s"+sep+tab
retlst += (self.pixelFirstNoiseValue,)
retstr += "stepSize=%s"+sep+tab
retlst += (self.stepSize,)
retstr += "numberOfNoiseLevelValues=%s"+sep+tab
retlst += (self.numberOfNoiseLevelValues,)
retstr += "noiseLevelValues=%s"+sep+tab
retlst += (self.noiseLevelValues,)
retstr += sep+":ReferenceNoiseLevel"
return retstr % retlst
class _RawDataAttributes(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.numberOfInputDataGaps = None
self.gapSize = None
self.numberOfMissingLines = None
self.rawDataAnalysis = [_RawDataAnalysis]*4
def set_from_etnode(self,node):
pass
def __str__(self):
return ""
class _RawDataAnalysis(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
def set_from_etnode(self,node):
pass
class _OrbitAndAttitude(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.orbitInformation = _OrbitInformation()
self.attitudeInformation = _AttitudeInformation()
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('orbitInformation'):
self.orbitInformation.set_from_etnode(z)
elif z.tag == self.elementName('attitudeInformation'):
self.attitudeInformation.set_from_etnode(z)
def __str__(self):
retstr = "OrbitAndAttitude:"+sep
retlst = ()
retstr += "%s"
retlst += (str(self.orbitInformation),)
retstr += "%s"
retlst += (str(self.attitudeInformation),)
retstr += ":OrbitAndAttitude"+sep
return retstr % retlst
class _OrbitInformation(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.passDirection = None
self.orbitDataSource = None
self.orbitDataFile = None
self.stateVectors = []
def set_from_etnode(self,node): # orbit
self.orbitDataSource=node.find("orbitHeader").find("sensor").text
for z in node.getchildren():
if z.tag == self.elementName('stateVec'):
sv = _StateVector()
sv.set_from_etnode(z)
self.stateVectors.append(sv)
def __str__(self):
retstr = "OrbitInformation:"+sep+tab
retlst = ()
retstr += "passDirection=%s"+sep+tab
retlst += (self.passDirection,)
retstr += "orbitDataSource=%s"+sep+tab
retlst += (self.orbitDataSource,)
retstr += "orbitDataFile=%s"+sep
retlst += (self.orbitDataFile,)
retstr += ":OrbitInformation"+sep
return retstr % retlst
class _StateVector(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.timeStamp = None
self.xPosition = None
self.yPosition = None
self.zPosition = None
self.xVelocity = None
self.yVelocity = None
self.zVelocity = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('timeUTC'):
self.timeStamp = self.convertToDateTime(z.text)
elif z.tag == self.elementName('posX'):
self.xPosition = float(z.text)
elif z.tag == self.elementName('posY'):
self.yPosition = float(z.text)
elif z.tag == self.elementName('posZ'):
self.zPosition = float(z.text)
elif z.tag == self.elementName('velX'):
self.xVelocity = float(z.text)
elif z.tag == self.elementName('velY'):
self.yVelocity = float(z.text)
elif z.tag == self.elementName('velZ'):
self.zVelocity = float(z.text)
def __str__(self):
retstr = "StateVector:"+sep+tab
retlst = ()
retstr += "timeStamp=%s"+sep+tab
retlst += (self.timeStamp,)
retstr += "xPosition=%s"+sep+tab
retlst += (self.xPosition,)
retstr += "yPosition=%s"+sep+tab
retlst += (self.yPosition,)
retstr += "zPosition=%s"+sep+tab
retlst += (self.zPosition,)
retstr += "xVelocity=%s"+sep+tab
retlst += (self.xVelocity,)
retstr += "yVelocity=%s"+sep+tab
retlst += (self.yVelocity,)
retstr += "zVelocity=%s"+sep+tab
retlst += (self.zVelocity,)
retstr += sep+":StateVector"
return retstr % retlst
class _AttitudeInformation(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.attitudeDataSource = None
self.attitudeOffsetApplied = None
self.attitudeAngles = []
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('attitudeDataSource'):
self.attitudeDataSource = z.text
elif z.tag == self.elementName('attitudeOffsetApplied'):
self.attitudeOffsetApplied = z.text
elif z.tag == self.elementName('attitudeAngles'):
aa = _AttitudeAngles()
aa.set_from_etnode(z)
self.attitudeAngles.append(aa)
def __str__(self):
retstr = "AttitudeInformation:"+sep+tab
retlst = ()
retstr += "attitudeDataSource=%s"+sep+tab
retlst += (self.attitudeDataSource,)
retstr += "attitudeOffsetApplied=%s"+sep+tab
retlst += (self.attitudeOffsetApplied,)
retstr += "%s"+sep+tab
retlst += (map(str,self.attitudeAngles),)
retstr += ":AttitudeInformation"+sep
return retstr % retlst
class _AttitudeAngles(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.timeStamp = None
self.yaw = None
self.roll = None
self.pitch = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('timeStamp'):
self.timeStamp = self.convertToDateTime(z.text)
elif z.tag == self.elementName('yaw'):
self.yaw = float(z.text)
elif z.tag == self.elementName('roll'):
self.roll = float(z.text)
elif z.tag == self.elementName('pitch'):
self.pitch = float(z.text)
def __str__(self):
retstr = "AttitudeAngles:"+sep+tab
retlst = ()
retstr += "timeStamp=%s"+sep+tab
retlst += (self.timeStamp,)
retstr += "yaw=%s"+sep+tab
retlst += (self.yaw,)
retstr += "roll=%s"+sep+tab
retlst += (self.roll,)
retstr += "pitch=%s"+sep+tab
retlst += (self.pitch,)
retstr += sep+":AttitudeAngles"
return retstr % retlst
class _ImageGenerationParameters(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.generalProcessingInformation = _GeneralProcessingInformation()
self.sarProcessingInformation = _SarProcessingInformation()
self.dopplerCentroid = _DopplerCentroid()
self.dopplerRateValues = _DopplerRateValues()
self.chirp = []
self.slantRangeToGroundRange = _SlantRangeToGroundRange()
self.payloadCharacteristicsFile = []
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('generalProcessingInformation'):
self.generalProcessingInformation.set_from_etnode(z)
elif z.tag == self.elementName('sarProcessingInformation'):
self.sarProcessingInformation.set_from_etnode(z)
elif z.tag == self.elementName('dopplerCentroid'):
self.dopplerCentroid.set_from_etnode(z)
elif z.tag == self.elementName('dopplerRateValues'):
self.dopplerRateValues.set_from_etnode(z)
elif z.tag == self.elementName('slantRangeToGroundRange'):
self.slantRangeToGroundRange.set_from_etnode(z)
def __str__(self):
retstr = "ImageGenerationParameters:"+sep
retlst = ()
retstr += "%s"
retlst += (str(self.generalProcessingInformation),)
retstr += "%s"
retlst += (str(self.sarProcessingInformation),)
retstr += "%s"
retlst += (str(self.dopplerCentroid),)
retstr += "%s"
retlst += (str(self.dopplerRateValues),)
retstr += ":ImageGenerationParameters"
return retstr % retlst
class _GeneralProcessingInformation(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.productType = None
self._processingFacility = None
self.processingTime = None
self.softwareVersion = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('productType'):
self.productType = z.text
elif z.tag == self.elementName('_processingFacility'):
self._processingFacility = z.text
elif z.tag == self.elementName('processingTime'):
self.processingTime = self.convertToDateTime(z.text)
elif z.tag == self.elementName('softwareVersion'):
self.softwareVersion = z.text
def __str__(self):
retstr = "GeneralProcessingInformation:"+sep+tab
retlst = ()
retstr += "productType=%s"+sep+tab
retlst += (self.productType,)
retstr += "_processingFacility=%s"+sep+tab
retlst += (self._processingFacility,)
retstr += "processingTime=%s"+sep+tab
retlst += (self.processingTime,)
retstr += "softwareVersion=%s"+sep
retlst += (self.softwareVersion,)
retstr += ":GeneralProcessingInformation"+sep
return retstr % retlst
class _SarProcessingInformation(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.lutApplied = None
self.elevationPatternCorrection = None
self.rangeSpreadingLossCorrection = None
self.pulseDependantGainCorrection = None
self.receiverSettableGain = None
self.rawDataCorrection = None
self.rangeReferenceFunctionSource = None
self.interPolarizationMatricesCorrection = None
self.zeroDopplerTimeFirstLine = None
self.zeroDopplerTimeLastLine = None
self.numberOfLinesProcessed = None
self.samplingWindowStartTimeFirstRawLine = None
self.samplingWindowStartTimeLastRawLine = None
self.numberOfSwstChanges = None
self.numberOfRangeLooks = None
self.rangeLookBandwidth = None
self.totalProcessedRangeBandwidth = None
self.numberOfAzimuthLooks = None
self.scalarLookWeights = None
self.azimuthLookBandwidth = None
self.totalProcessedAzimuthBandwidth = None
self.azimuthWindow = _Window('Azimuth')
self.rangeWindow = _Window('Range')
self.incidenceAngleNearRange = None
self.incidenceAngleFarRange = None
self.slantRangeNearEdge = None
self._satelliteHeight = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('lutApplied'):
self.lutApplied = z.text
elif z.tag == self.elementName('numberOfLinesProcessed'):
self.numberOfLinesProcessed = int(z.text)
elif z.tag == self.elementName('azimuthWindow'):
self.azimuthWindow.set_from_etnode(z)
elif z.tag == self.elementName('rangeWindow'):
self.rangeWindow.set_from_etnode(z)
elif z.tag == self.elementName('incidenceAngleNearRange'):
self.incidenceAngleNearRange = float(z.text)
elif z.tag == self.elementName('incidenceAngleFarRange'):
self.incidenceAngleFarRange = float(z.text)
elif z.tag == self.elementName('slantRangeNearEdge'):
self.slantRangeNearEdge = float(z.text)
elif z.tag == self.elementName('totalProcessedAzimuthBandwidth'):
self.totalProcessedAzimuthBandwidth = float(z.text)
elif z.tag == self.elementName('_satelliteHeight'):
self._satelliteHeight = float(z.text)
elif z.tag == self.elementName('zeroDopplerTimeFirstLine'):
self.zeroDopplerTimeFirstLine = self.convertToDateTime(z.text)
elif z.tag == self.elementName('zeroDopplerTimeLastLine'):
self.zeroDopplerTimeLastLine = self.convertToDateTime(z.text)
def __str__(self):
retstr = "sarProcessingInformation:"+sep+tab
retlst = ()
retstr += "lutApplied=%s"+sep+tab
retlst += (self.lutApplied,)
retstr += "numberOfLineProcessed=%s"+sep
retlst += (self.numberOfLinesProcessed,)
retstr += "%s"+sep
retlst += (str(self.azimuthWindow),)
retstr += "%s"+sep
retlst += (str(self.rangeWindow),)
retstr += ":sarProcessingInformation"+sep
return retstr % retlst
class _Window(LT1ABNamespace):
def __init__(self,type):
LT1ABNamespace.__init__(self)
self.type = type
self.windowName = None
self.windowCoefficient = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('windowName'):
self.windowName = z.text
elif z.tag == self.elementName('windowCoefficient'):
self.windowCoefficient = float(z.text)
def __str__(self):
retstr = "%sWindow:"+sep+tab
retlst = (self.type,)
retstr += "windowName=%s"+sep+tab
retlst += (self.windowName,)
retstr += "windowCoefficient=%s"+sep
retlst += (self.windowCoefficient,)
retstr += ":%sWindow"
retlst += (self.type,)
return retstr % retlst
class _DopplerCentroid(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.timeOfDopplerCentroidEstimate = None
self.dopplerAmbiguity = None
self.dopplerAmbiguityConfidence= None
self.dopplerCentroidReferenceTime = None
self.dopplerCentroidPolynomialPeriod = None
self.dopplerCentroidCoefficients = []
self.dopplerCentroidConfidence = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('timeOfDopplerCentroidEstimate'):
self.timeOfDopplerCentroidEstimate = True
elif z.tag == self.elementName('dopplerAmbiguity'):
self.dopplerAmbiguity = True
elif z.tag == self.elementName('referencePoint'):
self.dopplerCentroidReferenceTime = float(z.text)
coefficientlist=node.findall('coefficient')
for i in range(len(coefficientlist)):
self.dopplerCentroidCoefficients.append(0)
for i in range(len(coefficientlist)):
if int(coefficientlist[i].attrib["exponent"])==i:
self.dopplerCentroidCoefficients[i]=float(coefficientlist[i].text)
def __str__(self):
retstr = "DopplerCentroid:"+sep+tab
retlst = ()
retstr += "dopplerAmbiguity=%s"+sep+tab
retlst += (self.dopplerAmbiguity,)
retstr += "dopplerCentroidCoefficients=%s"+sep
retlst += (self.dopplerCentroidCoefficients,)
retstr += ":DopplerCentroid"+sep
return retstr % retlst
class _DopplerRateValues(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.dopplerRateReferenceTime = None
self.dopplerRateValuesCoefficients = []
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('referencePoint'):
self.dopplerRateReferenceTime = float(z.text)
coefficientlist=node.findall('coefficient')
for i in range(len(coefficientlist)):
self.dopplerRateValuesCoefficients.append(0)
for i in range(len(coefficientlist)):
if int(coefficientlist[i].attrib["exponent"])==i:
self.dopplerRateValuesCoefficients[i]=float(coefficientlist[i].text)
def __str__(self):
retstr = "DopplerRateValues:"+sep+tab
retlst = ()
retstr += "dopplerRateReferenceTime=%s"+sep+tab
retlst += (self.dopplerRateReferenceTime,)
retstr += "dopplerRateValuesCoefficients=%s"+sep+tab
retlst += (self.dopplerRateValuesCoefficients,)
retstr += ":DopplerRateValues"
return retstr % retlst
class _Chirp(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
class _SlantRangeToGroundRange(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.zeroDopplerAzimuthTime = None
self.slantRangeTimeToFirstRangeSample = None
self.groundRangeOrigin = None
self.groundToSlantRangeCoefficients = []
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('zeroDopplerAzimuthTime'):
self.zeroDopplerAzimuthTime = self.convertToDateTime(z.text)
elif z.tag == self.elementName('slantRangeTimeToFirstRangeSample'):
self.slantRangeTimeToFirstRangeSample = float(z.text)
class _ImageAttributes(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.productFormat = None
self.outputMediaInterleaving = None
self.rasterAttributes = _RasterAttributes()
self.geographicInformation = _GeographicInformation()
self.fullResolutionImageData = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('productFormat'):
self.productFormat = z.text
elif z.tag == self.elementName('outputMediaInterleaving'):
self.outputMediaInterleaving = z.text
elif z.tag == self.elementName('rasterAttributes'):
self.rasterAttributes.set_from_etnode(z)
elif z.tag == self.elementName('geographicInformation'):
self.geographicInformation.set_from_etnode(z)
elif z.tag == self.elementName('fullResolutionImageData'):
self.fullResolutionImageData = z.text
class _RasterAttributes(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.dataType = None
self.bitsPerSample = []
self.numberOfSamplesPerLine = None
self.numberOfLines = None
self.sampledPixelSpacing = None
self.sampledLineSpacing = None
self.lineTimeOrdering = None
self.pixelTimeOrdering = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('dataType'):
self.dataType = z.text
elif z.tag == self.elementName('bitsPerSample'):
self.bitsPerSample.append(z.text) # TODO: Make this a dictionary with keys of 'imaginary' and 'real'
elif z.tag == self.elementName('numberOfSamplesPerLine'):
self.numberOfSamplesPerLine = int(z.text)
elif z.tag == self.elementName('numberOfLines'):
self.numberOfLines = int(z.text)
elif z.tag == self.elementName('sampledPixelSpacing'):
self.sampledPixelSpacing = float(z.text)
elif z.tag == self.elementName('sampledLineSpacing'):
self.sampledLineSpacing = float(z.text)
elif z.tag == self.elementName('lineTimeOrdering'):
self.lineTimeOrdering = z.text
elif z.tag == self.elementName('pixelTimeOrdering'):
self.pixelTimeOrdering = z.text
class _GeographicInformation(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.geolocationGrid = _GeolocationGrid()
self.rationalFunctions = _RationalFunctions()
self.referenceEllipsoidParameters = _ReferenceEllipsoidParameters()
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('geolocationGrid'):
self.geolocationGrid.set_from_etnode(z)
class _GeolocationGrid(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.imageTiePoint = []
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('imageTiePoint'):
tp = _ImageTiePoint()
tp.set_from_etnode(z)
self.imageTiePoint.append(tp)
class _ImageTiePoint(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.imageCoordinates = _ImageCoordinates()
self.geodeticCoordinates = _GeodeticCoordinates()
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('imageCoordinate'):
self.imageCoordinates.set_from_etnode(z)
elif z.tag == self.elementName('geodeticCoordinate'):
self.geodeticCoordinates.set_from_etnode(z)
def __str__(self):
retstr = "ImageTiePoint:"+sep+tab
retlst = ()
retstr += "%s"
retlst += (str(self.imageCoordinates),)
retstr += "%s"
retlst += (str(self.geodeticCoordinates),)
retstr += ":ImageTiePoint"
return retstr % retlst
class _ImageCoordinates(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.line = None
self.pixel = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('line'):
self.line = float(z.text)
elif z.tag == self.elementName('pixel'):
self.pixel = float(z.text)
def __str__(self):
retstr = "ImageCoordinate:"+sep+tab
retlst = ()
retstr += "line=%s"+sep+tab
retlst += (self.line,)
retstr += "pixel=%s"+sep+tab
retlst += (self.pixel,)
retstr += ":ImageCoordinate"
return retstr % retlst
class _GeodeticCoordinates(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.latitude = None
self.longitude = None
self.height = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('latitude'):
self.latitude = float(z.text)
elif z.tag == self.elementName('longitude'):
self.longitude = float(z.text)
elif z.tag == self.elementName('height'):
self.height = float(z.text)
def __str__(self):
retstr = "GeodeticCoordinate:"+sep+tab
retlst = ()
retstr += "latitude=%s"+sep+tab
retlst += (self.latitude,)
retstr += "longitude=%s"+sep+tab
retlst += (self.longitude,)
retstr += "height=%s"+sep+tab
retlst += (self.height,)
retstr += ":GeodeticCoordinate"
return retstr % retlst
class _ReferenceEllipsoidParameters(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
self.ellipsoidName = None
self.semiMajorAxis = None
self.semiMinorAxis = None
self.geodeticTerrainHeight = None
def set_from_etnode(self,node):
for z in node.getchildren():
if z.tag == self.elementName('ellipsoidName'):
self.ellipsoidName = z.text
elif z.tag == self.elementName('semiMajorAxis'):
self.semiMajorAxis = float(z.text)
elif z.tag == self.elementName('semiMinorAxis'):
self.semiMinorAxis = float(z.text)
elif z.tag == self.elementName('geodeticTerrainHeight'):
self.geodeticTerrainHeight = float(z.text)
def __str__(self):
return ""
class _RationalFunctions(LT1ABNamespace):
def __init__(self):
LT1ABNamespace.__init__(self)
def set_from_etnode(self,node):
pass
def __str__(self):
return ""
def findPreciseOrbit(dirname, fname, year):
'''
Find precise orbit file in given folder.
'''
import glob
###First try root folder itself
res = glob.glob( os.path.join(dirname, fname.lower()))
if len(res) == 0:
res = glob.glob( os.path.join(dirname, "{0}".format(year), fname.lower()))
if len(res) == 0:
raise Exception('Orbit Dirname provided but no suitable orbit file found in {0}'.format(dirname))
if len(res) > 1:
print('More than one matching result found. Using first result.')
return res[0]
def convertRSTimeToDateTime(instr):
'''
Convert RS2 orbit time string to datetime.
'''
parts = instr.strip().split('-')
tparts = parts[-1].split(':')
secs = float(tparts[2])
intsecs = int(secs)
musecs = int((secs - intsecs)*1e6)
timestamp = datetime.datetime(int(parts[0]),1,1, int(tparts[0]), int(tparts[1]), intsecs, musecs) + datetime.timedelta(days = int(parts[1])-1)
return timestamp