option for discarding burst properties in TOPS ionosphere estimation

applications/topsApp.py

@@ -438,6 +438,20 @@ ION_DO_ION = Application.Parameter('ION_doIon',
     mandatory = False,
     doc = '')
 
+ION_APPLY_ION = Application.Parameter('ION_applyIon',
+    public_name = 'apply ionosphere correction',
+    default = False,
+    type = bool,
+    mandatory = False,
+    doc = '')
+
+ION_CONSIDER_BURST_PROPERTIES = Application.Parameter('ION_considerBurstProperties',
+    public_name = 'consider burst properties in ionosphere computation',
+    default = False,
+    type = bool,
+    mandatory = False,
+    doc = '')
+
 ION_START_STEP = Application.Parameter(
     'ION_startStep',
     public_name='start ionosphere step',
@@ -649,6 +663,8 @@ class TopsInSAR(Application):
                       ########################################################
                       #for ionospheric correction
                       ION_DO_ION,
+                      ION_APPLY_ION,
+                      ION_CONSIDER_BURST_PROPERTIES,
                       ION_START_STEP,
                       ION_END_STEP,
                       ION_ION_HEIGHT,
@@ -734,6 +750,9 @@ class TopsInSAR(Application):
         if(self.geocode_list is None):
             #if not provided by the user use the list from InsarProc
             self.geocode_list = self.insar.geocode_list
+            #for ionosphere
+            if 'topophase.ion' not in self.geocode_list:
+                self.geocode_list.append('topophase.ion')
         else:
             #if geocode_list defined here, then give it to InsarProc
             #for consistency between insarApp and InsarProc and warn the user

components/isceobj/TopsProc/runIon.py

@@ -42,6 +42,7 @@ def setup(self):
     #SECTION 1. PROCESSING CONTROL PARAMETERS
     #1. suggested default values of the parameters
     ionParam.doIon = False
+    ionParam.considerBurstProperties = False
     ionParam.startStep = ionParam.allSteps[0]
     ionParam.endStep = ionParam.allSteps[-1]
 
@@ -77,6 +78,7 @@ def setup(self):
 
     #2. accept the above parameters from topsApp.py
     ionParam.doIon = self.ION_doIon
+    ionParam.considerBurstProperties = self.ION_considerBurstProperties
     ionParam.startStep = self.ION_startStep
     ionParam.endStep = self.ION_endStep
 
@@ -2637,16 +2639,17 @@ def runIon(self):
     if run_step('filt_gaussian', ionParam):
         filt_gaussian(self, ionParam)
 
+    #only do the following steps when considering burst properties
     #ionosphere shift
-    if run_step('ionosphere_shift', ionParam):
+    if run_step('ionosphere_shift', ionParam) and ionParam.considerBurstProperties:
         ionosphere_shift(self, ionParam)
 
     #resample from ionospheric layer to ground layer, get ionosphere for each burst
-    if run_step('ion2grd', ionParam):
+    if run_step('ion2grd', ionParam) and ionParam.considerBurstProperties:
         ion2grd(self, ionParam)
 
     #esd
-    if run_step('esd', ionParam):
+    if run_step('esd', ionParam) and ionParam.considerBurstProperties:
         esd(self, ionParam)
 
     #pure esd without applying ionospheric correction

components/isceobj/TopsProc/runMergeBursts.py

@@ -20,6 +20,38 @@ import logging
 from isceobj.Util.ImageUtil import ImageLib as IML
 
 
+def interpolateDifferentNumberOfLooks(inputData, lengtho, widtho, nrli, nali, nrlo, nalo):
+    '''
+    inputData:   input numpy 2-d array
+    lengtho:     length of output array
+    widtho:      width of output array
+    nrli:        number of range looks input
+    nali:        number of azimuth looks input
+    nrlo:        number of range looks output
+    nalo:        number of azimuth looks output
+    '''
+    import numpy as np
+    from scipy.interpolate import interp1d
+
+    (lengthi, widthi) = inputData.shape
+
+    indexi = np.linspace(0, widthi-1, num=widthi, endpoint=True)
+    indexo = np.linspace(0, widtho-1, num=widtho, endpoint=True) * nrli/nrlo + (nrli-nrlo)/(2.0*nrlo)
+    outputData0 = np.zeros((lengthi, widtho), dtype=inputData.dtype)
+    for i in range(lengthi):
+        f = interp1d(indexi, inputData[i,:], kind='cubic', fill_value="extrapolate")
+        outputData0[i, :] = f(indexo)
+    
+    indexi = np.linspace(0, lengthi-1, num=lengthi, endpoint=True)
+    indexo = np.linspace(0, lengtho-1, num=lengtho, endpoint=True) * nali/nalo + (nali-nalo)/(2.0*nalo)
+    outputData = np.zeros((lengtho, widtho), dtype=inputData.dtype)
+    for j in range(widtho):
+        f = interp1d(indexi, outputData0[:, j], kind='cubic', fill_value="extrapolate")
+        outputData[:, j] = f(indexo)
+
+    return outputData
+
+
 def mergeBox(frame):
     '''
     Merging using VRTs.
@@ -666,14 +698,16 @@ def runMergeBursts(self, adjust=1):
     #totalLooksThreshold = 9
     totalLooksThreshold = 99999999999999
     #if doing ionospheric correction
-    ionCorrection = self.ION_doIon
+    doIon = self.ION_doIon
+    applyIon = self.ION_applyIon
+    considerBurstProperties = self.ION_considerBurstProperties
     ionDirname = 'ion/ion_burst'
     mergedIonname = 'topophase.ion'
     originalIfgname = 'topophase_ori.flat'
     #########################################
 
     # backing out the tigher constraints for ionosphere as it could itnroduce gabs between along track products produced seperately  
-    if not ionCorrection:
+    if not (doIon and considerBurstProperties):
         adjust=0
 
     #########################################
@@ -712,7 +746,7 @@ def runMergeBursts(self, adjust=1):
             #restore
             frames = frames_bak
         else:
-            validOnly==True
+            validOnly=True
 
 
     #########################################
@@ -738,7 +772,7 @@ def runMergeBursts(self, adjust=1):
         mergeBursts2(frames, os.path.join(self._insar.fineIfgDirname, 'IW%d',  'burst_%02d.int'), burstIndex, box, os.path.join(mergedir, self._insar.mergedIfgname+suffix), virtual=virtual, validOnly=True)
         if self.numberAzimuthLooks * self.numberRangeLooks < totalLooksThreshold:
             mergeBursts2(frames, os.path.join(self._insar.fineIfgDirname, 'IW%d',  'burst_%02d.cor'), burstIndex, box, os.path.join(mergedir, self._insar.correlationFilename+suffix), virtual=virtual, validOnly=True)
-        if ionCorrection == True:
+        if doIon and considerBurstProperties:
             mergeBursts2(frames, os.path.join(ionDirname, 'IW%d',  'burst_%02d.ion'), burstIndex, box, os.path.join(mergedir, mergedIonname+suffix), virtual=virtual, validOnly=True)
 
 
@@ -782,13 +816,61 @@ def runMergeBursts(self, adjust=1):
                 os.remove(os.path.join(mergedir, pwrfile+'.xml'))
                 os.remove(os.path.join(mergedir, pwrfile+'.vrt'))
 
-            if ionCorrection:
+            if doIon:
-                multilook(os.path.join(mergedir, mergedIonname+suffix),
+                if considerBurstProperties:
-                  outname = os.path.join(mergedir, mergedIonname),
+                    multilook(os.path.join(mergedir, mergedIonname+suffix),
-                  alks = self.numberAzimuthLooks, rlks=self.numberRangeLooks)
+                      outname = os.path.join(mergedir, mergedIonname),
+                      alks = self.numberAzimuthLooks, rlks=self.numberRangeLooks)
+                else:
+                    ionFilt = 'ion/ion_cal/filt.ion'
+                    img = isceobj.createImage()
+                    img.load(ionFilt+'.xml')
+                    ionFiltImage = (np.fromfile(ionFilt, dtype=np.float32).reshape(img.length*2, img.width))[1:img.length*2:2, :]
+                    img = isceobj.createImage()
+                    img.load(os.path.join(mergedir, self._insar.mergedIfgname)+'.xml')
+
+                    #interpolate original
+                    ionFiltImageOut = interpolateDifferentNumberOfLooks(ionFiltImage, img.length, img.width, self.numberRangeLooks, self.numberAzimuthLooks, self.ION_numberRangeLooks, self.ION_numberAzimuthLooks)
+                    ionFiltOut = os.path.join(mergedir, mergedIonname)
+                    ionFiltImageOut.astype(np.float32).tofile(ionFiltOut)
+
+                    image = isceobj.createImage()
+                    image.setDataType('FLOAT')
+                    image.setFilename(ionFiltOut)
+                    image.extraFilename = ionFiltOut + '.vrt'
+                    image.setWidth(img.width)
+                    image.setLength(img.length)
+                    #image.setAccessMode('read')
+                    #image.createImage()
+                    image.renderHdr()
+                    #image.finalizeImage()
     else:
         print('Skipping multi-looking ....')
 
+        if self.doInSAR and doIon and (not considerBurstProperties):
+            ionFilt = 'ion/ion_cal/filt.ion'
+            img = isceobj.createImage()
+            img.load(ionFilt+'.xml')
+            ionFiltImage = (np.fromfile(ionFilt, dtype=np.float32).reshape(img.length*2, img.width))[1:img.length*2:2, :]
+            img = isceobj.createImage()
+            img.load(os.path.join(mergedir, self._insar.mergedIfgname+suffix)+'.xml')
+
+            #interpolate original
+            ionFiltImageOut = interpolateDifferentNumberOfLooks(ionFiltImage, img.length, img.width, self.numberRangeLooks, self.numberAzimuthLooks, self.ION_numberRangeLooks, self.ION_numberAzimuthLooks)
+            ionFiltOut = os.path.join(mergedir, mergedIonname)
+            ionFiltImageOut.astype(np.float32).tofile(ionFiltOut)
+
+            image = isceobj.createImage()
+            image.setDataType('FLOAT')
+            image.setFilename(ionFiltOut)
+            image.extraFilename = ionFiltOut + '.vrt'
+            image.setWidth(img.width)
+            image.setLength(img.length)
+            #image.setAccessMode('read')
+            #image.createImage()
+            image.renderHdr()
+            #image.finalizeImage()
+
 
     #########################################
     # STEP 4. APPLY CORRECTIONS
@@ -796,8 +878,8 @@ def runMergeBursts(self, adjust=1):
     #do ionospheric and other corrections here
     #should also consider suffix, but usually we use multiple looks, so I ignore it for now.
     if self.doInSAR:
-        if ionCorrection:
+        if doIon and applyIon:
-            print('user choose to do ionospheric correction')
+            print('user choose to apply ionospheric correction')
 
             #define file names
             interferogramFilename = os.path.join(mergedir, self._insar.mergedIfgname)

examples/input_files/topsApp.xml

@@ -183,6 +183,8 @@ By default, ESD is always performed.
 #ionospheric correction module parameters
 #the values below are the default values used by the module
 <property name="do ionosphere correction">False</property>
+<property name="apply ionosphere correction">False</property>
+<property name="consider burst properties in ionosphere computation">False</property>
 #choose from: ['subband', 'rawion', 'grd2ion', 'filt_gaussian', 'ionosphere_shift', 'ion2grd', 'esd']
 <property name="start ionosphere step">subband</property>
 <property name="end ionosphere step">esd</property>