ionospheric phase estimation module for topsStack (#326)

* ionosphere correction for topsStack

* Update README.md

* Update Alos2ProcPublic.py

* Update Stack.py

Currently the script writes the `invertIon.py` command pointing to 'ion' and writing to 'ion_dates'. However we need the absolute paths in order to find the 'ion' directory from within the run_files directory (similar to how the `computeIon.py` function is written at l.1488)

* Add files via upload

* Add files via upload

* Add files via upload

* Add files via upload

* Add files via upload

Co-authored-by: Oliver Stephenson <[email protected]>

Author: CunrenLiang

components/isceobj/Alos2Proc/Alos2ProcPublic.py

@@ -222,8 +222,8 @@ def readOffset(filename):
     offsets = OffsetField()
     for linex in lines:
         #linexl = re.split('\s+', linex)
-        #detect blank lines with only spaces and tabs
-        if linex.strip() == '':
+        #detect blank lines with only spaces and tabs, lines with invalid numbers
+        if (linex.strip() == '') or ('*' in linex):
             continue
 
         linexl = linex.split()

components/isceobj/Scene/Track.py

@@ -237,7 +237,8 @@ def reAdjustStartLine(self, sortedList, width):
             indx = self.findOverlapLine(sortedList[i-1][1],sortedList[i][1],endLine,width,i-1,i)
             #if indx is not None than indx is the new start line
             #otherwise we use startLine  computed from acquisition time
-            if (indx is not None) and (indx + sortedList[i-1][0] != sortedList[i][0]):
+            #no need to do this for ALOS; otherwise there will be problems when there are multiple prfs and the data are interpolated. C. Liang, 20-dec-2021
+            if (self._frames[0].instrument.platform._mission != 'ALOS') and (indx is not None) and (indx + sortedList[i-1][0] != sortedList[i][0]):
                 startLine.append(indx + sortedList[i-1][0])
                 outputs.append(sortedList[i][1])
                 self.logger.info("Changing starting line for frame %d from %d to %d"%(i,endLine,indx))

components/isceobj/Sensor/ALOS.py

@@ -185,8 +185,10 @@ def populateMetadata(self):
 
         productLevel = float(self.leaderFile.sceneHeaderRecord.metadata[
             'Product level code'])
-        if productLevel == 1.0:
-            self.updateRawParameters()
+
+        #this proves to be not accurate. Adjusting first frame is OK, but adjusting following frames would cause discontinuities between frames. C. Liang, 20-dec-2021
+        # if productLevel == 1.0:
+        #     self.updateRawParameters()
         pass
 
     def _populatePlatform(self):

components/isceobj/Sensor/src/ALOS_pre_process/ALOS_pre_process.c

@@ -409,7 +409,7 @@ char   	date[8];
 	    //line number of first line of this prf file
 	    line_number_first = (rspi_new.SC_clock_start[i] - rspi_pre[0].SC_clock_start[0]) * d2s / (1.0 / rspi_pre[0].prf[0]);
 	    //unit: pri of first prf of first image
-	    num_lines_out = (rspi_pre[image_i].frame_counter_end[i] - rspi_pre[image_i].frame_counter_start[i] + 1) * (1.0/rspi_pre[image_i].prf[i]) / (1.0/rspi_pre[0].prf[0]);
+	    num_lines_out = (int)((rspi_pre[image_i].frame_counter_end[i] - rspi_pre[image_i].frame_counter_start[i] + 1) * (1.0/rspi_pre[image_i].prf[i]) / (1.0/rspi_pre[0].prf[0]));
 
         if((fabs(roundfi(line_number_first)-line_number_first)<0.1) && (rspi_pre[image_i].prf[i]==rspi_pre[0].prf[0]))
         	continue;
@@ -469,7 +469,7 @@ char   	date[8];
 	    //append prf i
 	    for(i = 1; i < rspi_new.nPRF; i++){
             //number of lines to be appended between frames if there are gaps
-            num_lines_append = (rspi_new.SC_clock_start[i] - rspi_new.SC_clock_start[0]) * d2s / (1.0/rspi_pre[0].prf[0]) - rspi_new.num_lines[0];
+            num_lines_append = roundfi((rspi_new.SC_clock_start[i] - rspi_new.SC_clock_start[0]) * d2s / (1.0/rspi_pre[0].prf[0])) - rspi_new.num_lines[0];
             if(num_lines_append >= 1){
             	for(j = 0; j < num_lines_append; j++){
             		for(k = 0; k < 2*rspi_new.num_bins[i]; k++)
@@ -485,7 +485,7 @@ char   	date[8];
 	        	die("can't open", rspi_new.input_file[i]);
             num_lines_append = 0;
             for(j = 0; j < rspi_new.num_lines[i]; j++){
-                if((rspi_new.SC_clock_start[i] + j * (1.0/rspi_pre[0].prf[0]) / d2s -  rspi_new.SC_clock_start[0]) * d2s / (1.0/rspi_pre[0].prf[0]) >= rspi_new.num_lines[0]){
+                if(roundfi((rspi_new.SC_clock_start[i] + j * (1.0/rspi_pre[0].prf[0]) / d2s -  rspi_new.SC_clock_start[0]) * d2s / (1.0/rspi_pre[0].prf[0])) >= rspi_new.num_lines[0]){
 			        if(fread((char *)data, 2*sizeof(char)*rspi_new.num_bins[i], 1, next_prf_fp) != 1)
 			        	die("can't read data from", rspi_new.input_file[i]);
 			        if(fwrite((char *)data, 2*sizeof(char)*rspi_new.num_bins[i], 1, first_prf_fp) != 1)

components/isceobj/TopsProc/runIon.py

@@ -64,6 +64,7 @@ def setup(self):
     #1: use mean value of a burst
     #2: use full burst
     ionParam.azshiftFlag = 1
+    ionParam.maskedAreas = None
 
     #better NOT try changing the following two parameters, since they are related
     #to the filtering parameters above
@@ -1043,8 +1044,8 @@ def multilook_unw(self, ionParam, mergedDirname):
     lowerint = np.fromfile(lowerIntName, dtype=np.complex64).reshape(lengthNew, widthNew)
     upperint = np.fromfile(upperIntName, dtype=np.complex64).reshape(lengthNew, widthNew)
     cor = np.zeros((lengthNew*2, widthNew), dtype=np.float32)
-    cor[0:length*2:2, :] = np.sqrt( (np.absolute(lowerint)+np.absolute(upperint))/2.0 )
-    cor[1:length*2:2, :] = cal_coherence(lowerint*np.conjugate(upperint), win=3, edge=4)
+    cor[0:lengthNew*2:2, :] = np.sqrt( (np.absolute(lowerint)+np.absolute(upperint))/2.0 )
+    cor[1:lengthNew*2:2, :] = cal_coherence(lowerint*np.conjugate(upperint), win=3, edge=4)
     cor.astype(np.float32).tofile(corName)
 
 

contrib/stack/topsStack/README.md

@@ -165,3 +165,124 @@ stackSentinel.py -s ../SLC/ -d ../../MexicoCity/demLat_N18_N20_Lon_W100_W097.dem
 This workflow is basically similar to the previous one. The difference is that the interferograms are not unwrapped.
 
 #### 5. Execute the commands in run files (run_01*, run_02*, etc) in the "run_files" folder ####
+
+
+
+-----------------------------------
+### Ionospheric Phase Estimation
+
+Ionospheric phase estimation can be performed for each of the workflow introduced above. Generally, we should do ionospheric phase estimation for data acquired at low latitudes on ascending tracks. However, ionospheric phase estimation only works well in areas with high coherence since it requires phase unwrapping. Details about Sentinel-1 ionospheric correction can be found in
+
++ C. Liang, P. Agram, M. Simons, and E. J. Fielding, “Ionospheric correction of InSAR time series analysis of C-band Sentinel-1 TOPS data,” IEEE Transactions on Geoscience and Remote Sensing, vol. 57, no. 9, pp. 6755-6773, Sep. 2019.
+
+Ionospheric phase estimation has more requirements than regular InSAR processing. The most important two requirements include
+
+-	All the acquistions need to be connected in the entire network.
+
+-	The swath starting ranges need to be the same among all acquistions; otherwise a phase offset between adjacent swaths need to be estimated in order to maintain consistency among the swaths, which might not be accurate enough.
+
+#### 1. select the usable acquistions ####
+
+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
+```
+
+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.
+
+#### 2. generate configure and run files ####
+
+In stackSentinel.py, two options are for ionospheric phase estimation
+
+-	--param_ion
+-	--num_connections_ion
+
+An example --param_ion file 'ion_param.txt' is provided in the code directory. For example, we want to do ionospheric phase estimation when processing a stack of interferograms
+
+```
+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
+```
+
+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
+-	run_ns+2_generateIgram_ion
+-	run_ns+3_mergeBurstsIon
+-	run_ns+4_unwrap_ion
+-	run_ns+5_look_ion
+-	run_ns+6_computeIon
+-	run_ns+7_filtIon
+-	run_ns+8_invertIon
+
+Note about **'areas masked out in ionospheric phase estimation'** in ion_param.txt. Seperated islands or areas usually lead to phase unwrapping errors and therefore significantly affect ionospheric phase estimation. It's better to mask them out. Check ion/date1_date2/ion_cal/raw_no_projection.ion for areas to be masked out. However, we don't have this file before processing the data. To quickly get this file, we can process a stack of two acquistions to get this file. NOTE that the reference of this two-acquisition stack should be the same as that of the full stack we want to process.
+
+**run_ns+1_subband_and_resamp**
+
+Two subband burst SLCs are generated for each burst. For secondary acquistions, the subband burst SLCs are also resampled to match reference burst SLCs. If the subband burst SLCs already exists, the program simply skips the burst.
+
+**run_ns+2_generateIgram_ion**
+
+Generate subband burst interferograms.
+
+**run_ns+3_mergeBurstsIon**
+
+Merge subband burst interferograms, and create a unique coherence file used in ionospheric phase estimation. This will be done swath by swath if the two acquistions of the pair has different swath starting ranges.
+
+**run_ns+4_unwrap_ion**
+
+Unwrap merged subband interferograms. This will be done swath by swath if the two acquistions of the pair has different swath starting ranges.
+
+**run_ns+5_look_ion**
+
+Take further looks on the unwrapped interferograms, and create the unique coherence file based on the further number of looks. This will be done swath by swath if the two acquistions of the pair has different swath starting ranges.
+
+**run_ns+6_computeIon**
+
+Compute ionospheric phase. This will be done swath by swath if the two acquistions of the pair has different swath starting ranges, and then the swath ionospheric phase will be merged.
+
+**run_ns+7_filtIon**
+
+Filter ionospheric phase.
+
+**run_ns+8_invertIon**
+
+Estimate ionospheric phase for each date. We highly recommend inspecting all pair ionospheric phases ion/date1_date2/ion_cal/filt.ion, and exclude those with anomalies in this command.
+
+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).
+
+#### 3. run command files generated ####
+
+Run the commands sequentially.
+
+#### 4. check results ####
+
+Results from ionospheric phase estimation.
+
+-	reference and coreg_secondarys: now contains also subband burst SLCs
+-	ion: original ionospheric phase estimation results
+-	ion_dates: ionospheric phase for each acquistion
+-	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
+-	ion/date1_date2/upper/merged/fine_look.unw: unwrapped upper band interferogram
+
+If ionospheric phase estimation processing is swath by swath because of different swath starting ranges, there will be swath processing directories including
+
+-	ion/date1_date2/ion_cal_IW*
+-	ion/date1_date2/lower/merged_IW*
+-	ion/date1_date2/upper/merged_IW*
+
+After processing, we can plot ionospheric phase estimation results using plotIonPairs.py and plotIonDates.py. For example
+
+```
+plotIonPairs.py -idir ion -odir ion_plot
+plotIonDates.py -idir ion_dates -odir ion_dates_plot
+```
+
+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
+```