1194 lines
43 KiB
Python
1194 lines
43 KiB
Python
#!/usr/bin/env python3
|
|
|
|
#
|
|
# Author: Cunren Liang
|
|
# Copyright 2015-present, NASA-JPL/Caltech
|
|
#
|
|
|
|
|
|
def runCmd(cmd, silent=0):
|
|
import os
|
|
|
|
if silent == 0:
|
|
print("{}".format(cmd))
|
|
status = os.system(cmd)
|
|
if status != 0:
|
|
raise Exception('error when running:\n{}\n'.format(cmd))
|
|
|
|
|
|
def find_vrt_keyword(xmlfile, keyword):
|
|
from xml.etree.ElementTree import ElementTree
|
|
|
|
value = None
|
|
xmlx = ElementTree(file=open(xmlfile,'r')).getroot()
|
|
#try 10 times
|
|
for i in range(10):
|
|
path=''
|
|
for j in range(i):
|
|
path += '*/'
|
|
value0 = xmlx.find(path+keyword)
|
|
if value0 != None:
|
|
value = value0.text
|
|
break
|
|
|
|
return value
|
|
|
|
|
|
def find_vrt_file(xmlfile, keyword, relative_path=True):
|
|
'''
|
|
find file in vrt in another directory
|
|
xmlfile: vrt file
|
|
relative_path: True: return relative (to current directory) path of the file
|
|
False: return absolute path of the file
|
|
'''
|
|
import os
|
|
#get absolute directory of xmlfile
|
|
xmlfile_dir = os.path.dirname(os.path.abspath(xmlfile))
|
|
#find source file path
|
|
file = find_vrt_keyword(xmlfile, keyword)
|
|
#get absolute path of source file
|
|
file = os.path.abspath(os.path.join(xmlfile_dir, file))
|
|
#get relative path of source file
|
|
if relative_path:
|
|
file = os.path.relpath(file, './')
|
|
return file
|
|
|
|
|
|
def create_xml(fileName, width, length, fileType):
|
|
import isceobj
|
|
|
|
if fileType == 'slc':
|
|
image = isceobj.createSlcImage()
|
|
elif fileType == 'int':
|
|
image = isceobj.createIntImage()
|
|
elif fileType == 'amp':
|
|
image = isceobj.createAmpImage()
|
|
elif fileType == 'cor':
|
|
image = isceobj.createOffsetImage()
|
|
elif fileType == 'rmg' or fileType == 'unw':
|
|
image = isceobj.Image.createUnwImage()
|
|
elif fileType == 'byte':
|
|
image = isceobj.createImage()
|
|
image.setDataType('BYTE')
|
|
elif fileType == 'float':
|
|
image = isceobj.createImage()
|
|
image.setDataType('FLOAT')
|
|
elif fileType == 'double':
|
|
image = isceobj.createImage()
|
|
image.setDataType('DOUBLE')
|
|
|
|
else:
|
|
raise Exception('format not supported yet!\n')
|
|
|
|
image.setFilename(fileName)
|
|
image.extraFilename = fileName + '.vrt'
|
|
image.setWidth(width)
|
|
image.setLength(length)
|
|
|
|
#image.setAccessMode('read')
|
|
#image.createImage()
|
|
image.renderHdr()
|
|
#image.finalizeImage()
|
|
|
|
|
|
def multilook_v1(data, nalks, nrlks):
|
|
'''
|
|
doing multiple looking
|
|
ATTENSION: original array changed after running this function
|
|
'''
|
|
|
|
(length, width)=data.shape
|
|
width2 = int(width/nrlks)
|
|
length2 = int(length/nalks)
|
|
|
|
for i in range(1, nalks):
|
|
data[0:length2*nalks:nalks, :] += data[i:length2*nalks:nalks, :]
|
|
for i in range(1, nrlks):
|
|
data[0:length2*nalks:nalks, 0:width2*nrlks:nrlks] += data[0:length2*nalks:nalks, i:width2*nrlks:nrlks]
|
|
|
|
return data[0:length2*nalks:nalks, 0:width2*nrlks:nrlks] / nrlks / nalks
|
|
|
|
|
|
def multilook(data, nalks, nrlks):
|
|
'''
|
|
doing multiple looking
|
|
'''
|
|
import numpy as np
|
|
|
|
(length, width)=data.shape
|
|
width2 = int(width/nrlks)
|
|
length2 = int(length/nalks)
|
|
|
|
data2=np.zeros((length2, width), dtype=data.dtype)
|
|
for i in range(0, nalks):
|
|
data2 += data[i:length2*nalks:nalks, :]
|
|
for i in range(1, nrlks):
|
|
data2[:, 0:width2*nrlks:nrlks] += data2[:, i:width2*nrlks:nrlks]
|
|
|
|
return data2[:, 0:width2*nrlks:nrlks] / nrlks / nalks
|
|
|
|
|
|
def cal_coherence_1(inf, win=5):
|
|
'''
|
|
Compute coherence using scipy convolve 2D. Same as "def cal_coherence(inf, win=5):" in funcs.py in insarzd
|
|
|
|
#still use standard coherence estimation equation, but with magnitude removed.
|
|
#for example, equation (2) in
|
|
#H. Zebker and K. Chen, Accurate Estimation of Correlation in InSAR Observations,
|
|
#IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, VOL. 2, NO. 2, APRIL 2005.
|
|
'''
|
|
import numpy as np
|
|
import scipy.signal as ss
|
|
|
|
filt = np.ones((win,win))/ (1.0*win*win)
|
|
flag = ss.convolve2d((inf!=0), filt, mode='same')
|
|
angle = inf / (np.absolute(inf)+(inf==0))
|
|
cor = ss.convolve2d(angle, filt, mode='same')
|
|
cor = np.absolute(cor)
|
|
#remove incomplete convolution result
|
|
cor[np.nonzero(flag < 0.999)] = 0.0
|
|
#print(np.max(cor), np.min(cor))
|
|
#cor.astype(np.float32).tofile(f)
|
|
|
|
return cor
|
|
|
|
|
|
|
|
def computeOffsetFromOrbit(masterSwath, masterTrack, slaveSwath, slaveTrack, masterSample, masterLine):
|
|
'''
|
|
compute range and azimuth offsets using orbit. all range/azimuth indexes start with 0
|
|
masterSample: master sample where offset is computed, no need to be integer
|
|
masterLine: master line where offset is computed, no need to be integer
|
|
'''
|
|
import datetime
|
|
|
|
pointingDirection = {'right': -1, 'left' :1}
|
|
|
|
#compute a pair of range and azimuth offsets using geometry
|
|
#using Piyush's code for computing range and azimuth offsets
|
|
midRange = masterSwath.startingRange + masterSwath.rangePixelSize * masterSample
|
|
midSensingStart = masterSwath.sensingStart + datetime.timedelta(seconds = masterLine / masterSwath.prf)
|
|
llh = masterTrack.orbit.rdr2geo(midSensingStart, midRange, side=pointingDirection[masterTrack.pointingDirection])
|
|
slvaz, slvrng = slaveTrack.orbit.geo2rdr(llh, side=pointingDirection[masterTrack.pointingDirection])
|
|
###Translate to offsets
|
|
#at this point, slave range pixel size and prf should be the same as those of master
|
|
rgoff = ((slvrng - slaveSwath.startingRange) / masterSwath.rangePixelSize) - masterSample
|
|
azoff = ((slvaz - slaveSwath.sensingStart).total_seconds() * masterSwath.prf) - masterLine
|
|
|
|
return (rgoff, azoff)
|
|
|
|
|
|
def overlapFrequency(centerfreq1, bandwidth1, centerfreq2, bandwidth2):
|
|
startfreq1 = centerfreq1 - bandwidth1 / 2.0
|
|
endingfreq1 = centerfreq1 + bandwidth1 / 2.0
|
|
|
|
startfreq2 = centerfreq2 - bandwidth2 / 2.0
|
|
endingfreq2 = centerfreq2 + bandwidth2 / 2.0
|
|
|
|
overlapfreq = []
|
|
if startfreq2 <= startfreq1 <= endingfreq2:
|
|
overlapfreq.append(startfreq1)
|
|
if startfreq2 <= endingfreq1 <= endingfreq2:
|
|
overlapfreq.append(endingfreq1)
|
|
|
|
if startfreq1 < startfreq2 < endingfreq1:
|
|
overlapfreq.append(startfreq2)
|
|
if startfreq1 < endingfreq2 < endingfreq1:
|
|
overlapfreq.append(endingfreq2)
|
|
|
|
if len(overlapfreq) != 2:
|
|
#no overlap bandwidth
|
|
return None
|
|
else:
|
|
startfreq = min(overlapfreq)
|
|
endingfreq = max(overlapfreq)
|
|
return [startfreq, endingfreq]
|
|
|
|
|
|
def readOffset(filename):
|
|
from isceobj.Location.Offset import OffsetField,Offset
|
|
|
|
with open(filename, 'r') as f:
|
|
lines = f.readlines()
|
|
# 0 1 2 3 4 5 6 7
|
|
#retstr = "%s %s %s %s %s %s %s %s" % (self.x,self.dx,self.y,self.dy,self.snr, self.sigmax, self.sigmay, self.sigmaxy)
|
|
|
|
offsets = OffsetField()
|
|
for linex in lines:
|
|
#linexl = re.split('\s+', linex)
|
|
#detect blank lines with only spaces and tabs
|
|
if linex.strip() == '':
|
|
continue
|
|
|
|
linexl = linex.split()
|
|
offset = Offset()
|
|
#offset.setCoordinate(int(linexl[0]),int(linexl[2]))
|
|
offset.setCoordinate(float(linexl[0]),float(linexl[2]))
|
|
offset.setOffset(float(linexl[1]),float(linexl[3]))
|
|
offset.setSignalToNoise(float(linexl[4]))
|
|
offset.setCovariance(float(linexl[5]),float(linexl[6]),float(linexl[7]))
|
|
offsets.addOffset(offset)
|
|
|
|
return offsets
|
|
|
|
|
|
def writeOffset(offset, fileName):
|
|
|
|
offsetsPlain = ''
|
|
for offsetx in offset:
|
|
offsetsPlainx = "{}".format(offsetx)
|
|
offsetsPlainx = offsetsPlainx.split()
|
|
offsetsPlain = offsetsPlain + "{:8d} {:10.3f} {:8d} {:12.3f} {:11.5f} {:11.6f} {:11.6f} {:11.6f}\n".format(
|
|
int(float(offsetsPlainx[0])),
|
|
float(offsetsPlainx[1]),
|
|
int(float(offsetsPlainx[2])),
|
|
float(offsetsPlainx[3]),
|
|
float(offsetsPlainx[4]),
|
|
float(offsetsPlainx[5]),
|
|
float(offsetsPlainx[6]),
|
|
float(offsetsPlainx[7])
|
|
)
|
|
|
|
offsetFile = fileName
|
|
with open(offsetFile, 'w') as f:
|
|
f.write(offsetsPlain)
|
|
|
|
|
|
def reformatGeometricalOffset(rangeOffsetFile, azimuthOffsetFile, reformatedOffsetFile, rangeStep=1, azimuthStep=1, maximumNumberOfOffsets=10000):
|
|
'''
|
|
reformat geometrical offset as ampcor output format
|
|
'''
|
|
import numpy as np
|
|
import isceobj
|
|
|
|
img = isceobj.createImage()
|
|
img.load(rangeOffsetFile+'.xml')
|
|
width = img.width
|
|
length = img.length
|
|
|
|
step = int(np.sqrt(width*length/maximumNumberOfOffsets) + 0.5)
|
|
if step == 0:
|
|
step = 1
|
|
|
|
rgoff = np.fromfile(rangeOffsetFile, dtype=np.float32).reshape(length, width)
|
|
azoff = np.fromfile(azimuthOffsetFile, dtype=np.float32).reshape(length, width)
|
|
|
|
offsetsPlain = ''
|
|
for i in range(0, length, step):
|
|
for j in range(0, width, step):
|
|
if (rgoff[i][j] == -999999.0) or (azoff[i][j] == -999999.0):
|
|
continue
|
|
|
|
offsetsPlain = offsetsPlain + "{:8d} {:10.3f} {:8d} {:12.3f} {:11.5f} {:11.6f} {:11.6f} {:11.6f}\n".format(
|
|
int(j*rangeStep+1),
|
|
float(rgoff[i][j]),
|
|
int(i*azimuthStep+1),
|
|
float(azoff[i][j]),
|
|
float(22.00015),
|
|
float(0.000273),
|
|
float(0.002126),
|
|
float(0.000013)
|
|
)
|
|
with open(reformatedOffsetFile, 'w') as f:
|
|
f.write(offsetsPlain)
|
|
|
|
return
|
|
|
|
|
|
def cullOffsets(offsets):
|
|
import isceobj
|
|
from iscesys.StdOEL.StdOELPy import create_writer
|
|
|
|
distances = (10,5,3,3,3,3,3,3)
|
|
#numCullOffsetsLimits = (100, 75, 50, 50, 50, 50, 50, 50)
|
|
numCullOffsetsLimits = (50, 40, 30, 30, 30, 30, 30, 30)
|
|
|
|
refinedOffsets = offsets
|
|
for i, (distance, numCullOffsetsLimit) in enumerate(zip(distances, numCullOffsetsLimits)):
|
|
|
|
cullOff = isceobj.createOffoutliers()
|
|
cullOff.wireInputPort(name='offsets', object=refinedOffsets)
|
|
cullOff.setSNRThreshold(2.0)
|
|
cullOff.setDistance(distance)
|
|
|
|
#set the tag used in the outfile. each message is precided by this tag
|
|
#is the writer is not of "file" type the call has no effect
|
|
stdWriter = create_writer("log", "", True, filename="offoutliers.log")
|
|
stdWriter.setFileTag("offoutliers", "log")
|
|
stdWriter.setFileTag("offoutliers", "err")
|
|
stdWriter.setFileTag("offoutliers", "out")
|
|
cullOff.setStdWriter(stdWriter)
|
|
|
|
#run it
|
|
cullOff.offoutliers()
|
|
|
|
refinedOffsets = cullOff.getRefinedOffsetField()
|
|
numLeft = len(refinedOffsets._offsets)
|
|
print('Number of offsets left after %2dth culling: %5d'%(i, numLeft))
|
|
if numLeft < numCullOffsetsLimit:
|
|
refinedOffsets = None
|
|
|
|
stdWriter.finalize()
|
|
|
|
return refinedOffsets
|
|
|
|
|
|
def cullOffsetsRoipac(offsets, numThreshold=50):
|
|
'''
|
|
cull offsets using fortran program from ROI_PAC
|
|
numThreshold: minmum number of offsets left
|
|
'''
|
|
import os
|
|
from contrib.alos2proc_f.alos2proc_f import fitoff
|
|
from isceobj.Alos2Proc.Alos2ProcPublic import readOffset
|
|
from isceobj.Alos2Proc.Alos2ProcPublic import writeOffset
|
|
|
|
offsetFile = 'offset.off'
|
|
cullOffsetFile = 'cull.off'
|
|
writeOffset(offsets, offsetFile)
|
|
|
|
#try different parameters to cull offsets
|
|
breakFlag = 0
|
|
for maxrms in [0.08, 0.16, 0.24]:
|
|
for nsig in [1.5, 1.4, 1.3, 1.2, 1.1, 1.0, 0.9]:
|
|
fitoff(offsetFile, cullOffsetFile, nsig, maxrms, numThreshold)
|
|
|
|
#check number of matching points left
|
|
with open(cullOffsetFile, 'r') as ff:
|
|
numCullOffsets = sum(1 for linex in ff)
|
|
if numCullOffsets < numThreshold:
|
|
print('offsets culling with nsig {} maxrms {}: {} left after culling, too few points'.format(nsig, maxrms, numCullOffsets))
|
|
else:
|
|
print('offsets culling with nsig {} maxrms {}: {} left after culling, success'.format(nsig, maxrms, numCullOffsets))
|
|
breakFlag = 1
|
|
break
|
|
|
|
if breakFlag == 1:
|
|
break
|
|
|
|
if numCullOffsets < numThreshold:
|
|
refinedOffsets = None
|
|
else:
|
|
refinedOffsets = readOffset(cullOffsetFile)
|
|
|
|
os.remove(offsetFile)
|
|
os.remove(cullOffsetFile)
|
|
|
|
return refinedOffsets
|
|
|
|
|
|
def meanOffset(offsets):
|
|
|
|
rangeOffset = 0.0
|
|
azimuthOffset = 0.0
|
|
i = 0
|
|
for offsetx in offsets:
|
|
i += 1
|
|
rangeOffset += offsetx.dx
|
|
azimuthOffset += offsetx.dy
|
|
|
|
rangeOffset /= i
|
|
azimuthOffset /= i
|
|
|
|
return (rangeOffset, azimuthOffset)
|
|
|
|
|
|
def fitOffset(inputOffset, order=1, axis='range'):
|
|
'''fit a polynomial to the offset
|
|
order=0 also works, output is mean offset
|
|
'''
|
|
import numpy as np
|
|
index = []
|
|
offset = []
|
|
for a in inputOffset:
|
|
if axis=='range':
|
|
index.append(a.x)
|
|
offset.append(a.dx)
|
|
else:
|
|
index.append(a.y)
|
|
offset.append(a.dy)
|
|
|
|
p = np.polyfit(index, offset, order)
|
|
|
|
return list(p[::-1])
|
|
|
|
|
|
def topo(swath, track, demFile, latFile, lonFile, hgtFile, losFile=None, incFile=None, mskFile=None, numberRangeLooks=1, numberAzimuthLooks=1, multilookTimeOffset=True):
|
|
import datetime
|
|
import isceobj
|
|
from zerodop.topozero import createTopozero
|
|
from isceobj.Planet.Planet import Planet
|
|
|
|
pointingDirection = {'right': -1, 'left' :1}
|
|
|
|
demImage = isceobj.createDemImage()
|
|
demImage.load(demFile + '.xml')
|
|
demImage.setAccessMode('read')
|
|
|
|
#####Run Topo
|
|
planet = Planet(pname='Earth')
|
|
topo = createTopozero()
|
|
topo.slantRangePixelSpacing = numberRangeLooks * swath.rangePixelSize
|
|
topo.prf = 1.0 / (numberAzimuthLooks * swath.azimuthLineInterval)
|
|
topo.radarWavelength = track.radarWavelength
|
|
topo.orbit = track.orbit
|
|
topo.width = int(swath.numberOfSamples/numberRangeLooks)
|
|
topo.length = int(swath.numberOfLines/numberAzimuthLooks)
|
|
topo.wireInputPort(name='dem', object=demImage)
|
|
topo.wireInputPort(name='planet', object=planet)
|
|
topo.numberRangeLooks = 1 #must be set as 1
|
|
topo.numberAzimuthLooks = 1 #must be set as 1 Cunren
|
|
topo.lookSide = pointingDirection[track.pointingDirection]
|
|
if multilookTimeOffset == True:
|
|
topo.sensingStart = swath.sensingStart + datetime.timedelta(seconds=(numberAzimuthLooks-1.0)/2.0/swath.prf)
|
|
topo.rangeFirstSample = swath.startingRange + (numberRangeLooks-1.0)/2.0 * swath.rangePixelSize
|
|
else:
|
|
topo.sensingStart = swath.sensingStart
|
|
topo.rangeFirstSample = swath.startingRange
|
|
topo.demInterpolationMethod='BIQUINTIC'
|
|
|
|
topo.latFilename = latFile
|
|
topo.lonFilename = lonFile
|
|
topo.heightFilename = hgtFile
|
|
if losFile != None:
|
|
topo.losFilename = losFile
|
|
if incFile != None:
|
|
topo.incFilename = incFile
|
|
if mskFile != None:
|
|
topo.maskFilename = mskFile
|
|
|
|
topo.topo()
|
|
|
|
return list(topo.snwe)
|
|
|
|
|
|
def geo2rdr(swath, track, latFile, lonFile, hgtFile, rangeOffsetFile, azimuthOffsetFile, numberRangeLooks=1, numberAzimuthLooks=1, multilookTimeOffset=True):
|
|
import datetime
|
|
import isceobj
|
|
from zerodop.geo2rdr import createGeo2rdr
|
|
from isceobj.Planet.Planet import Planet
|
|
|
|
pointingDirection = {'right': -1, 'left' :1}
|
|
|
|
latImage = isceobj.createImage()
|
|
latImage.load(latFile + '.xml')
|
|
latImage.setAccessMode('read')
|
|
|
|
lonImage = isceobj.createImage()
|
|
lonImage.load(lonFile + '.xml')
|
|
lonImage.setAccessMode('read')
|
|
|
|
hgtImage = isceobj.createDemImage()
|
|
hgtImage.load(hgtFile + '.xml')
|
|
hgtImage.setAccessMode('read')
|
|
|
|
planet = Planet(pname='Earth')
|
|
|
|
topo = createGeo2rdr()
|
|
topo.configure()
|
|
topo.slantRangePixelSpacing = numberRangeLooks * swath.rangePixelSize
|
|
topo.prf = 1.0 / (numberAzimuthLooks * swath.azimuthLineInterval)
|
|
topo.radarWavelength = track.radarWavelength
|
|
topo.orbit = track.orbit
|
|
topo.width = int(swath.numberOfSamples/numberRangeLooks)
|
|
topo.length = int(swath.numberOfLines/numberAzimuthLooks)
|
|
topo.demLength = hgtImage.length
|
|
topo.demWidth = hgtImage.width
|
|
topo.wireInputPort(name='planet', object=planet)
|
|
topo.numberRangeLooks = 1
|
|
topo.numberAzimuthLooks = 1 #must be set to be 1
|
|
topo.lookSide = pointingDirection[track.pointingDirection]
|
|
if multilookTimeOffset == True:
|
|
topo.sensingStart = swath.sensingStart + datetime.timedelta(seconds=(numberAzimuthLooks-1.0)/2.0*swath.azimuthLineInterval)
|
|
topo.rangeFirstSample = swath.startingRange + (numberRangeLooks-1.0)/2.0*swath.rangePixelSize
|
|
else:
|
|
topo.setSensingStart(swath.sensingStart)
|
|
topo.rangeFirstSample = swath.startingRange
|
|
topo.dopplerCentroidCoeffs = [0.] #we are using zero doppler geometry
|
|
topo.demImage = hgtImage
|
|
topo.latImage = latImage
|
|
topo.lonImage = lonImage
|
|
topo.rangeOffsetImageName = rangeOffsetFile
|
|
topo.azimuthOffsetImageName = azimuthOffsetFile
|
|
topo.geo2rdr()
|
|
|
|
return
|
|
|
|
|
|
def waterBodyRadar(latFile, lonFile, wbdFile, wbdOutFile):
|
|
'''
|
|
create water boday in radar coordinates
|
|
'''
|
|
import numpy as np
|
|
import isceobj
|
|
from isceobj.Alos2Proc.Alos2ProcPublic import create_xml
|
|
|
|
demImage = isceobj.createDemImage()
|
|
demImage.load(wbdFile + '.xml')
|
|
#demImage.setAccessMode('read')
|
|
wbd=np.memmap(wbdFile, dtype='byte', mode='r', shape=(demImage.length, demImage.width))
|
|
|
|
image = isceobj.createImage()
|
|
image.load(latFile+'.xml')
|
|
width = image.width
|
|
length = image.length
|
|
|
|
latFp = open(latFile, 'rb')
|
|
lonFp = open(lonFile, 'rb')
|
|
wbdOutFp = open(wbdOutFile, 'wb')
|
|
print("create water body in radar coordinates...")
|
|
for i in range(length):
|
|
if (((i+1)%200) == 0):
|
|
print("processing line %6d of %6d" % (i+1, length), end='\r', flush=True)
|
|
wbdOut = np.zeros(width, dtype='byte')-2
|
|
lat = np.fromfile(latFp, dtype=np.float64, count=width)
|
|
lon = np.fromfile(lonFp, dtype=np.float64, count=width)
|
|
#indexes start with zero
|
|
lineIndex = np.int32((lat - demImage.firstLatitude) / demImage.deltaLatitude + 0.5)
|
|
sampleIndex = np.int32((lon - demImage.firstLongitude) / demImage.deltaLongitude + 0.5)
|
|
inboundIndex = np.logical_and(
|
|
np.logical_and(lineIndex>=0, lineIndex<=demImage.length-1),
|
|
np.logical_and(sampleIndex>=0, sampleIndex<=demImage.width-1)
|
|
)
|
|
#keep SRTM convention. water body. (0) --- land; (-1) --- water; (-2 or other value) --- no data.
|
|
wbdOut = wbd[(lineIndex[inboundIndex], sampleIndex[inboundIndex])]
|
|
wbdOut.astype(np.int8).tofile(wbdOutFp)
|
|
print("processing line %6d of %6d" % (length, length))
|
|
#create_xml(wbdOutFile, width, length, 'byte')
|
|
|
|
image = isceobj.createImage()
|
|
image.setDataType('BYTE')
|
|
image.addDescription('water body. (0) --- land; (-1) --- water; (-2) --- no data.')
|
|
image.setFilename(wbdOutFile)
|
|
image.extraFilename = wbdOutFile + '.vrt'
|
|
image.setWidth(width)
|
|
image.setLength(length)
|
|
image.renderHdr()
|
|
|
|
del wbd, demImage, image
|
|
latFp.close()
|
|
lonFp.close()
|
|
wbdOutFp.close()
|
|
|
|
|
|
def renameFile(oldname, newname):
|
|
import os
|
|
import isceobj
|
|
img = isceobj.createImage()
|
|
img.load(oldname + '.xml')
|
|
img.setFilename(newname)
|
|
img.extraFilename = newname+'.vrt'
|
|
img.renderHdr()
|
|
|
|
os.rename(oldname, newname)
|
|
os.remove(oldname + '.xml')
|
|
os.remove(oldname + '.vrt')
|
|
|
|
|
|
def cal_coherence(inf, win=5, edge=0):
|
|
'''
|
|
compute coherence uisng only interferogram (phase).
|
|
This routine still follows the regular equation for computing coherence,
|
|
but assumes the amplitudes of master and slave are one, so that coherence
|
|
can be computed using phase only.
|
|
|
|
inf: interferogram
|
|
win: window size
|
|
edge: 0: remove all non-full convolution samples
|
|
|
|
1: remove samples computed from less than half convolution
|
|
(win=5 used to illustration below)
|
|
* * *
|
|
* * *
|
|
* * *
|
|
* * *
|
|
* * *
|
|
|
|
2: remove samples computed from less than quater convolution
|
|
(win=5 used to illustration below)
|
|
* * *
|
|
* * *
|
|
* * *
|
|
|
|
3: remove non-full convolution samples on image edges
|
|
|
|
4: keep all samples
|
|
'''
|
|
import numpy as np
|
|
import scipy.signal as ss
|
|
|
|
if win % 2 != 1:
|
|
raise Exception('window size must be odd!')
|
|
hwin = np.int(np.around((win - 1) / 2))
|
|
|
|
filt = np.ones((win, win))
|
|
amp = np.absolute(inf)
|
|
|
|
cnt = ss.convolve2d((amp!=0), filt, mode='same')
|
|
cor = ss.convolve2d(inf/(amp + (amp==0)), filt, mode='same')
|
|
cor = (amp!=0) * np.absolute(cor) / (cnt + (cnt==0))
|
|
|
|
#trim edges
|
|
if edge == 0:
|
|
num = win * win
|
|
cor[np.nonzero(cnt < num)] = 0.0
|
|
elif edge == 1:
|
|
num = win * (hwin+1)
|
|
cor[np.nonzero(cnt < num)] = 0.0
|
|
elif edge == 2:
|
|
num = (hwin+1) * (hwin+1)
|
|
cor[np.nonzero(cnt < num)] = 0.0
|
|
elif edge == 3:
|
|
cor[0:hwin, :] = 0.0
|
|
cor[-hwin:, :] = 0.0
|
|
cor[:, 0:hwin] = 0.0
|
|
cor[:, -hwin:] = 0.0
|
|
else:
|
|
pass
|
|
|
|
#print("coherence, max: {} min: {}".format(np.max(cor[np.nonzero(cor!=0)]), np.min(cor[np.nonzero(cor!=0)])))
|
|
return cor
|
|
|
|
|
|
def snaphuUnwrap(track, t, wrapName, corName, unwrapName, nrlks, nalks, costMode = 'DEFO',initMethod = 'MST', defomax = 4.0, initOnly = False):
|
|
#runUnwrap(self, costMode = 'SMOOTH',initMethod = 'MCF', defomax = 2, initOnly = True)
|
|
'''
|
|
track: track object
|
|
t: time for computing earth radius and altitude, normally mid azimuth time
|
|
wrapName: input interferogram
|
|
corName: input coherence file
|
|
unwrapName: output unwrapped interferogram
|
|
nrlks: number of range looks of the interferogram
|
|
nalks: number of azimuth looks of the interferogram
|
|
'''
|
|
import datetime
|
|
import numpy as np
|
|
import isceobj
|
|
from contrib.Snaphu.Snaphu import Snaphu
|
|
from isceobj.Planet.Planet import Planet
|
|
|
|
corImg = isceobj.createImage()
|
|
corImg.load(corName + '.xml')
|
|
width = corImg.width
|
|
length = corImg.length
|
|
|
|
#get altitude
|
|
orbit = track.orbit
|
|
peg = orbit.interpolateOrbit(t, method='hermite')
|
|
refElp = Planet(pname='Earth').ellipsoid
|
|
llh = refElp.xyz_to_llh(peg.getPosition())
|
|
hdg = orbit.getENUHeading(t)
|
|
refElp.setSCH(llh[0], llh[1], hdg)
|
|
earthRadius = refElp.pegRadCur
|
|
altitude = llh[2]
|
|
|
|
rangeLooks = nrlks
|
|
azimuthLooks = nalks
|
|
azfact = 0.8
|
|
rngfact = 0.8
|
|
corrLooks = rangeLooks * azimuthLooks/(azfact*rngfact)
|
|
maxComponents = 20
|
|
|
|
snp = Snaphu()
|
|
snp.setInitOnly(initOnly)
|
|
snp.setInput(wrapName)
|
|
snp.setOutput(unwrapName)
|
|
snp.setWidth(width)
|
|
snp.setCostMode(costMode)
|
|
snp.setEarthRadius(earthRadius)
|
|
snp.setWavelength(track.radarWavelength)
|
|
snp.setAltitude(altitude)
|
|
snp.setCorrfile(corName)
|
|
snp.setInitMethod(initMethod)
|
|
snp.setCorrLooks(corrLooks)
|
|
snp.setMaxComponents(maxComponents)
|
|
snp.setDefoMaxCycles(defomax)
|
|
snp.setRangeLooks(rangeLooks)
|
|
snp.setAzimuthLooks(azimuthLooks)
|
|
if corImg.bands == 1:
|
|
snp.setCorFileFormat('FLOAT_DATA')
|
|
snp.prepare()
|
|
snp.unwrap()
|
|
|
|
######Render XML
|
|
outImage = isceobj.Image.createUnwImage()
|
|
outImage.setFilename(unwrapName)
|
|
outImage.setWidth(width)
|
|
outImage.setAccessMode('read')
|
|
outImage.renderVRT()
|
|
outImage.createImage()
|
|
outImage.finalizeImage()
|
|
outImage.renderHdr()
|
|
|
|
#####Check if connected components was created
|
|
if snp.dumpConnectedComponents:
|
|
connImage = isceobj.Image.createImage()
|
|
connImage.setFilename(unwrapName+'.conncomp')
|
|
connImage.setWidth(width)
|
|
connImage.setAccessMode('read')
|
|
connImage.setDataType('BYTE')
|
|
connImage.renderVRT()
|
|
connImage.createImage()
|
|
connImage.finalizeImage()
|
|
connImage.renderHdr()
|
|
del connImage
|
|
|
|
del corImg
|
|
del snp
|
|
del outImage
|
|
|
|
#remove wired things in no-data area
|
|
amp=np.memmap(unwrapName, dtype='float32', mode='r+', shape=(length*2, width))
|
|
wrap = np.fromfile(wrapName, dtype=np.complex64).reshape(length, width)
|
|
(amp[0:length*2:2, :])[np.nonzero(wrap==0)]=0
|
|
(amp[1:length*2:2, :])[np.nonzero(wrap==0)]=0
|
|
del amp
|
|
del wrap
|
|
|
|
return
|
|
|
|
|
|
def snaphuUnwrapOriginal(wrapName, corName, ampName, unwrapName, costMode = 's', initMethod = 'mcf'):
|
|
'''
|
|
unwrap interferogram using original snaphu program
|
|
'''
|
|
import numpy as np
|
|
import isceobj
|
|
|
|
corImg = isceobj.createImage()
|
|
corImg.load(corName + '.xml')
|
|
width = corImg.width
|
|
length = corImg.length
|
|
|
|
#specify coherence file format in configure file
|
|
snaphuConfFile = 'snaphu.conf'
|
|
if corImg.bands == 1:
|
|
snaphuConf = '''CORRFILEFORMAT FLOAT_DATA
|
|
CONNCOMPFILE {}
|
|
MAXNCOMPS 20'''.format(unwrapName+'.conncomp')
|
|
|
|
else:
|
|
snaphuConf = '''CORRFILEFORMAT FLOAT_DATA
|
|
CONNCOMPFILE {}
|
|
MAXNCOMPS 20'''.format(unwrapName+'.conncomp')
|
|
with open(snaphuConfFile, 'w') as f:
|
|
f.write(snaphuConf)
|
|
cmd = 'snaphu {} {} -f {} -{} -o {} -a {} -c {} -v --{}'.format(
|
|
wrapName,
|
|
width,
|
|
snaphuConfFile,
|
|
costMode,
|
|
unwrapName,
|
|
ampName,
|
|
corName,
|
|
initMethod
|
|
)
|
|
runCmd(cmd)
|
|
create_xml(unwrapName, width, length, 'unw')
|
|
|
|
connImage = isceobj.Image.createImage()
|
|
connImage.setFilename(unwrapName+'.conncomp')
|
|
connImage.setWidth(width)
|
|
connImage.setAccessMode('read')
|
|
connImage.setDataType('BYTE')
|
|
connImage.renderVRT()
|
|
connImage.createImage()
|
|
connImage.finalizeImage()
|
|
connImage.renderHdr()
|
|
del connImage
|
|
|
|
#remove wired things in no-data area
|
|
amp=np.memmap(unwrapName, dtype='float32', mode='r+', shape=(length*2, width))
|
|
wrap = np.fromfile(wrapName, dtype=np.complex64).reshape(length, width)
|
|
(amp[0:length*2:2, :])[np.nonzero(wrap==0)]=0
|
|
(amp[1:length*2:2, :])[np.nonzero(wrap==0)]=0
|
|
del amp
|
|
del wrap
|
|
|
|
return
|
|
|
|
|
|
def getBboxGeo(track):
|
|
'''
|
|
get bounding box in geo-coordinate
|
|
'''
|
|
import numpy as np
|
|
|
|
pointingDirection = {'right': -1, 'left' :1}
|
|
|
|
bboxRdr = getBboxRdr(track)
|
|
|
|
rangeMin = bboxRdr[0]
|
|
rangeMax = bboxRdr[1]
|
|
azimuthTimeMin = bboxRdr[2]
|
|
azimuthTimeMax = bboxRdr[3]
|
|
|
|
#get bounding box using Piyush's code
|
|
hgtrange=[-500,9000]
|
|
ts = [azimuthTimeMin, azimuthTimeMax]
|
|
rngs = [rangeMin, rangeMax]
|
|
pos = []
|
|
for ht in hgtrange:
|
|
for tim in ts:
|
|
for rng in rngs:
|
|
llh = track.orbit.rdr2geo(tim, rng, height=ht, side=pointingDirection[track.pointingDirection])
|
|
pos.append(llh)
|
|
pos = np.array(pos)
|
|
# S N W E
|
|
bbox = [np.min(pos[:,0]), np.max(pos[:,0]), np.min(pos[:,1]), np.max(pos[:,1])]
|
|
|
|
return bbox
|
|
|
|
|
|
def getBboxRdr(track):
|
|
'''
|
|
get bounding box in radar-coordinate
|
|
'''
|
|
import datetime
|
|
|
|
numberOfFrames = len(track.frames)
|
|
numberOfSwaths = len(track.frames[0].swaths)
|
|
|
|
sensingStartList = []
|
|
sensingEndList = []
|
|
startingRangeList = []
|
|
endingRangeList = []
|
|
for i in range(numberOfFrames):
|
|
for j in range(numberOfSwaths):
|
|
swath = track.frames[i].swaths[j]
|
|
sensingStartList.append(swath.sensingStart)
|
|
sensingEndList.append(swath.sensingStart + datetime.timedelta(seconds=(swath.numberOfLines-1) * swath.azimuthLineInterval))
|
|
startingRangeList.append(swath.startingRange)
|
|
endingRangeList.append(swath.startingRange + (swath.numberOfSamples - 1) * swath.rangePixelSize)
|
|
azimuthTimeMin = min(sensingStartList)
|
|
azimuthTimeMax = max(sensingEndList)
|
|
azimuthTimeMid = azimuthTimeMin+datetime.timedelta(seconds=(azimuthTimeMax-azimuthTimeMin).total_seconds()/2.0)
|
|
rangeMin = min(startingRangeList)
|
|
rangeMax = max(endingRangeList)
|
|
rangeMid = (rangeMin + rangeMax) / 2.0
|
|
|
|
bbox = [rangeMin, rangeMax, azimuthTimeMin, azimuthTimeMax]
|
|
|
|
return bbox
|
|
|
|
|
|
def filterInterferogram(data, alpha, windowSize, stepSize):
|
|
'''
|
|
a filter wrapper
|
|
'''
|
|
import os
|
|
import numpy as np
|
|
from contrib.alos2filter.alos2filter import psfilt1
|
|
|
|
(length, width)=data.shape
|
|
data.astype(np.complex64).tofile('tmp1234.int')
|
|
psfilt1('tmp1234.int', 'filt_tmp1234.int', width, alpha, windowSize, stepSize)
|
|
|
|
data2 = np.fromfile('filt_tmp1234.int', dtype=np.complex64).reshape(length, width)
|
|
os.remove('tmp1234.int')
|
|
os.remove('filt_tmp1234.int')
|
|
|
|
return data2
|
|
|
|
|
|
|
|
###################################################################
|
|
# these are routines for burst-by-burst ScanSAR interferometry
|
|
###################################################################
|
|
|
|
def mosaicBurstInterferogram(swath, burstPrefix, outputFile, numberOfLooksThreshold=1):
|
|
'''
|
|
take a burst sequence and output mosaicked file
|
|
'''
|
|
import numpy as np
|
|
|
|
interferogram = np.zeros((swath.numberOfLines, swath.numberOfSamples), dtype=np.complex64)
|
|
cnt = np.zeros((swath.numberOfLines, swath.numberOfSamples), dtype=np.int8)
|
|
for i in range(swath.numberOfBursts):
|
|
burstFile = burstPrefix + '_%02d.int'%(i+1)
|
|
burstInterferogram = np.fromfile(burstFile, dtype=np.complex64).reshape(swath.burstSlcNumberOfLines, swath.burstSlcNumberOfSamples)
|
|
interferogram[0+swath.burstSlcFirstLineOffsets[i]:swath.burstSlcNumberOfLines+swath.burstSlcFirstLineOffsets[i], :] += burstInterferogram
|
|
cnt[0+swath.burstSlcFirstLineOffsets[i]:swath.burstSlcNumberOfLines+swath.burstSlcFirstLineOffsets[i], :] += (burstInterferogram!=0)
|
|
|
|
#trim upper and lower edges with less number of looks
|
|
#############################################################################
|
|
firstLine = 0
|
|
for i in range(swath.numberOfLines):
|
|
if np.sum(cnt[i,:]>=numberOfLooksThreshold) > swath.numberOfSamples/2:
|
|
firstLine = i
|
|
break
|
|
lastLine = swath.numberOfLines - 1
|
|
for i in range(swath.numberOfLines):
|
|
if np.sum(cnt[swath.numberOfLines-1-i,:]>=numberOfLooksThreshold) > swath.numberOfSamples/2:
|
|
lastLine = swath.numberOfLines-1-i
|
|
break
|
|
|
|
interferogram[:firstLine,:]=0
|
|
interferogram[lastLine+1:,:]=0
|
|
|
|
# if numberOfLooksThreshold!= None:
|
|
# interferogram[np.nonzero(cnt<numberOfLooksThreshold)] = 0
|
|
#############################################################################
|
|
|
|
interferogram.astype(np.complex64).tofile(outputFile)
|
|
create_xml(outputFile, swath.numberOfSamples, swath.numberOfLines, 'int')
|
|
|
|
|
|
def mosaicBurstAmplitude(swath, burstPrefix, outputFile, numberOfLooksThreshold=1):
|
|
'''
|
|
take a burst sequence and output the magnitude
|
|
'''
|
|
import numpy as np
|
|
|
|
amp = np.zeros((swath.numberOfLines, swath.numberOfSamples), dtype=np.float32)
|
|
cnt = np.zeros((swath.numberOfLines, swath.numberOfSamples), dtype=np.int8)
|
|
for i in range(swath.numberOfBursts):
|
|
burstFile = burstPrefix + '_%02d.slc'%(i+1)
|
|
#azLineOffset = round((swath.burstSlcStartTimes[i] - swath.burstSlcStartTimes[0]).total_seconds() / swath.azimuthLineInterval)
|
|
burstMag = np.absolute(np.fromfile(burstFile, dtype=np.complex64).reshape(swath.burstSlcNumberOfLines, swath.burstSlcNumberOfSamples))
|
|
burstPwr = burstMag * burstMag
|
|
amp[0+swath.burstSlcFirstLineOffsets[i]:swath.burstSlcNumberOfLines+swath.burstSlcFirstLineOffsets[i], :] += burstPwr
|
|
cnt[0+swath.burstSlcFirstLineOffsets[i]:swath.burstSlcNumberOfLines+swath.burstSlcFirstLineOffsets[i], :] += (burstPwr!=0)
|
|
|
|
#trim upper and lower edges with less number of looks
|
|
#############################################################################
|
|
firstLine = 0
|
|
for i in range(swath.numberOfLines):
|
|
if np.sum(cnt[i,:]>=numberOfLooksThreshold) > swath.numberOfSamples/2:
|
|
firstLine = i
|
|
break
|
|
lastLine = swath.numberOfLines - 1
|
|
for i in range(swath.numberOfLines):
|
|
if np.sum(cnt[swath.numberOfLines-1-i,:]>=numberOfLooksThreshold) > swath.numberOfSamples/2:
|
|
lastLine = swath.numberOfLines-1-i
|
|
break
|
|
|
|
amp[:firstLine,:]=0
|
|
amp[lastLine+1:,:]=0
|
|
|
|
# if numberOfLooksThreshold!= None:
|
|
# amp[np.nonzero(cnt<numberOfLooksThreshold)] = 0
|
|
#############################################################################
|
|
|
|
np.sqrt(amp).astype(np.float32).tofile(outputFile)
|
|
create_xml(outputFile, swath.numberOfSamples, swath.numberOfLines, 'float')
|
|
|
|
|
|
def resampleBursts(masterSwath, slaveSwath,
|
|
masterBurstDir, slaveBurstDir, slaveBurstResampledDir, interferogramDir,
|
|
masterBurstPrefix, slaveBurstPrefix, slaveBurstResampledPrefix, interferogramPrefix,
|
|
rangeOffset, azimuthOffset, rangeOffsetResidual=0, azimuthOffsetResidual=0):
|
|
|
|
import os
|
|
import datetime
|
|
import numpy as np
|
|
import numpy.matlib
|
|
from contrib.alos2proc.alos2proc import resamp
|
|
|
|
os.makedirs(slaveBurstResampledDir, exist_ok=True)
|
|
os.makedirs(interferogramDir, exist_ok=True)
|
|
|
|
#get burst file names
|
|
masterBurstSlc = [masterBurstPrefix+'_%02d.slc'%(i+1) for i in range(masterSwath.numberOfBursts)]
|
|
slaveBurstSlc = [slaveBurstPrefix+'_%02d.slc'%(i+1) for i in range(slaveSwath.numberOfBursts)]
|
|
slaveBurstSlcResampled = [slaveBurstPrefix+'_%02d.slc'%(i+1) for i in range(masterSwath.numberOfBursts)]
|
|
interferogram = [interferogramPrefix+'_%02d.int'%(i+1) for i in range(masterSwath.numberOfBursts)]
|
|
|
|
length = masterSwath.burstSlcNumberOfLines
|
|
width = masterSwath.burstSlcNumberOfSamples
|
|
lengthSlave = slaveSwath.burstSlcNumberOfLines
|
|
widthSlave = slaveSwath.burstSlcNumberOfSamples
|
|
|
|
#slave burst slc start times
|
|
slaveBurstStartTimesSlc = [slaveSwath.firstBurstSlcStartTime + \
|
|
datetime.timedelta(seconds=slaveSwath.burstSlcFirstLineOffsets[i]*slaveSwath.azimuthLineInterval) \
|
|
for i in range(slaveSwath.numberOfBursts)]
|
|
#slave burst raw start times
|
|
slaveBurstStartTimesRaw = [slaveSwath.firstBurstRawStartTime + \
|
|
datetime.timedelta(seconds=i*slaveSwath.burstCycleLength/slaveSwath.prf) \
|
|
for i in range(slaveSwath.numberOfBursts)]
|
|
|
|
|
|
for i in range(masterSwath.numberOfBursts):
|
|
|
|
##########################################################################
|
|
# 1. get offsets and corresponding slave burst
|
|
##########################################################################
|
|
#range offset
|
|
with open(rangeOffset, 'rb') as f:
|
|
f.seek(masterSwath.burstSlcFirstLineOffsets[i] * width * np.dtype(np.float32).itemsize, 0)
|
|
rgoffBurst = np.fromfile(f, dtype=np.float32, count=length*width).reshape(length,width)
|
|
if type(rangeOffsetResidual) == np.ndarray:
|
|
residual = rangeOffsetResidual[0+masterSwath.burstSlcFirstLineOffsets[i]:length+masterSwath.burstSlcFirstLineOffsets[i],:]
|
|
rgoffBurst[np.nonzero(rgoffBurst!=-999999.0)] += residual[np.nonzero(rgoffBurst!=-999999.0)]
|
|
else:
|
|
rgoffBurst[np.nonzero(rgoffBurst!=-999999.0)] += rangeOffsetResidual
|
|
#azimuth offset
|
|
with open(azimuthOffset, 'rb') as f:
|
|
f.seek(masterSwath.burstSlcFirstLineOffsets[i] * width * np.dtype(np.float32).itemsize, 0)
|
|
azoffBurst = np.fromfile(f, dtype=np.float32, count=length*width).reshape(length,width)
|
|
if type(azimuthOffsetResidual) == np.ndarray:
|
|
residual = azimuthOffsetResidual[0+masterSwath.burstSlcFirstLineOffsets[i]:length+masterSwath.burstSlcFirstLineOffsets[i],:]
|
|
azoffBurst[np.nonzero(azoffBurst!=-999999.0)] += residual[np.nonzero(azoffBurst!=-999999.0)]
|
|
else:
|
|
azoffBurst[np.nonzero(azoffBurst!=-999999.0)] += azimuthOffsetResidual
|
|
|
|
#find the corresponding slave burst
|
|
#get mean offset to use
|
|
#remove BAD_VALUE = -999999.0 as defined in geo2rdr.f90
|
|
#single precision is not accurate enough to compute mean
|
|
azoffBurstMean = np.mean(azoffBurst[np.nonzero(azoffBurst!=-999999.0)], dtype=np.float64)
|
|
iSlave = -1
|
|
for j in range(slaveSwath.numberOfBursts):
|
|
if abs(masterSwath.burstSlcFirstLineOffsets[i] + azoffBurstMean - slaveSwath.burstSlcFirstLineOffsets[j]) < (masterSwath.burstLength / masterSwath.prf * 2.0) / masterSwath.azimuthLineInterval:
|
|
iSlave = j
|
|
break
|
|
|
|
#output zero resampled burst/interferogram if no slave burst found
|
|
if iSlave == -1:
|
|
print('\nburst pair, master: %2d, slave: no'%(i+1))
|
|
#output an interferogram with all pixels set to zero
|
|
os.chdir(interferogramDir)
|
|
np.zeros((length, width), dtype=np.complex64).astype(np.complex64).tofile(interferogram[i])
|
|
create_xml(interferogram[i], width, length, 'int')
|
|
os.chdir('../')
|
|
#output a resampled slave image with all pixels set to zero
|
|
os.chdir(slaveBurstResampledDir)
|
|
np.zeros((length, width), dtype=np.complex64).astype(np.complex64).tofile(slaveBurstSlcResampled[i])
|
|
create_xml(slaveBurstSlcResampled[i], width, length, 'slc')
|
|
os.chdir('../')
|
|
continue
|
|
else:
|
|
print('\nburst pair, master: %2d, slave: %3d'%(i+1, iSlave+1))
|
|
|
|
#adjust azimuth offset accordingly, since original azimuth offset assumes master and slave start with sensingStart
|
|
azoffBurst -= (slaveSwath.burstSlcFirstLineOffsets[iSlave]-masterSwath.burstSlcFirstLineOffsets[i])
|
|
|
|
|
|
##########################################################################
|
|
# 2. compute deramp and reramp signals
|
|
##########################################################################
|
|
cj = np.complex64(1j)
|
|
tbase = (slaveBurstStartTimesSlc[iSlave] - (slaveBurstStartTimesRaw[iSlave] + \
|
|
datetime.timedelta(seconds=(slaveSwath.burstLength - 1.0) / 2.0 / slaveSwath.prf))).total_seconds()
|
|
|
|
#compute deramp signal
|
|
index1 = np.matlib.repmat(np.arange(widthSlave), lengthSlave, 1)
|
|
index2 = np.matlib.repmat(np.arange(lengthSlave).reshape(lengthSlave, 1), 1, widthSlave)
|
|
ka = slaveSwath.azimuthFmrateVsPixel[3] * index1**3 + slaveSwath.azimuthFmrateVsPixel[2] * index1**2 + \
|
|
slaveSwath.azimuthFmrateVsPixel[1] * index1 + slaveSwath.azimuthFmrateVsPixel[0]
|
|
#use the convention that ka > 0
|
|
ka = -ka
|
|
t = tbase + index2*slaveSwath.azimuthLineInterval
|
|
deramp = np.exp(cj * np.pi * (-ka) * t**2)
|
|
|
|
#compute reramp signal
|
|
index1 = np.matlib.repmat(np.arange(width), length, 1) + rgoffBurst
|
|
index2 = np.matlib.repmat(np.arange(length).reshape(length, 1), 1, width) + azoffBurst
|
|
ka = slaveSwath.azimuthFmrateVsPixel[3] * index1**3 + slaveSwath.azimuthFmrateVsPixel[2] * index1**2 + \
|
|
slaveSwath.azimuthFmrateVsPixel[1] * index1 + slaveSwath.azimuthFmrateVsPixel[0]
|
|
#use the convention that ka > 0
|
|
ka = -ka
|
|
t = tbase + index2*slaveSwath.azimuthLineInterval
|
|
reramp = np.exp(cj * np.pi * (ka) * t**2)
|
|
|
|
|
|
##########################################################################
|
|
# 3. resample slave burst
|
|
##########################################################################
|
|
#go to slave directory to do resampling
|
|
os.chdir(slaveBurstDir)
|
|
|
|
#output offsets
|
|
rgoffBurstFile = "burst_rg.off"
|
|
azoffBurstFile = "burst_az.off"
|
|
rgoffBurst.astype(np.float32).tofile(rgoffBurstFile)
|
|
azoffBurst.astype(np.float32).tofile(azoffBurstFile)
|
|
|
|
#deramp slave burst
|
|
slaveBurstDerampedFile = "slave.slc"
|
|
sburst = np.fromfile(slaveBurstSlc[iSlave], dtype=np.complex64).reshape(lengthSlave, widthSlave)
|
|
(deramp * sburst).astype(np.complex64).tofile(slaveBurstDerampedFile)
|
|
create_xml(slaveBurstDerampedFile, widthSlave, lengthSlave, 'slc')
|
|
|
|
#resampled slave burst
|
|
slaveBurstResampFile = 'slave_resamp.slc'
|
|
|
|
#resample slave burst
|
|
#now doppler has bigger impact now, as it's value is about 35 Hz (azimuth resampling frequency is now only 1/20 * PRF)
|
|
#we don't know if this doppler value is accurate or not, so we set it to zero, which seems to give best resampling result
|
|
#otherwise if it is not accurate and we still use it, it will significantly affect resampling result
|
|
dopplerVsPixel = slaveSwath.dopplerVsPixel
|
|
dopplerVsPixel = [0.0, 0.0, 0.0, 0.0]
|
|
|
|
resamp(slaveBurstDerampedFile, slaveBurstResampFile, rgoffBurstFile, azoffBurstFile, width, length, 1.0/slaveSwath.azimuthLineInterval, dopplerVsPixel,
|
|
rgcoef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
|
azcoef=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
|
|
azpos_off=0.0)
|
|
|
|
#read resampled slave burst and reramp
|
|
sburstResamp = reramp * (np.fromfile(slaveBurstResampFile, dtype=np.complex64).reshape(length, width))
|
|
|
|
#clear up
|
|
os.remove(rgoffBurstFile)
|
|
os.remove(azoffBurstFile)
|
|
os.remove(slaveBurstDerampedFile)
|
|
os.remove(slaveBurstDerampedFile+'.vrt')
|
|
os.remove(slaveBurstDerampedFile+'.xml')
|
|
os.remove(slaveBurstResampFile)
|
|
os.remove(slaveBurstResampFile+'.vrt')
|
|
os.remove(slaveBurstResampFile+'.xml')
|
|
|
|
os.chdir('../')
|
|
|
|
|
|
##########################################################################
|
|
# 4. dump results
|
|
##########################################################################
|
|
#dump resampled slave burst
|
|
os.chdir(slaveBurstResampledDir)
|
|
sburstResamp.astype(np.complex64).tofile(slaveBurstSlcResampled[i])
|
|
create_xml(slaveBurstSlcResampled[i], width, length, 'slc')
|
|
os.chdir('../')
|
|
|
|
#dump burst interferogram
|
|
mburst = np.fromfile(os.path.join(masterBurstDir, masterBurstSlc[i]), dtype=np.complex64).reshape(length, width)
|
|
os.chdir(interferogramDir)
|
|
(mburst * np.conj(sburstResamp)).astype(np.complex64).tofile(interferogram[i])
|
|
create_xml(interferogram[i], width, length, 'int')
|
|
os.chdir('../')
|
|
|
|
|
|
def create_multi_index(width, rgl):
|
|
import numpy as np
|
|
#create index after multilooking
|
|
#assuming original index start with 0
|
|
#applies to both range and azimuth direction
|
|
|
|
widthm = int(width/rgl)
|
|
|
|
#create range index: This applies to both odd and even cases, "rgl = 1" case, and "rgl = 2" case
|
|
start_rgindex = (rgl - 1.0) / 2.0
|
|
rgindex0 = start_rgindex + np.arange(widthm) * rgl
|
|
|
|
return rgindex0
|
|
|
|
|
|
def create_multi_index2(width2, l1, l2):
|
|
import numpy as np
|
|
#for number of looks of l1 and l2
|
|
#calculate the correponding index number of l2 in the l1 array
|
|
#applies to both range and azimuth direction
|
|
|
|
return ((l2 - l1) / 2.0 + np.arange(width2) * l2) / l1
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|