removing azimuth burst ramps caused by ionosphere in topsStack (#600)

* removing azimuth burst ramps caused by ionosphere in topsStack These updates are for removing azimuth burst ramps caused by ionosphere in the ionospheric correction in topsStack. The updates are only within this module. * Update Stack.py * update argument for s1_select_ion.py * Update runFrameOffset.py * Update contrib/stack/topsStack/README.md Co-authored-by: Zhang Yunjun <[email protected]> * Update contrib/stack/topsStack/stackSentinel.py Co-authored-by: Zhang Yunjun <[email protected]> * Update contrib/stack/topsStack/stackSentinel.py Co-authored-by: Zhang Yunjun <[email protected]> * Update contrib/stack/topsStack/stackSentinel.py Co-authored-by: Zhang Yunjun <[email protected]> * Update contrib/stack/topsStack/stackSentinel.py Co-authored-by: Zhang Yunjun <[email protected]> --------- Co-authored-by: Zhang Yunjun <[email protected]>
isce-framework · Jul 1, 2024 · b231d6e · b231d6e
1 parent ca7649a
commit b231d6e
Show file tree

Hide file tree

Showing 9 changed files with 631 additions and 17 deletions.
diff --git a/components/isceobj/Alos2Proc/runFrameOffset.py b/components/isceobj/Alos2Proc/runFrameOffset.py
@@ -97,18 +97,23 @@ def frameOffset(track, image, outputfile, crossCorrelation=True, matchingMode=0)
         swath1 = track.frames[j-1].swaths[0]
         swath2 = track.frames[j].swaths[0]
 
+        #consider frame/swath azimuth sensing start differences caused by swath mosaicking
+        # tested with alos2App.py, alos2burstApp.py and alosStack
+        frame1 = track.frames[j-1]
+        frame2 = track.frames[j]
+        delta_az = -((frame2.sensingStart - frame1.sensingStart).total_seconds() - (swath2.sensingStart - swath1.sensingStart).total_seconds()) / swath1.azimuthLineInterval
 
         #offset from geometry
         offsetGeometrical = computeFrameOffset(swath1, swath2)
         rangeOffsetGeometrical.append(offsetGeometrical[0])
-        azimuthOffsetGeometrical.append(offsetGeometrical[1])
+        azimuthOffsetGeometrical.append(offsetGeometrical[1]+delta_az)
 
         #offset from cross-correlation
         if crossCorrelation:
             offsetMatching = estimateFrameOffset(swath1, swath2, image1, image2, matchingMode=matchingMode)
             if offsetMatching != None:
                 rangeOffsetMatching.append(offsetMatching[0])
-                azimuthOffsetMatching.append(offsetMatching[1])
+                azimuthOffsetMatching.append(offsetMatching[1]+delta_az)
             else:
                 print('******************************************************************')
                 print('WARNING: bad matching offset, we are forced to use')

diff --git a/contrib/stack/topsStack/README.md b/contrib/stack/topsStack/README.md
@@ -186,7 +186,7 @@ Ionospheric phase estimation has more requirements than regular InSAR processing
 In stack ionospheric phase estimation, acquistions with same swath starting ranges are put in a group. A network is formed within a group. Extra pairs are also processed to connect the different groups so that all acquistions are connected. But we need to estimate a phase offset for these extra pairs, which might not be accurate. Therefore, for a particualr swath starting ranges, if there are only a few acquistions, it's better to just discard them so that we don't have to estimate the phase offsets.
 
 ```
-s1_select_ion.py -dir data/slc -sn 33.550217/37.119545 -nr 10
+s1_select_ion.py -dir data/slc -sn 33.550217 37.119545 -nr 10
 ```
 
 Acquistions to be used need to fully cover the south/north bounds. After running this command, acquistion not to be used will be put in a folder named 'not_used'. It's OK to run this command multiple times.
@@ -204,6 +204,8 @@ An example --param_ion file 'ion_param.txt' is provided in the code directory. F
 stackSentinel.py -s ../data/slc -d ../data/dem/dem_1_arcsec/demLat_N32_N41_Lon_W113_W107.dem.wgs84 -b '33.550217 37.119545 -111.233932 -107.790451' -a ../data/s1_aux_cal -o ../data/orbit -C geometry -c 2 --param_ion ../code/ion_param.txt --num_connections_ion 3
 ```
 
+Note in 'ion_param.txt', if 'consider burst properties in ionosphere computation' is True, Coregistration options '-C', '--coregistration' in stackSentinel.py must be NESD.
+
 If ionospheric phase estimation is enabled in stackSentinel.py, it will generate the following run files. Here ***ns*** means number of steps in the original stack processing, which depends on the type of stack (slc, correlation, interferogram, and offset).
 
 -	run_ns+1_subband_and_resamp
@@ -251,6 +253,22 @@ Estimate ionospheric phase for each date. We highly recommend inspecting all pai
 
 Typical anomalies include dense fringes caused by phase unwrapping errors, and a range ramp as a result of errors in estimating phase offsets for pairs with different swath starting ranges (check pairs_diff_starting_ranges.txt).
 
+**run_ns+9_filtIonShift**
+
+Filter azimuth ionospheric shift.
+
+**run_ns+10_invertIonShift**
+
+Estimate azimuth ionospheric shift for each date. As in step **run_ns+8_invertIon**, check if there are anamolies.
+
+**run_ns+11_burstRampIon**
+
+Compute azimuth burst ramps as a result of ionosphere for each date.
+
+**run_ns+12_mergeBurstRampIon**
+
+Merge azimuth burst ramps.
+
 #### 3. run command files generated ####
 
 Run the commands sequentially.
@@ -261,7 +279,11 @@ Results from ionospheric phase estimation.
 
 -	reference and coreg_secondarys: now contains also subband burst SLCs
 -	ion: original ionospheric phase estimation results
+-       ion_azshift_dates: azimuth ionospheric shift for each acquistion
+-       ion_burst_ramp_dates: azimuth burst ramps caused by ionosphere for each acquistion
+-       ion_burst_ramp_merged_dates: merged azimuth burst ramps caused by ionosphere for each acquistion
 -	ion_dates: ionospheric phase for each acquistion
+-	ion/date1_date2/ion_cal/azshift.ion: azimuth ionospheric shift
 -	ion/date1_date2/ion_cal/filt.ion: filtered ionospheric phase
 -	ion/date1_date2/ion_cal/raw_no_projection.ion: original ionospheric phase
 -	ion/date1_date2/lower/merged/fine_look.unw: unwrapped lower band interferogram
@@ -273,6 +295,7 @@ If ionospheric phase estimation processing is swath by swath because of differen
 -	ion/date1_date2/lower/merged_IW*
 -	ion/date1_date2/upper/merged_IW*
 
+Unit of azimuth ionospheric shift is number of single look azimuth lines.
 After processing, we can plot ionospheric phase estimation results using plotIonPairs.py and plotIonDates.py. For example
 
 ```
@@ -285,4 +308,5 @@ Relationships of the ionospheric phases:
 ```
 ion_dates/date1.ion - ion_dates/date2.ion = ion/date1_date2/ion_cal/filt.ion
 ion_dates/date1.ion - ion_dates/date2.ion = ionospheric phase in merged/interferograms/date1_date2/filt_fine.unw
+ion_burst_ramp_merged_dates/date1.float - ion_burst_ramp_merged_dates/date2.float = azimuth burst ramps caused by ionosphere in merged/interferograms/date1_date2/filt_fine.unw
 ```
diff --git a/contrib/stack/topsStack/Stack.py b/contrib/stack/topsStack/Stack.py
@@ -339,8 +339,33 @@ def filtIon(self, function):
         self.f.write('win_min : ' + '{}'.format(self.win_min) + '\n')
         self.f.write('win_max : ' + '{}'.format(self.win_max) + '\n')
 
+    def filtIonShift(self, function):
+        self.f.write('###################################'+'\n')
+        self.f.write(function + '\n')
+        self.f.write('filtIonShift : ' + '\n')
+        self.f.write('reference_stack : ' + self.reference_stack + '\n')
+        self.f.write('reference : ' + self.reference + '\n')
+        self.f.write('secondary : ' + self.secondary + '\n')
+        self.f.write('input : ' + self.input + '\n')
+        self.f.write('coherence : ' + self.coherence + '\n')
+        self.f.write('output : ' + self.output + '\n')
+        self.f.write('win_min : ' + '{}'.format(self.win_min) + '\n')
+        self.f.write('win_max : ' + '{}'.format(self.win_max) + '\n')
+        self.f.write('nrlks : ' + '{}'.format(self.nrlks) + '\n')
+        self.f.write('nalks : ' + '{}'.format(self.nalks) + '\n')
 
 
+    def burstRampIon(self, function):
+        self.f.write('###################################'+'\n')
+        self.f.write(function + '\n')
+        self.f.write('burstRampIon : ' + '\n')
+        self.f.write('reference_stack : ' + self.reference_stack + '\n')
+        self.f.write('input : ' + self.input + '\n')
+        self.f.write('output : ' + self.output + '\n')
+        self.f.write('nrlks : ' + '{}'.format(self.nrlks) + '\n')
+        self.f.write('nalks : ' + '{}'.format(self.nalks) + '\n')
+        self.f.write('ion_height : ' + '{}'.format(self.ion_height) + '\n')
+
 
     def write_wrapper_config2run_file(self, configName, line_cnt, numProcess = 1):
         # dispassionate list of commands for single process
@@ -1330,7 +1355,7 @@ def unwrap_ion(self, pairs_same_starting_ranges_update, pairs_diff_starting_rang
                     configObj.defoMax = '2'
                     configObj.rangeLooks = '{}'.format(ionParamUsrObj.ION_numberRangeLooks0)
                     configObj.azimuthLooks = '{}'.format(ionParamUsrObj.ION_numberAzimuthLooks0)
-                    configObj.rmfilter = False
+                    configObj.rmFilter = False
                     configObj.unwMethod = 'snaphu'
                     configObj.unwrap(function)
 
@@ -1348,7 +1373,7 @@ def unwrap_ion(self, pairs_same_starting_ranges_update, pairs_diff_starting_rang
                     configObj.defoMax = '2'
                     configObj.rangeLooks = '{}'.format(ionParamUsrObj.ION_numberRangeLooks0)
                     configObj.azimuthLooks = '{}'.format(ionParamUsrObj.ION_numberAzimuthLooks0)
-                    configObj.rmfilter = False
+                    configObj.rmFilter = False
                     configObj.unwMethod = 'snaphu'
                     configObj.unwrap(function)
 
@@ -1542,6 +1567,105 @@ def invertIon(self):
         self.runf.write(self.text_cmd + cmd + '\n')
 
 
+    def filtIonShift(self, pairs):
+
+        ionParamUsrObj = ionParamUsr(self.param_ion)
+        ionParamUsrObj.configure()
+
+        line_cnt = 0
+        for p in pairs:
+            configName = os.path.join(self.config_path,'config_filtIonShift_{}_{}'.format(p[0], p[1]))
+            configObj = config(configName)
+            configObj.configure(self)
+            configObj.reference_stack = os.path.join(self.work_dir, 'reference')
+            configObj.reference = os.path.join(self.work_dir, 'coreg_secondarys', '{}'.format(p[0]))
+            configObj.secondary = os.path.join(self.work_dir, 'coreg_secondarys', '{}'.format(p[1]))
+            configObj.input = os.path.join(self.work_dir, 'ion', '{}_{}'.format(p[0], p[1]), 'ion_cal', 'filt.ion')
+            configObj.coherence = os.path.join(self.work_dir, 'ion', '{}_{}'.format(p[0], p[1]), 'ion_cal', 'raw_no_projection.cor')
+            configObj.output = os.path.join(self.work_dir, 'ion', '{}_{}'.format(p[0], p[1]), 'ion_cal', 'azshift.ion')
+            configObj.win_min = ionParamUsrObj.ION_ionshiftFilteringWinsizeMin
+            configObj.win_max = ionParamUsrObj.ION_ionshiftFilteringWinsizeMax
+            configObj.nrlks = ionParamUsrObj.ION_numberRangeLooks
+            configObj.nalks = ionParamUsrObj.ION_numberAzimuthLooks
+            configObj.filtIonShift('[Function-1]')
+            configObj.finalize()
+
+            line_cnt += 1
+            #line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt, self.numProcess)
+            line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt)
+            del configObj
+
+
+    def invertIonShift(self):
+
+        ionParamUsrObj = ionParamUsr(self.param_ion)
+        ionParamUsrObj.configure()
+
+        ion_in = os.path.join(self.work_dir,'ion')
+        ion_out = os.path.join(self.work_dir,'ion_azshift_dates')
+        hgt = os.path.join(self.work_dir,'merged/geom_reference/hgt.rdr')
+
+        cmd = 'invertIon.py --idir {} --filename azshift.ion --odir {} --nrlks1 {} --nalks1 {} --nrlks2 {} --nalks2 {} --merged_geom {} --interp --msk_overlap'.format(ion_in,ion_out,ionParamUsrObj.ION_numberRangeLooks, ionParamUsrObj.ION_numberAzimuthLooks, self.rangeLooks, self.azimuthLooks,hgt)
+
+        self.runf.write(self.text_cmd + cmd + '\n')
+
+
+    def burstRampIon(self, dates):
+
+        ionParamUsrObj = ionParamUsr(self.param_ion)
+        ionParamUsrObj.configure()
+
+        line_cnt = 0
+        for p in dates:
+            configName = os.path.join(self.config_path,'config_burstRampIon_{}'.format(p))
+            configObj = config(configName)
+            configObj.configure(self)
+            configObj.reference_stack = os.path.join(self.work_dir, 'reference')
+            configObj.input = os.path.join(self.work_dir, 'ion_azshift_dates', '{}.ion'.format(p))
+            configObj.output = os.path.join(self.work_dir, 'ion_burst_ramp_dates', '{}'.format(p))
+            configObj.nrlks = self.rangeLooks
+            configObj.nalks = self.azimuthLooks
+            configObj.ion_height = ionParamUsrObj.ION_ionHeight
+            configObj.burstRampIon('[Function-1]')
+            configObj.finalize()
+
+            line_cnt += 1
+            #line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt, self.numProcess)
+            line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt)
+            del configObj
+
+
+    def mergeBurstRampIon(self, dates, stackReferenceDate):
+
+        line_cnt = 0
+        for d in dates:
+            configName = os.path.join(self.config_path,'config_mergeBurstRampIon_' + d)
+            configObj = config(configName)
+            configObj.configure(self)
+            configObj.stack = os.path.join(self.work_dir, 'stack')
+            if d == stackReferenceDate:
+                configObj.reference = os.path.join(self.work_dir, 'reference')
+            else:
+                configObj.reference = os.path.join(self.work_dir, 'coreg_secondarys/' + d)
+            configObj.dirName = os.path.join(self.work_dir, 'ion_burst_ramp_dates/' + d)
+            configObj.namePattern = 'burst*float'
+            configObj.mergedFile = os.path.join(self.work_dir, 'ion_burst_ramp_merged_dates/' + d + '.float')
+            configObj.mergeBurstsMethod = 'top'
+            if d == stackReferenceDate:
+                configObj.aligned = 'False'
+            else:
+                configObj.aligned = 'True'
+            configObj.validOnly = 'True'
+            configObj.useVirtualFiles = 'True'
+            configObj.multiLook = 'True'
+            configObj.mergeBurst('[Function-1]')
+            configObj.finalize()
+
+            line_cnt += 1
+            line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt)
+            del configObj
+
+
     def finalize(self):
         self.runf.close()
         #writeJobFile(self.run_outname)