267 lines
9.7 KiB
Python
267 lines
9.7 KiB
Python
#!/usr/bin/env python3
|
|
from __future__ import print_function
|
|
import argparse
|
|
import isce
|
|
from make_raw import makeRawApp
|
|
import numpy as np
|
|
import os
|
|
import itertools
|
|
from isceobj.XmlUtil.XmlUtil import XmlUtil
|
|
from isceobj.Orbit.Orbit import Orbit, StateVector
|
|
from iscesys.StdOEL.StdOELPy import create_writer
|
|
#import sarxml
|
|
import stdproc
|
|
import datetime
|
|
|
|
|
|
|
|
stdWriter = create_writer("log", "", True, filename="prepareStack.log")
|
|
|
|
def pulseTiming(frame):
|
|
#From runPulseTiming() in InsarProc
|
|
numberOfLines = frame.getNumberOfLines()
|
|
prf = frame.getInstrument().getPulseRepetitionFrequency()
|
|
pri = 1.0 / prf
|
|
startTime = frame.getSensingStart()
|
|
orbit = frame.getOrbit()
|
|
|
|
pulseOrbit = Orbit()
|
|
startTimeUTC0 = (startTime - datetime.datetime(startTime.year,startTime.month,startTime.day))
|
|
timeVec = [pri*i + startTimeUTC0.seconds + 10**-6*startTimeUTC0.microseconds for i in xrange(numberOfLines)]
|
|
for i in range(numberOfLines):
|
|
dt = i * pri
|
|
time = startTime + datetime.timedelta(seconds=dt)
|
|
sv = orbit.interpolateOrbit(time, method='hermite')
|
|
pulseOrbit.addStateVector(sv)
|
|
|
|
return pulseOrbit
|
|
|
|
def getPeg(planet, orbit):
|
|
#Returns relevant peg point. From runSetMocompPath.py
|
|
|
|
objPeg = stdproc.createGetpeg()
|
|
objPeg.wireInputPort(name='planet', object=planet)
|
|
objPeg.wireInputPort(name='Orbit', object=orbit)
|
|
|
|
stdWriter.setFileTag("getpeg", "log")
|
|
stdWriter.setFileTag("getpeg", "err")
|
|
stdWriter.setFileTag("getpeg", "out")
|
|
# objSetmocomppath.setStdWriter(self._stdWriter)
|
|
objPeg.setStdWriter(stdWriter)
|
|
objPeg.estimatePeg()
|
|
|
|
return objPeg.getPeg(), objPeg.getAverageHeight()
|
|
|
|
class orbit_info:
|
|
def __init__(self, sar, fname):
|
|
'''Initialize with a sarProc object and corresponding XML file name'''
|
|
orbit = pulseTiming(sar.make_raw.frame)
|
|
tim, pos, vel, offset = orbit._unpackOrbit()
|
|
planet = sar.make_raw.planet
|
|
self.tim = tim
|
|
self.pos = pos
|
|
self.vel = vel
|
|
self.dt = sar.make_raw.frame.sensingMid
|
|
self.prf = sar.make_raw.doppler.prf
|
|
self.fd = sar.make_raw.dopplerValues() * self.prf
|
|
self.nvec = len(self.tim)
|
|
self.peg, self.hgt = getPeg(planet, orbit)
|
|
self.rds = self.peg.getRadiusOfCurvature()
|
|
self.rng = sar.make_raw.frame.startingRange
|
|
self.clook = None
|
|
self.slook = None
|
|
self.filename = fname
|
|
self.computeLookAngle()
|
|
|
|
def computeLookAngle(self):
|
|
self.clook = (2*self.hgt*self.rds+self.hgt**2+self.rng**2)/(2*self.rng*(self.rds+self.hgt))
|
|
self.slook = np.sqrt(1-self.clook**2)
|
|
# print('Estimated Look Angle: %3.2f degrees'%(np.arccos(self.clook)*180.0/np.pi))
|
|
|
|
def getBaseline(self, secondary):
|
|
'''Compute baseline between current object and another orbit object.'''
|
|
|
|
ind = int(self.nvec/2)
|
|
|
|
mpos = np.array(self.pos[ind])
|
|
mvel = np.array(self.vel[ind])
|
|
|
|
#######From the ROI-PAC scripts
|
|
rvec = mpos/np.linalg.norm(mpos)
|
|
crp = np.cross(rvec, mvel)/np.linalg.norm(mvel)
|
|
crp = crp/np.linalg.norm(crp)
|
|
vvec = np.cross(crp, rvec)
|
|
mvel = np.linalg.norm(mvel)
|
|
|
|
ind = int(secondary.nvec/2) #First guess
|
|
spos = np.array(secondary.pos[ind])
|
|
svel = np.array(secondary.vel[ind])
|
|
svel = np.linalg.norm(svel)
|
|
|
|
dx = spos - mpos;
|
|
z_offset = secondary.prf*np.dot(dx, vvec)/mvel
|
|
|
|
ind = int(ind - z_offset) #Refined estimate
|
|
spos = secondary.pos[ind]
|
|
svel = secondary.vel[ind]
|
|
svel = np.linalg.norm(svel)
|
|
|
|
dx = spos-mpos
|
|
hb = np.dot(dx, crp)
|
|
vb = np.dot(dx, rvec)
|
|
|
|
csb = -1.0*hb*self.clook + vb*self.slook
|
|
|
|
# print('Estimated Baseline: %4.2f'%csb)
|
|
return csb
|
|
|
|
|
|
def parse():
|
|
|
|
# class RangeObj(object):
|
|
# '''Class to deal with input ranges.'''
|
|
# def __init__(self, start, end):
|
|
# self.start = start
|
|
# self.end = end
|
|
# def __eq__(self, other):
|
|
# return self.start <= other <= self.end
|
|
|
|
|
|
def Range(nmin, nmax):
|
|
class RangeObj(argparse.Action):
|
|
def __call__(self, parser, args, values, option_string=None):
|
|
if not nmin <= values <= nmax:
|
|
msg = 'Argument "{f}" requires value between {nmin} and {nmax}'.format(f=self.dest, nmin=nmin, nmax=nmax)
|
|
raise argparse.ArgumentTypeError(msg)
|
|
setattr(args, self.dest, values)
|
|
|
|
return RangeObj
|
|
|
|
#####Actual parser set up
|
|
parser = argparse.ArgumentParser(description='Computes the baseline plot for given set of SAR images.')
|
|
parser.add_argument('fnames', nargs='+', default=None, help = 'XML files corresponding to the SAR scenes.')
|
|
parser.add_argument('-Bcrit', dest='Bcrit', default=1200.0, help='Critical Geometric Baseline in meters [0., 10000.]', type=float, action=Range(0., 10000.))
|
|
parser.add_argument('-Tau', dest='Tau', default=1080.0, help='Temporal Decorrelation Time Constant in days [0., 3650.]', type=float, action=Range(0., 3650.))
|
|
parser.add_argument('-dop', dest='dop', default=0.5, help='Critical Doppler difference in fraction of PRF', type=float, action=Range(0., 1.))
|
|
parser.add_argument('-coh', dest='cThresh', default=0.3, help='Coherence Threshold to estimate viable interferograms. [0., 1.0]', type=float, action=Range(0., 1.))
|
|
parser.add_argument('-dir', dest='dirname', default='insar_XML', help='Directory in which the individual insar XML files are created.', type=str, action='store')
|
|
parser.add_argument('-base', dest='base', default='base.xml', help='Base XML for the insar.xml files.', type=str)
|
|
inps = parser.parse_args()
|
|
|
|
return inps
|
|
|
|
if __name__ == '__main__':
|
|
inps = parse()
|
|
nSar = len(inps.fnames)
|
|
print(inps.fnames)
|
|
print('Number of SAR Scenes = %d'%nSar)
|
|
|
|
Orbits = []
|
|
print('Reading in all the raw files and metadata.')
|
|
for k in xrange(nSar):
|
|
sar = makeRawApp()
|
|
sar.run(inps.fnames[k])
|
|
Orbits.append(orbit_info(sar, inps.fnames[k]))
|
|
|
|
##########We now have all the pegpoints to start processing.
|
|
Dopplers = np.zeros(nSar)
|
|
Bperp = np.zeros(nSar)
|
|
Days = np.zeros(nSar)
|
|
|
|
#######Setting the first scene as temporary reference.
|
|
reference = Orbits[0]
|
|
|
|
|
|
Dopplers[0] = reference.fd
|
|
Days[0] = reference.dt.toordinal()
|
|
for k in xrange(1,nSar):
|
|
secondary = Orbits[k]
|
|
Bperp[k] = reference.getBaseline(secondary)
|
|
Dopplers[k] = secondary.fd
|
|
Days[k] = secondary.dt.toordinal()
|
|
|
|
|
|
print("************************************")
|
|
print("Index Date Bperp Doppler")
|
|
print("************************************")
|
|
|
|
for k in xrange(nSar):
|
|
print('{0:>3} {1:>10} {2:4.2f} {3:4.2f}'.format(k+1, Orbits[k].dt.strftime('%Y-%m-%d'), Bperp[k],Dopplers[k]))
|
|
|
|
|
|
print("************************************")
|
|
|
|
geomRho = (1-np.clip(np.abs(Bperp[:,None]-Bperp[None,:])/inps.Bcrit, 0., 1.))
|
|
tempRho = np.exp(-1.0*np.abs(Days[:,None]-Days[None,:])/inps.Tau)
|
|
dopRho = (np.abs(Dopplers[:,None] - Dopplers[None,:])/ reference.prf) < inps.dop
|
|
|
|
Rho = geomRho * tempRho * dopRho
|
|
for kk in xrange(nSar):
|
|
Rho[kk,kk] = 0.
|
|
|
|
|
|
avgRho = np.mean(Rho, axis=1)*nSar/(nSar-1)
|
|
numViable = np.sum((Rho> inps.cThresh), axis=1)
|
|
|
|
####Currently sorting on average coherence.
|
|
|
|
referenceChoice = np.argsort(avgRho)
|
|
referenceOrbit = Orbits[referenceChoice[0]]
|
|
referenceBperp = Bperp[referenceChoice[0]]
|
|
|
|
|
|
print('*************************************')
|
|
print('Ranking for Reference Scene Selection: ')
|
|
print('**************************************')
|
|
print('Rank Index Date nViable Avg. Coh.' )
|
|
for kk in xrange(nSar):
|
|
ind = referenceChoice[kk]
|
|
print('{0:>3} {1:>3} {2:>10} {3:>4} {4:>2.3f}'.format(kk+1, ind+1, Orbits[ind].dt.strftime('%Y-%m-%d'), numViable[ind], avgRho[ind]))
|
|
|
|
print('***************************************')
|
|
|
|
print('***************************************')
|
|
print('List of Viable interferograms:')
|
|
print('***************************************')
|
|
|
|
# if not os.path.isdir(inps.dirname):
|
|
# try:
|
|
# os.mkdir(inps.dirname)
|
|
# except:
|
|
# raise OSError("%s Directory cannot be created"%(inps.dirname))
|
|
|
|
|
|
|
|
[ii,jj] = np.where(Rho > inps.cThresh)
|
|
|
|
print('Reference Secondary Bperp Deltat')
|
|
for mind, sind in itertools.izip(ii,jj):
|
|
reference = Orbits[mind]
|
|
secondary = Orbits[sind]
|
|
if reference.dt > secondary.dt:
|
|
print('{0:>10} {1:>10} {2:>4.2f} {3:>4.2f}'.format(reference.dt.strftime('%Y-%m-%d'), secondary.dt.strftime('%Y-%m-%d'), Bperp[mind]-Bperp[sind], Days[mind] - Days[sind]))
|
|
xmlname = '%s/insar_%s_%s.xml'%(inps.dirname, reference.dt.strftime('%Y%m%d'), secondary.dt.strftime('%Y%m%d'))
|
|
|
|
# sarxml.sartoinsarXML(reference.filename, secondary.filename, base=inps.base, out=xmlname)
|
|
|
|
|
|
print('***************************************')
|
|
|
|
#######Currently picks reference peg point.
|
|
print('***************************************')
|
|
commonPeg = referenceOrbit.peg
|
|
print('Common peg point: ')
|
|
print(commonPeg)
|
|
print('Bperp Range: [%f , %f] '%(Bperp.min()-referenceBperp, Bperp.max()-referenceBperp))
|
|
|
|
######Choose median doppler
|
|
commonDop = np.median(Dopplers)
|
|
maxDop = np.max(Dopplers)
|
|
minDop = np.min(Dopplers)
|
|
varDop = np.max(np.abs(Dopplers-commonDop))/referenceOrbit.prf
|
|
|
|
print('Common Doppler: ', commonDop)
|
|
print('Doppler Range: [%f, %f]'%(minDop, maxDop))
|
|
print('MAx Doppler Variation = %f %%'%(varDop*100))
|
|
print('******************************************')
|