diff --git a/components/isceobj/Alos2Proc/runGeocode.py b/components/isceobj/Alos2Proc/runGeocode.py index 9dfbbbf..386bb38 100644 --- a/components/isceobj/Alos2Proc/runGeocode.py +++ b/components/isceobj/Alos2Proc/runGeocode.py @@ -4,6 +4,7 @@ # import os +import glob import logging import numpy as np @@ -42,7 +43,9 @@ def runGeocode(self): if self.doIon: geocodeList.append(self._insar.multilookIon) else: - geocodeList = self.geocodeList + geocodeList = [] + for xxx in self.geocodeList: + geocodeList += glob.glob(xxx) numberRangeLooks = self._insar.numberRangeLooks1 * self._insar.numberRangeLooks2 numberAzimuthLooks = self._insar.numberAzimuthLooks1 * self._insar.numberAzimuthLooks2 diff --git a/components/isceobj/Alos2Proc/runIonSubband.py b/components/isceobj/Alos2Proc/runIonSubband.py index ef95806..ce7d16f 100644 --- a/components/isceobj/Alos2Proc/runIonSubband.py +++ b/components/isceobj/Alos2Proc/runIonSubband.py @@ -88,6 +88,23 @@ def runIonSubband(self): frameDir = 'f{}_{}'.format(i+1, frameNumber) for j, swathNumber in enumerate(range(self._insar.startingSwath, self._insar.endingSwath + 1)): swathDir = 's{}'.format(swathNumber) + + #skip this time consuming process, if interferogram already exists + if os.path.isfile(os.path.join(ionDir['subband'][0], frameDir, swathDir, self._insar.interferogram)) and \ + os.path.isfile(os.path.join(ionDir['subband'][0], frameDir, swathDir, self._insar.interferogram+'.vrt')) and \ + os.path.isfile(os.path.join(ionDir['subband'][0], frameDir, swathDir, self._insar.interferogram+'.xml')) and \ + os.path.isfile(os.path.join(ionDir['subband'][0], frameDir, swathDir, self._insar.amplitude)) and \ + os.path.isfile(os.path.join(ionDir['subband'][0], frameDir, swathDir, self._insar.amplitude+'.vrt')) and \ + os.path.isfile(os.path.join(ionDir['subband'][0], frameDir, swathDir, self._insar.amplitude+'.xml')) and \ + os.path.isfile(os.path.join(ionDir['subband'][1], frameDir, swathDir, self._insar.interferogram)) and \ + os.path.isfile(os.path.join(ionDir['subband'][1], frameDir, swathDir, self._insar.interferogram+'.vrt')) and \ + os.path.isfile(os.path.join(ionDir['subband'][1], frameDir, swathDir, self._insar.interferogram+'.xml')) and \ + os.path.isfile(os.path.join(ionDir['subband'][1], frameDir, swathDir, self._insar.amplitude)) and \ + os.path.isfile(os.path.join(ionDir['subband'][1], frameDir, swathDir, self._insar.amplitude+'.vrt')) and \ + os.path.isfile(os.path.join(ionDir['subband'][1], frameDir, swathDir, self._insar.amplitude+'.xml')): + print('interferogram already exists at swath {}, frame {}'.format(swathNumber, frameNumber)) + continue + #filter master and slave images for slcx in [self._insar.masterSlc, self._insar.slaveSlc]: slc = os.path.join('../', frameDir, swathDir, slcx) @@ -279,18 +296,56 @@ def runIonSubband(self): #list of input files inputInterferograms = [] inputAmplitudes = [] + phaseDiff = [None] for j, swathNumber in enumerate(range(self._insar.startingSwath, self._insar.endingSwath + 1)): swathDir = 's{}'.format(swathNumber) inputInterferograms.append(os.path.join('../', swathDir, self._insar.interferogram)) inputAmplitudes.append(os.path.join('../', swathDir, self._insar.amplitude)) + #compute phase needed to be compensated using startingRange + if j >= 1: + #phaseDiffSwath1 = -4.0 * np.pi * (masterTrack.frames[i].swaths[j-1].startingRange - slaveTrack.frames[i].swaths[j-1].startingRange)/subbandRadarWavelength[k] + #phaseDiffSwath2 = -4.0 * np.pi * (masterTrack.frames[i].swaths[j].startingRange - slaveTrack.frames[i].swaths[j].startingRange)/subbandRadarWavelength[k] + phaseDiffSwath1 = +4.0 * np.pi * masterTrack.frames[i].swaths[j-1].startingRange * (1.0/radarWavelength - 1.0/subbandRadarWavelength[k]) \ + -4.0 * np.pi * slaveTrack.frames[i].swaths[j-1].startingRange * (1.0/radarWavelength - 1.0/subbandRadarWavelength[k]) + phaseDiffSwath2 = +4.0 * np.pi * masterTrack.frames[i].swaths[j].startingRange * (1.0/radarWavelength - 1.0/subbandRadarWavelength[k]) \ + -4.0 * np.pi * slaveTrack.frames[i].swaths[j].startingRange * (1.0/radarWavelength - 1.0/subbandRadarWavelength[k]) + if masterTrack.frames[i].swaths[j-1].startingRange - slaveTrack.frames[i].swaths[j-1].startingRange == \ + masterTrack.frames[i].swaths[j].startingRange - slaveTrack.frames[i].swaths[j].startingRange: + #phaseDiff.append(phaseDiffSwath2 - phaseDiffSwath1) + #if master and slave versions are all before or after version 2.025 (starting range error < 0.5 m), + #it should be OK to do the above. + #see results in neom where it meets the above requirement, but there is still phase diff + #to be less risky, we do not input values here + phaseDiff.append(None) + else: + phaseDiff.append(None) + #note that frame parameters are updated after mosaicking, here no need to update parameters #mosaic amplitudes swathMosaic(masterTrack.frames[i], inputAmplitudes, self._insar.amplitude, rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, resamplingMethod=0) #mosaic interferograms - swathMosaic(masterTrack.frames[i], inputInterferograms, self._insar.interferogram, - rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateFrame=False, phaseCompensation=True, resamplingMethod=1) + #These are for ALOS-2, may need to change for ALOS-4! + phaseDiffFixed = [0.0, 0.4754024578084084, 0.9509913179406437, 1.4261648478671614, 2.179664007520499, 2.6766909968024932, 3.130810857] + + snapThreshold = 0.2 + (phaseDiffEst, phaseDiffUsed, phaseDiffSource) = swathMosaic(masterTrack.frames[i], inputInterferograms, self._insar.interferogram, + rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateFrame=False, + phaseCompensation=True, phaseDiff=phaseDiff, phaseDiffFixed=phaseDiffFixed, snapThreshold=snapThreshold, pcRangeLooks=1, pcAzimuthLooks=4, + filt=False, resamplingMethod=1) + + #the first item is meaningless for all the following list, so only record the following items + catalog.addItem('{} subswath phase difference input'.format(ionDir['subband'][k]), phaseDiff[1:], 'runIonSubband') + catalog.addItem('{} subswath phase difference estimated'.format(ionDir['subband'][k]), phaseDiffEst[1:], 'runIonSubband') + catalog.addItem('{} subswath phase difference used'.format(ionDir['subband'][k]), phaseDiffUsed[1:], 'runIonSubband') + catalog.addItem('{} subswath phase difference used source'.format(ionDir['subband'][k]), phaseDiffSource[1:], 'runIonSubband') + #check if there is value around 3.130810857, which may not be stable + phaseDiffUnstableExist = False + for xxx in phaseDiffUsed: + if abs(abs(xxx) - 3.130810857) < 0.2: + phaseDiffUnstableExist = True + catalog.addItem('{} subswath phase difference unstable exists'.format(ionDir['subband'][k]), phaseDiffUnstableExist, 'runIonSubband') create_xml(self._insar.amplitude, masterTrack.frames[i].numberOfSamples, masterTrack.frames[i].numberOfLines, 'amp') create_xml(self._insar.interferogram, masterTrack.frames[i].numberOfSamples, masterTrack.frames[i].numberOfLines, 'int') @@ -385,7 +440,8 @@ def runIonSubband(self): os.chdir(ionDir['subband'][k]) for i, frameNumber in enumerate(self._insar.masterFrames): frameDir = 'f{}_{}'.format(i+1, frameNumber) - shutil.rmtree(frameDir) + #keep subswath interferograms + #shutil.rmtree(frameDir) #cmd = 'rm -rf {}'.format(frameDir) #runCmd(cmd) os.chdir('../') diff --git a/components/isceobj/Alos2Proc/runSwathMosaic.py b/components/isceobj/Alos2Proc/runSwathMosaic.py index 9d6b175..3e450e6 100644 --- a/components/isceobj/Alos2Proc/runSwathMosaic.py +++ b/components/isceobj/Alos2Proc/runSwathMosaic.py @@ -162,21 +162,30 @@ def runSwathMosaic(self): self._insar.procDoc.addAllFromCatalog(catalog) -def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks, updateFrame=False, phaseCompensation=False, pcRangeLooks=1, pcAzimuthLooks=4, filt=False, resamplingMethod=0): +def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks, updateFrame=False, phaseCompensation=False, phaseDiff=None, phaseDiffFixed=None, snapThreshold=None, pcRangeLooks=1, pcAzimuthLooks=4, filt=False, resamplingMethod=0): ''' mosaic swaths + #PART 1. REGULAR INPUT PARAMTERS frame: frame inputFiles: input file list - output file: output mosaic file + outputfile: output mosaic file rangeOffsets: range offsets azimuthOffsets: azimuth offsets numberOfRangeLooks: number of range looks of the input files numberOfAzimuthLooks: number of azimuth looks of the input files + updateFrame: whether update frame parameters + + #PART 2. PARAMETERS FOR COMPUTING PHASE DIFFERENCE BETWEEN SUBSWATHS phaseCompensation: whether do phase compensation for each swath + phaseDiff: pre-computed compensation phase for each swath + phaseDiffFixed: if provided, the estimated value will snap to one of these values, which is nearest to the estimated one. + snapThreshold: this is used with phaseDiffFixed pcRangeLooks: number of range looks to take when compute swath phase difference pcAzimuthLooks: number of azimuth looks to take when compute swath phase difference filt: whether do filtering when compute swath phase difference + + #PART 3. RESAMPLING METHOD resamplingMethod: 0: amp resampling. 1: int resampling. ''' from contrib.alos2proc_f.alos2proc_f import rect_with_looks @@ -335,10 +344,33 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num #compute compensation phase for each swath diffMean2 = [0.0 for i in range(numberOfSwaths)] + phaseDiffEst = [None for i in range(numberOfSwaths)] + #True if: + # (1) used diff phase from input + # (2) used estimated diff phase after snapping to a fixed diff phase provided + #False if: + # (1) used purely estimated diff phase + phaseDiffSource = ['estimated' for i in range(numberOfSwaths)] + # 1. 'estimated': estimated from subswath overlap + # 2. 'estimated+snap': estimated from subswath overlap and snap to a fixed value + # 3. 'input': pre-computed + # confidence level: 3 > 2 > 1 if phaseCompensation: #compute swath phase offset diffMean = [0.0] for i in range(1, numberOfSwaths): + + #no need to estimate diff phase if provided from input + ##################################################################### + if phaseDiff!=None: + if phaseDiff[i]!=None: + diffMean.append(phaseDiff[i]) + phaseDiffSource[i] = 'input' + print('using pre-computed phase offset given from input') + print('phase offset: subswath{} - subswath{}: {}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, phaseDiff[i])) + continue + ##################################################################### + #all indexes start with zero, all the computed start/end sample/line indexes are included. #no need to add edge here, as we are going to find first/last nonzero sample/lines later @@ -467,6 +499,19 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num create_xml(filename, width7, length7, 'float') + #save purely estimated diff phase + phaseDiffEst[i] = diffMean0 + + #if fixed diff phase provided and the estimated diff phase is close enough to a fixed value, snap to it + ############################################################################################################ + if phaseDiffFixed != None: + phaseDiffTmp = np.absolute(np.absolute(np.array(phaseDiffFixed)) - np.absolute(diffMean0)) + phaseDiffTmpMinIndex = np.argmin(phaseDiffTmp) + if phaseDiffTmp[phaseDiffTmpMinIndex] < snapThreshold: + diffMean0 = np.sign(diffMean0) * np.absolute(phaseDiffFixed[phaseDiffTmpMinIndex]) + phaseDiffSource[i] = 'estimated+snap' + ############################################################################################################ + diffMean.append(diffMean0) print('phase offset: subswath{} - subswath{}: {}'.format(frame.swaths[i-1].swathNumber, frame.swaths[i].swathNumber, diffMean0)) @@ -503,6 +548,10 @@ def swathMosaic(frame, inputFiles, outputfile, rangeOffsets, azimuthOffsets, num frame.azimuthLineInterval = frame.swaths[0].azimuthLineInterval + if phaseCompensation: + # estimated phase diff, used phase diff, used phase diff source + return (phaseDiffEst, diffMean, phaseDiffSource) + def swathMosaicParameters(frame, rangeOffsets, azimuthOffsets, numberOfRangeLooks, numberOfAzimuthLooks): ''' mosaic swaths (this is simplified version of swathMosaic to update parameters only)