improve ALOS-2 ionospheric correciton and geocode

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

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, 
+            #These are for ALOS-2, may need to change for ALOS-4!
-                rangeOffsets, azimuthOffsets, self._insar.numberRangeLooks1, self._insar.numberAzimuthLooks1, updateFrame=False, phaseCompensation=True, resamplingMethod=1)
+            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('../')

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)