314 lines
12 KiB
Python
314 lines
12 KiB
Python
|
#!/usr/bin/env python3
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
import argparse
|
||
|
|
import os
|
||
|
|
import isce
|
||
|
|
import isceobj
|
||
|
|
import datetime
|
||
|
|
import sys
|
||
|
|
import s1a_isce_utils as ut
|
||
|
|
|
||
|
|
def createParser():
|
||
|
|
parser = argparse.ArgumentParser( description='Generate offset field between two Sentinel swaths')
|
||
|
|
parser.add_argument('-m', '--master', type=str, dest='master', required=True,
|
||
|
|
help='Directory with the master image')
|
||
|
|
parser.add_argument('-s', '--slave', type=str, dest='slave', required=True,
|
||
|
|
help='Directory with the slave image')
|
||
|
|
parser.add_argument('-g', '--geom_masterDir', type=str, dest='geom_masterDir', default='geom_master',
|
||
|
|
help='Directory for geometry files of the master')
|
||
|
|
parser.add_argument('-c', '--coregSLCdir', type=str, dest='coregdir', default='coreg_slaves',
|
||
|
|
help='Directory with coregistered SLC data')
|
||
|
|
parser.add_argument('-a', '--azimuth_misreg', type=str, dest='misreg_az', default='',
|
||
|
|
help='A text file that contains zimuth misregistration in subpixels')
|
||
|
|
parser.add_argument('-r', '--range_misreg', type=str, dest='misreg_rng', default='',
|
||
|
|
help='A text file that contains range misregistration in meters')
|
||
|
|
parser.add_argument('-v', '--overlap', dest='overlap', action='store_true', default=False,
|
||
|
|
help='Flatten the interferograms with offsets if needed')
|
||
|
|
parser.add_argument('-o', '--overlap_dir', type=str, dest='overlapDir', default='overlap',
|
||
|
|
help='overlap directory name')
|
||
|
|
parser.add_argument('-useGPU', '--useGPU', dest='useGPU',action='store_true', default=False,
|
||
|
|
help='Allow App to use GPU when available')
|
||
|
|
return parser
|
||
|
|
|
||
|
|
def cmdLineParse(iargs=None):
|
||
|
|
'''
|
||
|
|
Command line parser.
|
||
|
|
'''
|
||
|
|
parser = createParser()
|
||
|
|
return parser.parse_args(args=iargs)
|
||
|
|
|
||
|
|
def runGeo2rdrCPU(info, rdict, misreg_az=0.0, misreg_rg=0.0):
|
||
|
|
from zerodop.geo2rdr import createGeo2rdr
|
||
|
|
from isceobj.Planet.Planet import Planet
|
||
|
|
|
||
|
|
latImage = isceobj.createImage()
|
||
|
|
latImage.load(rdict['lat'] + '.xml')
|
||
|
|
latImage.setAccessMode('READ')
|
||
|
|
|
||
|
|
lonImage = isceobj.createImage()
|
||
|
|
lonImage.load(rdict['lon'] + '.xml')
|
||
|
|
lonImage.setAccessMode('READ')
|
||
|
|
|
||
|
|
demImage = isceobj.createImage()
|
||
|
|
demImage.load(rdict['hgt'] + '.xml')
|
||
|
|
demImage.setAccessMode('READ')
|
||
|
|
|
||
|
|
misreg_az = misreg_az*info.azimuthTimeInterval
|
||
|
|
delta = datetime.timedelta(seconds=misreg_az)
|
||
|
|
print('Additional time offset applied in geo2rdr: {0} secs'.format(misreg_az))
|
||
|
|
print('Additional range offset applied in geo2rdr: {0} m'.format(misreg_rg))
|
||
|
|
|
||
|
|
|
||
|
|
#####Run Geo2rdr
|
||
|
|
planet = Planet(pname='Earth')
|
||
|
|
grdr = createGeo2rdr()
|
||
|
|
grdr.configure()
|
||
|
|
|
||
|
|
grdr.slantRangePixelSpacing = info.rangePixelSize
|
||
|
|
grdr.prf = 1.0 / info.azimuthTimeInterval
|
||
|
|
grdr.radarWavelength = info.radarWavelength
|
||
|
|
grdr.orbit = info.orbit
|
||
|
|
grdr.width = info.numberOfSamples
|
||
|
|
grdr.length = info.numberOfLines
|
||
|
|
grdr.demLength = demImage.getLength()
|
||
|
|
grdr.demWidth = demImage.getWidth()
|
||
|
|
grdr.wireInputPort(name='planet', object=planet)
|
||
|
|
grdr.numberRangeLooks = 1
|
||
|
|
grdr.numberAzimuthLooks = 1
|
||
|
|
grdr.lookSide = -1
|
||
|
|
grdr.setSensingStart(info.sensingStart - delta)
|
||
|
|
grdr.rangeFirstSample = info.startingRange - misreg_rg
|
||
|
|
grdr.dopplerCentroidCoeffs = [0.] ###Zero doppler
|
||
|
|
|
||
|
|
grdr.rangeOffsetImageName = rdict['rangeOffName']
|
||
|
|
grdr.azimuthOffsetImageName = rdict['azOffName']
|
||
|
|
grdr.demImage = demImage
|
||
|
|
grdr.latImage = latImage
|
||
|
|
grdr.lonImage = lonImage
|
||
|
|
|
||
|
|
grdr.geo2rdr()
|
||
|
|
|
||
|
|
return
|
||
|
|
|
||
|
|
|
||
|
|
def runGeo2rdrGPU(info, rdict, misreg_az=0.0, misreg_rg=0.0):
|
||
|
|
from zerodop.GPUgeo2rdr.GPUgeo2rdr import PyGeo2rdr
|
||
|
|
from isceobj.Planet.Planet import Planet
|
||
|
|
from iscesys import DateTimeUtil as DTU
|
||
|
|
|
||
|
|
latImage = isceobj.createImage()
|
||
|
|
latImage.load(rdict['lat'] + '.xml')
|
||
|
|
latImage.setAccessMode('READ')
|
||
|
|
latImage.createImage()
|
||
|
|
|
||
|
|
lonImage = isceobj.createImage()
|
||
|
|
lonImage.load(rdict['lon'] + '.xml')
|
||
|
|
lonImage.setAccessMode('READ')
|
||
|
|
lonImage.createImage()
|
||
|
|
|
||
|
|
demImage = isceobj.createImage()
|
||
|
|
demImage.load(rdict['hgt'] + '.xml')
|
||
|
|
demImage.setAccessMode('READ')
|
||
|
|
demImage.createImage()
|
||
|
|
|
||
|
|
misreg_az = misreg_az*info.azimuthTimeInterval
|
||
|
|
delta = datetime.timedelta(seconds=misreg_az)
|
||
|
|
print('Additional time offset applied in geo2rdr: {0} secs'.format(misreg_az))
|
||
|
|
print('Additional range offset applied in geo2rdr: {0} m'.format(misreg_rg))
|
||
|
|
|
||
|
|
#####Run Geo2rdr
|
||
|
|
planet = Planet(pname='Earth')
|
||
|
|
grdr = PyGeo2rdr()
|
||
|
|
|
||
|
|
grdr.setRangePixelSpacing(info.rangePixelSize)
|
||
|
|
grdr.setPRF(1.0 / info.azimuthTimeInterval)
|
||
|
|
grdr.setRadarWavelength(info.radarWavelength)
|
||
|
|
|
||
|
|
grdr.createOrbit(0, len(info.orbit.stateVectors.list))
|
||
|
|
count = 0
|
||
|
|
for sv in info.orbit.stateVectors.list:
|
||
|
|
td = DTU.seconds_since_midnight(sv.getTime())
|
||
|
|
pos = sv.getPosition()
|
||
|
|
vel = sv.getVelocity()
|
||
|
|
|
||
|
|
grdr.setOrbitVector(count, td, pos[0], pos[1], pos[2], vel[0], vel[1], vel[2])
|
||
|
|
count += 1
|
||
|
|
|
||
|
|
grdr.setOrbitMethod(0)
|
||
|
|
grdr.setWidth(info.numberOfSamples)
|
||
|
|
grdr.setLength(info.numberOfLines)
|
||
|
|
grdr.setSensingStart(DTU.seconds_since_midnight(info.sensingStart -delta))
|
||
|
|
grdr.setRangeFirstSample(info.startingRange - misreg_rg)
|
||
|
|
grdr.setNumberRangeLooks(1)
|
||
|
|
grdr.setNumberAzimuthLooks(1)
|
||
|
|
grdr.setEllipsoidMajorSemiAxis(planet.ellipsoid.a)
|
||
|
|
grdr.setEllipsoidEccentricitySquared(planet.ellipsoid.e2)
|
||
|
|
grdr.createPoly(0, 0., 1.)
|
||
|
|
grdr.setPolyCoeff(0, 0.)
|
||
|
|
|
||
|
|
grdr.setDemLength(demImage.getLength())
|
||
|
|
grdr.setDemWidth(demImage.getWidth())
|
||
|
|
grdr.setBistaticFlag(0)
|
||
|
|
|
||
|
|
rangeOffsetImage = isceobj.createImage()
|
||
|
|
rangeOffsetImage.setFilename(rdict['rangeOffName'])
|
||
|
|
rangeOffsetImage.setAccessMode('write')
|
||
|
|
rangeOffsetImage.setDataType('FLOAT')
|
||
|
|
rangeOffsetImage.setCaster('write', 'DOUBLE')
|
||
|
|
rangeOffsetImage.setWidth(demImage.width)
|
||
|
|
rangeOffsetImage.createImage()
|
||
|
|
|
||
|
|
azimuthOffsetImage = isceobj.createImage()
|
||
|
|
azimuthOffsetImage.setFilename(rdict['azOffName'])
|
||
|
|
azimuthOffsetImage.setAccessMode('write')
|
||
|
|
azimuthOffsetImage.setDataType('FLOAT')
|
||
|
|
azimuthOffsetImage.setCaster('write', 'DOUBLE')
|
||
|
|
azimuthOffsetImage.setWidth(demImage.width)
|
||
|
|
azimuthOffsetImage.createImage()
|
||
|
|
|
||
|
|
grdr.setLatAccessor(latImage.getImagePointer())
|
||
|
|
grdr.setLonAccessor(lonImage.getImagePointer())
|
||
|
|
grdr.setHgtAccessor(demImage.getImagePointer())
|
||
|
|
grdr.setAzAccessor(0)
|
||
|
|
grdr.setRgAccessor(0)
|
||
|
|
grdr.setAzOffAccessor(azimuthOffsetImage.getImagePointer())
|
||
|
|
grdr.setRgOffAccessor(rangeOffsetImage.getImagePointer())
|
||
|
|
|
||
|
|
grdr.geo2rdr()
|
||
|
|
|
||
|
|
rangeOffsetImage.finalizeImage()
|
||
|
|
rangeOffsetImage.renderHdr()
|
||
|
|
|
||
|
|
azimuthOffsetImage.finalizeImage()
|
||
|
|
azimuthOffsetImage.renderHdr()
|
||
|
|
latImage.finalizeImage()
|
||
|
|
lonImage.finalizeImage()
|
||
|
|
demImage.finalizeImage()
|
||
|
|
|
||
|
|
return
|
||
|
|
pass
|
||
|
|
|
||
|
|
def main(iargs=None):
|
||
|
|
'''
|
||
|
|
Estimate offsets for the overlap regions of the bursts.
|
||
|
|
'''
|
||
|
|
inps = cmdLineParse(iargs)
|
||
|
|
|
||
|
|
# see if the user compiled isce with GPU enabled
|
||
|
|
run_GPU = False
|
||
|
|
try:
|
||
|
|
from zerodop.GPUtopozero.GPUtopozero import PyTopozero
|
||
|
|
from zerodop.GPUgeo2rdr.GPUgeo2rdr import PyGeo2rdr
|
||
|
|
run_GPU = True
|
||
|
|
except:
|
||
|
|
pass
|
||
|
|
|
||
|
|
if inps.useGPU and not run_GPU:
|
||
|
|
print("GPU mode requested but no GPU ISCE code found")
|
||
|
|
|
||
|
|
# setting the respective version of geo2rdr for CPU and GPU
|
||
|
|
if run_GPU and inps.useGPU:
|
||
|
|
print('GPU mode')
|
||
|
|
runGeo2rdr = runGeo2rdrGPU
|
||
|
|
else:
|
||
|
|
print('CPU mode')
|
||
|
|
runGeo2rdr = runGeo2rdrCPU
|
||
|
|
|
||
|
|
|
||
|
|
masterSwathList = ut.getSwathList(inps.master)
|
||
|
|
slaveSwathList = ut.getSwathList(inps.slave)
|
||
|
|
|
||
|
|
swathList = list(sorted(set(masterSwathList+slaveSwathList)))
|
||
|
|
|
||
|
|
for swath in swathList:
|
||
|
|
##Load slave metadata
|
||
|
|
slave = ut.loadProduct(os.path.join(inps.slave, 'IW{0}.xml'.format(swath)))
|
||
|
|
master = ut.loadProduct(os.path.join(inps.master, 'IW{0}.xml'.format(swath)))
|
||
|
|
|
||
|
|
### output directory
|
||
|
|
if inps.overlap:
|
||
|
|
outdir = os.path.join(inps.coregdir, inps.overlapDir, 'IW{0}'.format(swath))
|
||
|
|
else:
|
||
|
|
outdir = os.path.join(inps.coregdir, 'IW{0}'.format(swath))
|
||
|
|
|
||
|
|
if not os.path.isdir(outdir):
|
||
|
|
os.makedirs(outdir)
|
||
|
|
|
||
|
|
if os.path.exists(str(inps.misreg_az)):
|
||
|
|
with open(inps.misreg_az, 'r') as f:
|
||
|
|
misreg_az = float(f.readline())
|
||
|
|
else:
|
||
|
|
misreg_az = 0.0
|
||
|
|
|
||
|
|
if os.path.exists(str(inps.misreg_rng)):
|
||
|
|
with open(inps.misreg_rng, 'r') as f:
|
||
|
|
misreg_rg = float(f.readline())
|
||
|
|
else:
|
||
|
|
misreg_rg = 0.0
|
||
|
|
|
||
|
|
burstoffset, minBurst, maxBurst = master.getCommonBurstLimits(slave)
|
||
|
|
|
||
|
|
###Burst indices w.r.t master
|
||
|
|
if inps.overlap:
|
||
|
|
maxBurst = maxBurst - 1
|
||
|
|
geomDir = os.path.join(inps.geom_masterDir, inps.overlapDir, 'IW{0}'.format(swath))
|
||
|
|
|
||
|
|
else:
|
||
|
|
geomDir = os.path.join(inps.geom_masterDir, 'IW{0}'.format(swath))
|
||
|
|
|
||
|
|
|
||
|
|
slaveBurstStart = minBurst + burstoffset
|
||
|
|
|
||
|
|
for mBurst in range(minBurst, maxBurst):
|
||
|
|
|
||
|
|
###Corresponding slave burst
|
||
|
|
sBurst = slaveBurstStart + (mBurst - minBurst)
|
||
|
|
burstTop = slave.bursts[sBurst]
|
||
|
|
if inps.overlap:
|
||
|
|
burstBot = slave.bursts[sBurst+1]
|
||
|
|
|
||
|
|
print('Overlap pair {0}: Burst {1} of master matched with Burst {2} of slave'.format(mBurst-minBurst, mBurst, sBurst))
|
||
|
|
if inps.overlap:
|
||
|
|
####Generate offsets for top burst
|
||
|
|
rdict = {'lat': os.path.join(geomDir,'lat_%02d_%02d.rdr'%(mBurst+1,mBurst+2)),
|
||
|
|
'lon': os.path.join(geomDir,'lon_%02d_%02d.rdr'%(mBurst+1,mBurst+2)),
|
||
|
|
'hgt': os.path.join(geomDir,'hgt_%02d_%02d.rdr'%(mBurst+1,mBurst+2)),
|
||
|
|
'rangeOffName': os.path.join(outdir, 'range_top_%02d_%02d.off'%(mBurst+1,mBurst+2)),
|
||
|
|
'azOffName': os.path.join(outdir, 'azimuth_top_%02d_%02d.off'%(mBurst+1,mBurst+2))}
|
||
|
|
|
||
|
|
runGeo2rdr(burstTop, rdict, misreg_az=misreg_az, misreg_rg=misreg_rg)
|
||
|
|
|
||
|
|
print('Overlap pair {0}: Burst {1} of master matched with Burst {2} of slave'.format(mBurst-minBurst, mBurst+1, sBurst+1))
|
||
|
|
####Generate offsets for bottom burst
|
||
|
|
rdict = {'lat': os.path.join(geomDir,'lat_%02d_%02d.rdr'%(mBurst+1,mBurst+2)),
|
||
|
|
'lon': os.path.join(geomDir, 'lon_%02d_%02d.rdr'%(mBurst+1,mBurst+2)),
|
||
|
|
'hgt': os.path.join(geomDir, 'hgt_%02d_%02d.rdr'%(mBurst+1,mBurst+2)),
|
||
|
|
'rangeOffName': os.path.join(outdir, 'range_bot_%02d_%02d.off'%(mBurst+1,mBurst+2)),
|
||
|
|
'azOffName': os.path.join(outdir, 'azimuth_bot_%02d_%02d.off'%(mBurst+1,mBurst+2))}
|
||
|
|
|
||
|
|
runGeo2rdr(burstBot, rdict, misreg_az=misreg_az, misreg_rg=misreg_rg)
|
||
|
|
|
||
|
|
else:
|
||
|
|
print('Burst {1} of master matched with Burst {2} of slave'.format(mBurst-minBurst, mBurst, sBurst))
|
||
|
|
####Generate offsets for top burst
|
||
|
|
rdict = {'lat': os.path.join(geomDir,'lat_%02d.rdr'%(mBurst+1)),
|
||
|
|
'lon': os.path.join(geomDir,'lon_%02d.rdr'%(mBurst+1)),
|
||
|
|
'hgt': os.path.join(geomDir,'hgt_%02d.rdr'%(mBurst+1)),
|
||
|
|
'rangeOffName': os.path.join(outdir, 'range_%02d.off'%(mBurst+1)),
|
||
|
|
'azOffName': os.path.join(outdir, 'azimuth_%02d.off'%(mBurst+1))}
|
||
|
|
|
||
|
|
runGeo2rdr(burstTop, rdict, misreg_az=misreg_az, misreg_rg=misreg_rg)
|
||
|
|
|
||
|
|
|
||
|
|
|
||
|
|
if __name__ == '__main__':
|
||
|
|
'''
|
||
|
|
Generate burst-by-burst reverse geometry mapping for resampling.
|
||
|
|
'''
|
||
|
|
# Main Driver
|
||
|
|
main()
|
||
|
|
|
||
|
|
|
||
|
|
|