#!/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 = np.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 = np.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 = np.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('******************************************')