409 lines
15 KiB
Python
409 lines
15 KiB
Python
#!/usr/bin/env python3
|
|
|
|
import numpy as np
|
|
import argparse
|
|
import os
|
|
import isce
|
|
import isceobj
|
|
import copy
|
|
import logging
|
|
import scipy.signal as SS
|
|
from isceobj.Util.ImageUtil import ImageLib as IML
|
|
import s1a_isce_utils as ut
|
|
|
|
|
|
def createParser():
|
|
parser = argparse.ArgumentParser( description='Generate offset field between two Sentinel swaths')
|
|
|
|
parser.add_argument('-i', '--interferogram', type=str, dest='interferogram',required=True,
|
|
help='Directory with the overlap interferogram')
|
|
parser.add_argument('-m', '--reference_dir', type=str, dest='reference', required=True,
|
|
help='Directory with the secondary image')
|
|
parser.add_argument('-s', '--secondary_dir', type=str, dest='secondary', required=True,
|
|
help='Directory with the secondary image')
|
|
parser.add_argument('-d', '--overlap_dir', type=str, dest='overlap', required=True,
|
|
help='Directory with overlap products')
|
|
|
|
parser.add_argument('-a', '--esdAzimuthLooks', type=int, dest='esdAzimuthLooks', default = 5,
|
|
help='ESD azimuth looks')
|
|
|
|
parser.add_argument('-r', '--esdRangeLooks', type=int, dest='esdRangeLooks', default = 15,
|
|
help='ESD range looks')
|
|
return parser
|
|
|
|
|
|
|
|
def cmdLineParse(iargs=None):
|
|
'''
|
|
Command line parser.
|
|
'''
|
|
parser = createParser()
|
|
return parser.parse_args(args=iargs)
|
|
|
|
|
|
|
|
def multilook(intName, alks=5, rlks=15):
|
|
cmd = 'looks.py -i {0} -a {1} -r {2}'.format(intName,alks,rlks)
|
|
flag = os.system(cmd)
|
|
|
|
if flag:
|
|
raise Exception('Failed to multilook %s'%(intName))
|
|
|
|
spl = os.path.splitext(intName)
|
|
return '{0}.{1}alks_{2}rlks{3}'.format(spl[0],alks,rlks,spl[1])
|
|
|
|
|
|
|
|
def overlapSpectralSeparation(topBurstIfg, botBurstIfg, referenceTop, referenceBot, secondaryTop, secondaryBot, azMasTop, rgMasTop, azMasBot, rgMasBot, azSlvTop, rgSlvTop, azSlvBot, rgSlvBot , misreg=0.0):
|
|
# Added by Heresh Fattahi
|
|
'''
|
|
Estimate separation in frequency due to unit pixel misregistration.
|
|
'''
|
|
'''
|
|
dt = topBurstIfg.azimuthTimeInterval
|
|
topStart = int(np.round((topBurstIfg.sensingStart - referenceTop.sensingStart).total_seconds() / dt))
|
|
overlapLen = topBurstIfg.numberOfLines
|
|
botStart = int(np.round((botBurstIfg.sensingStart - referenceBot.sensingStart).total_seconds() / dt))
|
|
|
|
print(topBurstIfg.sensingStart, referenceTop.sensingStart)
|
|
print(botBurstIfg.sensingStart, referenceBot.sensingStart)
|
|
print(topStart, botStart, overlapLen)
|
|
'''
|
|
print ('++++++++++++++++++++++')
|
|
dt = topBurstIfg.azimuthTimeInterval
|
|
topStart = int ( np.round( (referenceBot.sensingStart - referenceTop.sensingStart).total_seconds()/dt))+ referenceBot.firstValidLine
|
|
overlapLen = topBurstIfg.numberOfLines
|
|
botStart = referenceBot.firstValidLine
|
|
print(topStart, botStart, overlapLen)
|
|
#print(Debug)
|
|
|
|
##############
|
|
# reference top : m1
|
|
|
|
|
|
|
|
y = np.arange(topStart, topStart+overlapLen)[:,None] * np.ones((overlapLen, topBurstIfg.numberOfSamples))
|
|
x = np.ones((overlapLen, topBurstIfg.numberOfSamples)) * np.arange(topBurstIfg.numberOfSamples)[None,:]
|
|
|
|
if os.path.exists(azMasTop) and os.path.exists(rgMasTop):
|
|
yy = np.memmap( azMasTop, dtype=np.float32, mode='r',
|
|
shape=(topBurstIfg.numberOfLines, topBurstIfg.numberOfSamples))
|
|
xx = np.memmap( rgMasTop, dtype=np.float32, mode='r',
|
|
shape=(topBurstIfg.numberOfLines, topBurstIfg.numberOfSamples))
|
|
else:
|
|
yy = 0.0
|
|
xx = 0.0
|
|
|
|
|
|
azi = y + yy
|
|
rng = x + xx
|
|
|
|
Vs = np.linalg.norm(referenceTop.orbit.interpolateOrbit(referenceTop.sensingMid, method='hermite').getVelocity())
|
|
Ks = 2 * Vs * referenceTop.azimuthSteeringRate / referenceTop.radarWavelength
|
|
rng = referenceTop.startingRange + referenceTop.rangePixelSize * rng
|
|
Ka = referenceTop.azimuthFMRate(rng)
|
|
|
|
Ktm1 = Ks / (1.0 - Ks / Ka)
|
|
tm1 = (azi - (referenceTop.numberOfLines//2)) * referenceTop.azimuthTimeInterval
|
|
|
|
fm1 = referenceTop.doppler(rng)
|
|
|
|
##############
|
|
# reference bottom : m2
|
|
y = np.arange(botStart, botStart + overlapLen)[:,None] * np.ones((overlapLen, botBurstIfg.numberOfSamples))
|
|
x = np.ones((overlapLen, botBurstIfg.numberOfSamples)) * np.arange(botBurstIfg.numberOfSamples)[None,:]
|
|
|
|
if os.path.exists(azMasBot) and os.path.exists(rgMasBot):
|
|
yy = np.memmap( azMasBot, dtype=np.float32, mode='r',
|
|
shape=(botBurstIfg.numberOfLines, botBurstIfg.numberOfSamples))
|
|
xx = np.memmap( rgMasBot, dtype=np.float32, mode='r',
|
|
shape=(botBurstIfg.numberOfLines, botBurstIfg.numberOfSamples))
|
|
else:
|
|
yy = 0.0
|
|
xx = 0.0
|
|
|
|
azi = y + yy
|
|
rng = x + xx
|
|
|
|
Vs = np.linalg.norm(referenceBot.orbit.interpolateOrbit(referenceBot.sensingMid, method='hermite').getVelocity())
|
|
Ks = 2 * Vs * referenceBot.azimuthSteeringRate / referenceBot.radarWavelength
|
|
rng = referenceBot.startingRange + referenceBot.rangePixelSize * rng
|
|
Ka = referenceBot.azimuthFMRate(rng)
|
|
|
|
Ktm2 = Ks / (1.0 - Ks / Ka)
|
|
tm2 = (azi - (referenceBot.numberOfLines//2)) * referenceBot.azimuthTimeInterval
|
|
fm2 = referenceBot.doppler(rng)
|
|
|
|
|
|
##############
|
|
# secondary top : s1
|
|
y = np.arange(topStart, topStart+overlapLen)[:,None] * np.ones((overlapLen, topBurstIfg.numberOfSamples))
|
|
x = np.ones((overlapLen, topBurstIfg.numberOfSamples)) * np.arange(topBurstIfg.numberOfSamples)[None,:]
|
|
|
|
if os.path.exists(azSlvTop) and os.path.exists(rgSlvTop):
|
|
yy = np.memmap( azSlvTop, dtype=np.float32, mode='r',
|
|
shape=(topBurstIfg.numberOfLines, topBurstIfg.numberOfSamples))
|
|
xx = np.memmap( rgSlvTop, dtype=np.float32, mode='r',
|
|
shape=(topBurstIfg.numberOfLines, topBurstIfg.numberOfSamples))
|
|
else:
|
|
yy = 0.0
|
|
xx = 0.0
|
|
|
|
|
|
azi = y + yy + misreg
|
|
rng = x + xx
|
|
|
|
# print('Azi top: ', azi[0,0], azi[-1,-1])
|
|
# print('YY top: ', yy[0,0], yy[-1,-1])
|
|
# print('Rng top: ', rng[0,0], azi[-1,-1])
|
|
# print('XX top: ', xx[0,0], xx[-1,-1])
|
|
|
|
Vs = np.linalg.norm(secondaryTop.orbit.interpolateOrbit(secondaryTop.sensingMid, method='hermite').getVelocity())
|
|
Ks = 2 * Vs * secondaryTop.azimuthSteeringRate / secondaryTop.radarWavelength
|
|
rng = secondaryTop.startingRange + secondaryTop.rangePixelSize * rng
|
|
Ka = secondaryTop.azimuthFMRate(rng)
|
|
|
|
Kts1 = Ks / (1.0 - Ks / Ka)
|
|
ts1 = (azi - (secondaryTop.numberOfLines//2)) * secondaryTop.azimuthTimeInterval
|
|
fs1 = secondaryTop.doppler(rng)
|
|
|
|
|
|
|
|
##############
|
|
# secondary bot : s2
|
|
y = np.arange(botStart, botStart + overlapLen)[:,None] * np.ones((overlapLen, botBurstIfg.numberOfSamples))
|
|
x = np.ones((overlapLen, botBurstIfg.numberOfSamples)) * np.arange(botBurstIfg.numberOfSamples)[None,:]
|
|
|
|
####Bottom secondary
|
|
if os.path.exists(azSlvBot) and os.path.exists(rgSlvBot):
|
|
yy = np.memmap( azSlvBot, dtype=np.float32, mode='r',
|
|
shape=(botBurstIfg.numberOfLines, botBurstIfg.numberOfSamples))
|
|
xx = np.memmap( rgSlvBot, dtype=np.float32, mode='r',
|
|
shape=(botBurstIfg.numberOfLines, botBurstIfg.numberOfSamples))
|
|
else:
|
|
yy = 0.0
|
|
xx = 0.0
|
|
|
|
azi = y + yy + misreg
|
|
rng = x + xx
|
|
|
|
# print('Azi bot: ', azi[0,0], azi[-1,-1])
|
|
# print('YY bot: ', yy[0,0], yy[-1,-1])
|
|
# print('Rng bot: ', rng[0,0], azi[-1,-1])
|
|
# print('XX bot: ', xx[0,0], xx[-1,-1])
|
|
|
|
Vs = np.linalg.norm(secondaryBot.orbit.interpolateOrbit(secondaryBot.sensingMid, method='hermite').getVelocity())
|
|
Ks = 2 * Vs * secondaryBot.azimuthSteeringRate / secondaryBot.radarWavelength
|
|
rng = secondaryBot.startingRange + secondaryBot.rangePixelSize * rng
|
|
Ka = secondaryBot.azimuthFMRate(rng)
|
|
Kts2 = Ks / (1.0 - Ks / Ka)
|
|
|
|
ts2 = (azi - (secondaryBot.numberOfLines//2)) * secondaryBot.azimuthTimeInterval
|
|
fs2 = secondaryBot.doppler(rng)
|
|
|
|
##############
|
|
frequencySeparation = -Ktm2*tm2 + Ktm1*tm1 + Kts1*ts1 - Kts2*ts2 + fm2 - fm1 + fs1 -fs2
|
|
#print(frequencySeparation)
|
|
#print(tm2)
|
|
#print(tm1)
|
|
#print('*********')
|
|
#print(ts1)
|
|
#print(ts2)
|
|
#print(Debug)
|
|
return frequencySeparation
|
|
|
|
|
|
def createCoherence(intfile, win=5):
|
|
'''
|
|
Compute coherence using scipy convolve 2D.
|
|
'''
|
|
|
|
corfile = os.path.splitext(intfile)[0] + '.cor'
|
|
filt = np.ones((win,win))/ (1.0*win*win)
|
|
|
|
inimg = IML.mmapFromISCE(intfile + '.xml', logging)
|
|
cJ = np.complex64(1.0j)
|
|
angle = np.exp(cJ * np.angle(inimg.bands[0]))
|
|
|
|
res = SS.convolve2d(angle, filt, mode='same')
|
|
res[0:win-1,:] = 0.0
|
|
res[-win+1:,:] = 0.0
|
|
res[:,0:win-1] = 0.0
|
|
res[:,-win+1:] = 0.0
|
|
|
|
res = np.abs(res)
|
|
|
|
with open(corfile, 'wb') as f:
|
|
res.astype(np.float32).tofile(f)
|
|
|
|
img = isceobj.createImage()
|
|
img.setFilename(corfile)
|
|
img.setWidth(res.shape[1])
|
|
img.dataType='FLOAT'
|
|
img.setAccessMode('READ')
|
|
img.renderHdr()
|
|
img.renderVRT()
|
|
# img.createImage()
|
|
# img.finalizeImage()
|
|
|
|
return corfile
|
|
|
|
def main(iargs=None):
|
|
'''
|
|
Create additional layers for performing ESD.
|
|
'''
|
|
|
|
inps = cmdLineParse(iargs)
|
|
inps.interferogram = os.path.join(inps.interferogram,'overlap')
|
|
inps.reference = os.path.join(inps.reference,'overlap')
|
|
inps.secondary = os.path.join(inps.secondary,'overlap')
|
|
|
|
referenceSwathList = ut.getSwathList(inps.reference)
|
|
secondarySwathList = ut.getSwathList(inps.secondary)
|
|
|
|
swathList = list(sorted(set(referenceSwathList+secondarySwathList)))
|
|
|
|
for swath in swathList:
|
|
IWstr = 'IW{0}'.format(swath)
|
|
referenceTop = ut.loadProduct(os.path.join(inps.reference, IWstr + '_top.xml'))
|
|
referenceBot = ut.loadProduct(os.path.join(inps.reference , IWstr + '_bottom.xml'))
|
|
|
|
secondaryTop = ut.loadProduct(os.path.join(inps.secondary, IWstr + '_top.xml'))
|
|
secondaryBot = ut.loadProduct(os.path.join(inps.secondary, IWstr + '_bottom.xml'))
|
|
|
|
|
|
####Load metadata for burst IFGs
|
|
ifgTop = ut.loadProduct(os.path.join(inps.interferogram , IWstr + '_top.xml'))
|
|
ifgBottom = ut.loadProduct(os.path.join(inps.interferogram, IWstr + '_bottom.xml'))
|
|
|
|
####Create ESD output directory
|
|
esddir = os.path.join(inps.overlap, IWstr)
|
|
os.makedirs(esddir, exist_ok=True)
|
|
|
|
####Overlap offsets directory
|
|
referenceOffdir = os.path.join(inps.reference, IWstr)
|
|
secondaryOffdir = os.path.join(inps.secondary,IWstr)
|
|
|
|
#########
|
|
minReference = referenceTop.bursts[0].burstNumber
|
|
maxReference = referenceTop.bursts[-1].burstNumber
|
|
|
|
minSecondary = secondaryTop.bursts[0].burstNumber
|
|
maxSecondary = secondaryTop.bursts[-1].burstNumber
|
|
|
|
minBurst = ifgTop.bursts[0].burstNumber
|
|
maxBurst = ifgTop.bursts[-1].burstNumber
|
|
print ('minSecondary,maxSecondary',minSecondary, maxSecondary)
|
|
print ('minReference,maxReference',minReference, maxReference)
|
|
print ('minBurst, maxBurst: ', minBurst, maxBurst)
|
|
|
|
#########
|
|
|
|
|
|
ifglist = []
|
|
factorlist = []
|
|
offsetlist = []
|
|
cohlist = []
|
|
|
|
for ii in range(minBurst, maxBurst + 1):
|
|
ind = ii - minBurst ###Index into overlaps
|
|
mind = ii - minReference ### Index into reference
|
|
sind = ii - minSecondary ###Index into secondary
|
|
|
|
topBurstIfg = ifgTop.bursts[ind]
|
|
botBurstIfg = ifgBottom.bursts[ind]
|
|
|
|
###############
|
|
'''stackReferenceTop = ifgTop.source.bursts[mind]
|
|
stackReferenceBot = ifgBottom.source.bursts[mind]
|
|
|
|
dt = stackReferenceTop.azimuthTimeInterval
|
|
topStart = int(np.round((stackReferenceBot.sensingStart - stackReferenceTop.sensingStart).total_seconds() / dt))
|
|
#overlapLen = .numberOfLines
|
|
botStart = stackReferenceBot.firstValidLine #int(np.round((.sensingStart - referenceBot.sensingStart).total_seconds() / dt))
|
|
print('+++++++++++++++++++')
|
|
print(topStart, botStart)
|
|
print('+++++++++++++++++++') '''
|
|
###############
|
|
|
|
####Double difference interferograms
|
|
topInt = np.memmap( topBurstIfg.image.filename,
|
|
dtype=np.complex64, mode='r',
|
|
shape = (topBurstIfg.numberOfLines, topBurstIfg.numberOfSamples))
|
|
|
|
botInt = np.memmap( botBurstIfg.image.filename,
|
|
dtype=np.complex64, mode='r',
|
|
shape = (botBurstIfg.numberOfLines, botBurstIfg.numberOfSamples))
|
|
|
|
intName = os.path.join(esddir, 'overlap_%02d.int'%(ii))
|
|
freqName = os.path.join(esddir, 'freq_%02d.bin'%(ii))
|
|
|
|
with open(intName, 'wb') as fid:
|
|
fid.write( topInt * np.conj(botInt))
|
|
|
|
img = isceobj.createIntImage()
|
|
img.setFilename(intName)
|
|
img.setWidth(topBurstIfg.numberOfSamples)
|
|
img.setLength(topBurstIfg.numberOfLines)
|
|
img.setAccessMode('READ')
|
|
img.renderHdr()
|
|
img.renderVRT()
|
|
img.createImage()
|
|
img.finalizeImage()
|
|
|
|
multIntName= multilook(intName, alks = inps.esdAzimuthLooks, rlks=inps.esdRangeLooks)
|
|
ifglist.append(multIntName)
|
|
|
|
|
|
####Estimate coherence of double different interferograms
|
|
multCor = createCoherence(multIntName)
|
|
cohlist.append(multCor)
|
|
|
|
####Estimate the frequency difference
|
|
azMasTop = os.path.join(referenceOffdir, 'azimuth_top_%02d_%02d.off'%(ii,ii+1))
|
|
rgMasTop = os.path.join(referenceOffdir, 'range_top_%02d_%02d.off'%(ii,ii+1))
|
|
azMasBot = os.path.join(referenceOffdir, 'azimuth_bot_%02d_%02d.off'%(ii,ii+1))
|
|
rgMasBot = os.path.join(referenceOffdir, 'range_bot_%02d_%02d.off'%(ii,ii+1))
|
|
|
|
azSlvTop = os.path.join(secondaryOffdir, 'azimuth_top_%02d_%02d.off'%(ii,ii+1))
|
|
rgSlvTop = os.path.join(secondaryOffdir, 'range_top_%02d_%02d.off'%(ii,ii+1))
|
|
azSlvBot = os.path.join(secondaryOffdir, 'azimuth_bot_%02d_%02d.off'%(ii,ii+1))
|
|
rgSlvBot = os.path.join(secondaryOffdir, 'range_bot_%02d_%02d.off'%(ii,ii+1))
|
|
|
|
mFullTop = referenceTop.source.bursts[mind]
|
|
mFullBot = referenceBot.source.bursts[mind+1]
|
|
sFullTop = secondaryTop.source.bursts[sind]
|
|
sFullBot = secondaryBot.source.bursts[sind+1]
|
|
|
|
freqdiff = overlapSpectralSeparation(topBurstIfg, botBurstIfg, mFullTop, mFullBot, sFullTop, sFullBot,
|
|
azMasTop, rgMasTop, azMasBot, rgMasBot, azSlvTop, rgSlvTop, azSlvBot, rgSlvBot)
|
|
|
|
with open(freqName, 'wb') as fid:
|
|
(freqdiff * 2 * np.pi * mFullTop.azimuthTimeInterval).astype(np.float32).tofile(fid)
|
|
|
|
img = isceobj.createImage()
|
|
img.setFilename(freqName)
|
|
img.setWidth(topBurstIfg.numberOfSamples)
|
|
img.setLength(topBurstIfg.numberOfLines)
|
|
img.setAccessMode('READ')
|
|
img.bands = 1
|
|
img.dataType = 'FLOAT'
|
|
# img.createImage()
|
|
img.renderHdr()
|
|
img.renderVRT()
|
|
img.createImage()
|
|
img.finalizeImage()
|
|
|
|
multConstName = multilook(freqName, alks = inps.esdAzimuthLooks, rlks = inps.esdRangeLooks)
|
|
factorlist.append(multConstName)
|
|
|
|
if __name__ == '__main__':
|
|
'''
|
|
Main driver.
|
|
'''
|
|
# Main Driver
|
|
main()
|
|
|
|
|
|
|