removing azimuth burst ramps caused by ionosphere in topsStack

contrib/stack/topsStack/README.md

@@ -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
 ```

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
@@ -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)

contrib/stack/topsStack/burstRampIon.py

@@ -0,0 +1,169 @@
+#!/usr/bin/env python3
+
+# Author: Cunren Liang
+# Copyright 2022
+
+import os
+import copy
+import glob
+import shutil
+import argparse
+import numpy as np
+
+import isce
+import isceobj
+import s1a_isce_utils as ut
+from VRTManager import Swath
+from isceobj.TopsProc.runIon import computeDopplerOffset
+
+
+def createParser():
+    parser = argparse.ArgumentParser(description='compute burst phase ramps using azimuth ionospheric shift')
+    parser.add_argument('-k', '--reference_stack', type=str, dest='reference_stack', required=True,
+            help='Directory with the reference image of the stack')
+    parser.add_argument('-i', '--input', dest='input', type=str, required=True,
+            help='input azimuth ionospheric shift')
+    parser.add_argument('-o', '--output', dest='output', type=str, required=True,
+            help='Directory with output')
+    parser.add_argument('-r', '--nrlks', dest='nrlks', type=int, default=1, 
+            help='number of range looks of azimuth ionospheric shift. Default: 1')
+    parser.add_argument('-a', '--nalks', dest='nalks', type=int, default=1, 
+            help='number of azimuth looks of azimuth ionospheric shift. Default: 1')
+    parser.add_argument('-t', '--ion_height', dest='ion_height', type=float, default=200.0, 
+            help='height of ionospheric layer above the Earth surface in km. Default: 200.0')
+
+    return parser
+
+
+def cmdLineParse(iargs = None):
+    parser = createParser()
+    return parser.parse_args(args=iargs)
+
+
+def main(iargs=None):
+    '''
+    compute burst phase ramps using azimuth ionospheric shift
+    both regular and subband bursts are merged in a box defined by reference of stack
+
+    '''
+    inps = cmdLineParse(iargs)
+
+    #ionospheric layer height (m)
+    ionHeight = inps.ion_height * 1000.0
+    #earth's radius (m)
+    earthRadius = 6371 * 1000.0
+
+    img = isceobj.createImage()
+    img.load(inps.input + '.xml')
+    width = img.width
+    length = img.length
+    ionShift = np.fromfile(inps.input, dtype=np.float32).reshape(length, width)
+
+    swathList = ut.getSwathList(inps.reference_stack)
+    frameReferenceAll = [ut.loadProduct(os.path.join(inps.reference_stack, 'IW{0}.xml'.format(swath))) for swath in swathList]
+
+    refSwaths = [Swath(x) for x in frameReferenceAll]
+    topSwath = min(refSwaths, key = lambda x: x.sensingStart)
+    botSwath = max(refSwaths, key = lambda x: x.sensingStop)
+    leftSwath = min(refSwaths, key = lambda x: x.nearRange)
+    rightSwath = max(refSwaths, key = lambda x: x.farRange)
+
+    totalWidth  = int(np.round((rightSwath.farRange - leftSwath.nearRange)/leftSwath.dr + 1))
+    totalLength = int(np.round((botSwath.sensingStop - topSwath.sensingStart).total_seconds()/topSwath.dt + 1 ))
+
+    #!!!should CHECK if input ionospheric shift is really in this geometry
+    nearRange = leftSwath.nearRange
+    sensingStart = topSwath.sensingStart
+    dr = leftSwath.dr
+    dt = topSwath.dt
+
+    for swath in swathList:
+        frameReference = ut.loadProduct(os.path.join(inps.reference_stack, 'IW{0}.xml'.format(swath)))
+        swathObj = Swath(frameReference)
+
+        minBurst = frameReference.bursts[0].burstNumber
+        maxBurst = frameReference.bursts[-1].burstNumber
+        #!!!should CHECK if input ionospheric shift is really in this geometry
+        if minBurst==maxBurst:
+            print('Skipping processing of swath {0}'.format(swath))
+            continue
+
+        midBurstIndex = round((minBurst + maxBurst) / 2.0) - minBurst
+        midBurst = frameReference.bursts[midBurstIndex]
+
+        outputDir = os.path.join(inps.output, 'IW{0}'.format(swath))
+        os.makedirs(outputDir, exist_ok=True)
+
+        #compute mean offset: indices start from 0
+        firstLine   = int(np.round((swathObj.sensingStart - sensingStart).total_seconds()/(dt*inps.nalks))) + 1
+        lastLine    = int(np.round((swathObj.sensingStop - sensingStart).total_seconds()/(dt*inps.nalks))) - 1
+        firstColumn = int(np.round((swathObj.nearRange - nearRange)/(dr*inps.nrlks)))     + 1
+        lastColumn  = int(np.round((swathObj.farRange - nearRange)/(dr*inps.nrlks)))     - 1
+        ionShiftSwath = ionShift[firstLine:lastLine+1, firstColumn:lastColumn+1]
+        ionShiftSwathValid = ionShiftSwath[np.nonzero(ionShiftSwath!=0)]
+
+        if ionShiftSwathValid.size == 0:
+            ionShiftSwathMean = 0.0
+            print('mean azimuth ionospheric shift of swath {}: {} single look azimuth lines'.format(swath, ionShiftSwathMean))
+        else:
+            ionShiftSwathMean = np.mean(ionShiftSwathValid, dtype=np.float64)
+            print('mean azimuth ionospheric shift of swath {}: {} single look azimuth lines'.format(swath, ionShiftSwathMean))
+
+        for burst in frameReference.bursts:
+            #compute mean offset: indices start from 0
+            firstLine   = int(np.round((burst.sensingStart - sensingStart).total_seconds()/(dt*inps.nalks))) + 1
+            lastLine    = int(np.round((burst.sensingStop - sensingStart).total_seconds()/(dt*inps.nalks))) - 1
+            firstColumn = int(np.round((burst.startingRange - nearRange)/(dr*inps.nrlks))) + 1
+            lastColumn  = int(np.round((burst.startingRange + (burst.numberOfSamples - 1) * dr - nearRange)/(dr*inps.nrlks))) - 1
+            ionShiftBurst = ionShift[firstLine:lastLine+1, firstColumn:lastColumn+1]
+
+            ionShiftBurstValid = ionShiftBurst[np.nonzero(ionShiftBurst!=0)]
+            if ionShiftBurstValid.size < (lastLine - firstLine + 1) * (lastColumn - firstColumn + 1) / 2.0:
+                ionShiftBurstMean = 0.0
+                print('mean azimuth ionospheric shift of burst {}: 0.0 single look azimuth lines'.format(burst.burstNumber))
+            else:
+                #ionShiftBurstMean should use both phaseRamp1 and phaseRamp2, while
+                #ionShiftSwathMean should use phaseRamp2 only as in (to be consistent with) previous ESD
+                #The above is tested against runIon.py in topsApp.py
+                #ionShiftBurstMean = np.mean(ionShiftBurstValid, dtype=np.float64) - ionShiftSwathMean
+                ionShiftBurstMean = np.mean(ionShiftBurstValid, dtype=np.float64)
+                print('mean azimuth ionospheric shift of burst {}: {} single look azimuth lines'.format(burst.burstNumber, ionShiftBurstMean))
+
+            #compute burst phase ramps
+            (dopplerOffset, Ka) = computeDopplerOffset(burst, 1, burst.numberOfLines, 1, burst.numberOfSamples, nrlks=1, nalks=1)
+
+            #satellite height
+            satHeight = np.linalg.norm(burst.orbit.interpolateOrbit(burst.sensingMid, method='hermite').getPosition())
+            #orgininal doppler offset should be multiplied by this ratio
+            ratio = ionHeight/(satHeight-earthRadius)
+
+            #phaseRamp1 and phaseRamp2 have same sign according to Liang et. al., 2019.
+            #phase ramp due to imaging geometry
+            phaseRamp1 = dopplerOffset * burst.azimuthTimeInterval * ratio * (burst.azimuthTimeInterval * Ka[None,:] * 4.0 * np.pi)
+            #phase ramp due to non-zero doppler centroid frequency, ESD
+            phaseRamp2 = dopplerOffset * burst.azimuthTimeInterval * Ka[None,:] * 2.0 * np.pi * burst.azimuthTimeInterval
+            phaseRamp = (phaseRamp1 + phaseRamp2) * ionShiftBurstMean - phaseRamp2 * ionShiftSwathMean
+
+            outfile = os.path.join(outputDir, '%s_%02d.float'%('burst', burst.burstNumber))
+
+            phaseRamp.astype(np.float32).tofile(outfile)
+
+            image = isceobj.createImage()
+            image.setDataType('FLOAT')
+            image.setFilename(outfile)
+            image.extraFilename = outfile + '.vrt'
+            image.setWidth(burst.numberOfSamples)
+            image.setLength(burst.numberOfLines)
+            image.renderHdr()
+
+
+
+if __name__ == '__main__':
+    '''
+    Main driver.
+    '''
+    # Main Driver
+    main()
+
+
+