From b231d6e80a9f9d40a0ea95a54f00f57b9a11c058 Mon Sep 17 00:00:00 2001 From: CunrenLiang <56097947+CunrenLiang@users.noreply.github.com> Date: Tue, 2 Jul 2024 07:19:55 +0800 Subject: [PATCH] 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 * Update contrib/stack/topsStack/stackSentinel.py Co-authored-by: Zhang Yunjun * Update contrib/stack/topsStack/stackSentinel.py Co-authored-by: Zhang Yunjun * Update contrib/stack/topsStack/stackSentinel.py Co-authored-by: Zhang Yunjun * Update contrib/stack/topsStack/stackSentinel.py Co-authored-by: Zhang Yunjun --------- Co-authored-by: Zhang Yunjun --- .../isceobj/Alos2Proc/runFrameOffset.py | 9 +- contrib/stack/topsStack/README.md | 26 +- contrib/stack/topsStack/Stack.py | 128 ++++++++- contrib/stack/topsStack/burstRampIon.py | 169 ++++++++++++ contrib/stack/topsStack/filtIonShift.py | 242 ++++++++++++++++++ contrib/stack/topsStack/invertIon.py | 6 +- contrib/stack/topsStack/ion_param.txt | 3 + contrib/stack/topsStack/s1_select_ion.py | 7 +- contrib/stack/topsStack/stackSentinel.py | 58 ++++- 9 files changed, 631 insertions(+), 17 deletions(-) create mode 100755 contrib/stack/topsStack/burstRampIon.py create mode 100755 contrib/stack/topsStack/filtIonShift.py diff --git a/components/isceobj/Alos2Proc/runFrameOffset.py b/components/isceobj/Alos2Proc/runFrameOffset.py index 635e94c9..59356688 100644 --- 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 index 54b6ac7f..75f4454b 100644 --- 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 index 4fd860f5..580ac6c7 100644 --- 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) diff --git a/contrib/stack/topsStack/burstRampIon.py b/contrib/stack/topsStack/burstRampIon.py new file mode 100755 index 00000000..d2c810fc --- /dev/null +++ b/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() + + + diff --git a/contrib/stack/topsStack/filtIonShift.py b/contrib/stack/topsStack/filtIonShift.py new file mode 100755 index 00000000..3bec81c9 --- /dev/null +++ b/contrib/stack/topsStack/filtIonShift.py @@ -0,0 +1,242 @@ +#!/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 +from isceobj.TopsProc.runIon import adaptive_gaussian +from isceobj.TopsProc.runIon import weight_fitting +from isceobj.TopsProc.runIon import fit_surface +from isceobj.TopsProc.runIon import cal_surface +import s1a_isce_utils as ut + + +def createParser(): + parser = argparse.ArgumentParser(description='compute and filter azimuth ionospheric shift [unit: masBurst.azimuthTimeInterval]') + 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('-f', '--reference', type=str, dest='reference', required=True, + help='Directory with the reference image (coregistered or reference of the stack)') + parser.add_argument('-s', '--secondary', type=str, dest='secondary', required=True, + help='Directory with the secondary image (coregistered or reference of the stack)') + parser.add_argument('-i', '--input', dest='input', type=str, required=True, + help='input ionosphere') + parser.add_argument('-c', '--coherence', dest='coherence', type=str, required=True, + help='coherence, e.g. raw_no_projection.cor') + parser.add_argument('-o', '--output', dest='output', type=str, required=True, + help='output azimuth ionospheric shift') + parser.add_argument('-b', '--win_min', dest='win_min', type=int, default=100, + help='minimum filtering window size') + parser.add_argument('-d', '--win_max', dest='win_max', type=int, default=200, + help='maximum filtering window size') + parser.add_argument('-r', '--nrlks', dest='nrlks', type=int, default=1, + help='number of range looks. Default: 1') + parser.add_argument('-a', '--nalks', dest='nalks', type=int, default=1, + help='number of azimuth looks. Default: 1') + #parser.add_argument('-m', '--masked_areas', dest='masked_areas', type=int, nargs='+', action='append', default=None, + # help='This is a 2-d list. Each element in the 2-D list is a four-element list: [firstLine, lastLine, firstColumn, lastColumn], with line/column numbers starting with 1. If one of the four elements is specified with -1, the program will use firstLine/lastLine/firstColumn/lastColumn instead. e.g. two areas masked out: --masked_areas 10 20 10 20 --masked_areas 110 120 110 120') + + return parser + + +def cmdLineParse(iargs = None): + parser = createParser() + return parser.parse_args(args=iargs) + + +def main(iargs=None): + ''' + check overlap among all acquistions, only keep the bursts that in the common overlap, + and then renumber the bursts. + ''' + inps = cmdLineParse(iargs) + + ''' + calculate azimuth shift caused by ionosphere using ionospheric phase + ''' + + ################################################# + #SET PARAMETERS HERE + #gaussian filtering window size + #size = np.int(np.around(width / 12.0)) + #size = ionParam.ionshiftFilteringWinsize + size_max = inps.win_max + size_min = inps.win_min + + #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self) + #if applying polynomial fitting + #0: no fitting, 1: with fitting + fit = 0 + corThresholdIonshift = 0.85 + ################################################# + + +#################################################################### + #STEP 1. GET DERIVATIVE OF IONOSPHERE +#################################################################### + + #get files + ionfile = inps.input + #we are using filtered ionosphere, so we should use coherence file that is not projected. + #corfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCor) + corfile = inps.coherence + img = isceobj.createImage() + img.load(ionfile + '.xml') + width = img.width + length = img.length + amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :] + ion = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :] + cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :] + + ######################################################################################## + #AFTER COHERENCE IS RESAMPLED AT grd2ion, THERE ARE SOME WIRED VALUES + cor[np.nonzero(cor<0)] = 0.0 + cor[np.nonzero(cor>1)] = 0.0 + ######################################################################################## + + #get the azimuth derivative of ionosphere + dion = np.diff(ion, axis=0) + dion = np.concatenate((dion, np.zeros((1,width))), axis=0) + + #remove the samples affected by zeros + flag_ion0 = (ion!=0) + #moving down by one line + flag_ion1 = np.roll(flag_ion0, 1, axis=0) + flag_ion1[0,:] = 0 + #moving up by one line + flag_ion2 = np.roll(flag_ion0, -1, axis=0) + flag_ion2[-1,:] = 0 + #now remove the samples affected by zeros + flag_ion = flag_ion0 * flag_ion1 * flag_ion2 + dion *= flag_ion + + flag = flag_ion * (cor>corThresholdIonshift) + index = np.nonzero(flag) + + +#################################################################### + #STEP 2. FIT A POLYNOMIAL TO THE DERIVATIVE OF IONOSPHERE +#################################################################### + + order = 3 + + #look for data to use + point_index = np.nonzero(flag) + m = point_index[0].shape[0] + + #calculate input index matrix + x0=np.matlib.repmat(np.arange(width), length, 1) + y0=np.matlib.repmat(np.arange(length).reshape(length, 1), 1, width) + + x = x0[point_index].reshape(m, 1) + y = y0[point_index].reshape(m, 1) + z = dion[point_index].reshape(m, 1) + w = cor[point_index].reshape(m, 1) + + #convert to higher precision type before use + x=np.asfarray(x,np.float64) + y=np.asfarray(y,np.float64) + z=np.asfarray(z,np.float64) + w=np.asfarray(w,np.float64) + coeff = fit_surface(x, y, z, w, order) + + rgindex = np.arange(width) + azindex = np.arange(length).reshape(length, 1) + #convert to higher precision type before use + rgindex=np.asfarray(rgindex,np.float64) + azindex=np.asfarray(azindex,np.float64) + dion_fit = cal_surface(rgindex, azindex, coeff, order) + + #no fitting + if fit == 0: + dion_fit *= 0 + dion_res = (dion - dion_fit)*(dion!=0) + + +#################################################################### + #STEP 3. FILTER THE RESIDUAL OF THE DERIVATIVE OF IONOSPHERE +#################################################################### + + #this will be affected by low coherence areas like water, so not use this. + #filter the derivation of ionosphere + #if size % 2 == 0: + # size += 1 + #sigma = size / 2.0 + + #g2d = gaussian(size, sigma, scale=1.0) + #scale = ss.fftconvolve((dion_res!=0), g2d, mode='same') + #dion_filt = ss.fftconvolve(dion_res, g2d, mode='same') / (scale + (scale==0)) + + #minimize the effect of low coherence pixels + cor[np.nonzero( (cor<0.85)*(cor!=0) )] = 0.00001 + dion_filt = adaptive_gaussian(dion_res, cor, size_max, size_min) + + dion = (dion_fit + dion_filt)*(dion!=0) + + #return dion + + +#################################################################### + #STEP 4. CONVERT TO AZIMUTH SHIFT +#################################################################### + + #!!! use number of swaths in reference for now + + #these are coregistered secondaries, so there should be <= swaths in reference of entire stack + referenceSwathList = ut.getSwathList(inps.reference) + secondarySwathList = ut.getSwathList(inps.secondary) + swathList = list(sorted(set(referenceSwathList+secondarySwathList))) + + #swathList = ut.getSwathList(inps.reference) + firstSwath = None + for swath in swathList: + frameReference = ut.loadProduct(os.path.join(inps.reference_stack, 'IW{0}.xml'.format(swath))) + + minBurst = frameReference.bursts[0].burstNumber + maxBurst = frameReference.bursts[-1].burstNumber + if minBurst==maxBurst: + print('Skipping processing of swath {0}'.format(swath)) + continue + else: + firstSwath = swath + break + if firstSwath is None: + raise Exception('no valid swaths!') + + midBurstIndex = round((minBurst + maxBurst) / 2.0) - minBurst + masBurst = frameReference.bursts[midBurstIndex] + + #shift casued by ionosphere [unit: masBurst.azimuthTimeInterval] + rng = masBurst.rangePixelSize * ((np.arange(width))*inps.nrlks + (inps.nrlks - 1.0) / 2.0) + masBurst.startingRange + Ka = masBurst.azimuthFMRate(rng) + ionShift = dion / (masBurst.azimuthTimeInterval * inps.nalks) / (4.0 * np.pi) / Ka[None, :] / masBurst.azimuthTimeInterval + + #output + outfile = inps.output + os.makedirs(os.path.dirname(outfile), exist_ok=True) + ion = np.zeros((length*2, width), dtype=np.float32) + ion[0:length*2:2, :] = amp + ion[1:length*2:2, :] = ionShift + ion.astype(np.float32).tofile(outfile) + img.filename = outfile + img.extraFilename = outfile + '.vrt' + img.renderHdr() + + +if __name__ == '__main__': + ''' + Main driver. + ''' + # Main Driver + main() + + + diff --git a/contrib/stack/topsStack/invertIon.py b/contrib/stack/topsStack/invertIon.py index b3e5f1e9..43abd5f1 100755 --- a/contrib/stack/topsStack/invertIon.py +++ b/contrib/stack/topsStack/invertIon.py @@ -96,6 +96,8 @@ def cmdLineParse(): parser = argparse.ArgumentParser(description='least squares estimation') parser.add_argument('--idir', dest='idir', type=str, required=True, help = 'input directory where each pair (YYYYMMDD_YYYYMMDD) is located. only folders are recognized') + parser.add_argument('--filename', dest='filename', type=str, default='filt.ion', + help = 'input file name. default: filt.ion') parser.add_argument('--odir', dest='odir', type=str, required=True, help = 'output directory for estimated result of each date') parser.add_argument('--zro_date', dest='zro_date', type=str, default=None, @@ -205,7 +207,7 @@ def cmdLineParse(): ndate = len(dates) npair = len(pairs) - ionfile = os.path.join(idir, pairs[0], 'ion_cal', 'filt.ion') + ionfile = os.path.join(idir, pairs[0], 'ion_cal', inps.filename) img = isceobj.createImage() img.load(ionfile+'.xml') @@ -219,7 +221,7 @@ def cmdLineParse(): wls = False stdPairs = np.ones((npair, length, width), dtype=np.float32) for i in range(npair): - ionfile = os.path.join(idir, pairs[i], 'ion_cal', 'filt.ion') + ionfile = os.path.join(idir, pairs[i], 'ion_cal', inps.filename) ionPairs[i, :, :] = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :] #flag of valid/invalid is defined by amplitde image amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :] diff --git a/contrib/stack/topsStack/ion_param.txt b/contrib/stack/topsStack/ion_param.txt index 3896804d..a01bbbbf 100644 --- a/contrib/stack/topsStack/ion_param.txt +++ b/contrib/stack/topsStack/ion_param.txt @@ -2,6 +2,9 @@ ###the values below are the default values used by the module ###remove # to set the parameters +#consider burst properties in ionosphere computation: False +#height of ionosphere layer in km: 200.0 + #maximum window size for filtering ionosphere phase: 200 #minimum window size for filtering ionosphere phase: 100 #maximum window size for filtering ionosphere azimuth shift: 150 diff --git a/contrib/stack/topsStack/s1_select_ion.py b/contrib/stack/topsStack/s1_select_ion.py index ec517d1a..75d70551 100755 --- a/contrib/stack/topsStack/s1_select_ion.py +++ b/contrib/stack/topsStack/s1_select_ion.py @@ -422,8 +422,8 @@ def cmdLineParse(): parser = argparse.ArgumentParser( description='select Sentinel-1A/B acquistions good for ionosphere correction. not used slices are moved to folder: not_used') parser.add_argument('-dir', dest='dir', type=str, required=True, help = 'directory containing the "S1*_IW_SLC_*.zip" files') - parser.add_argument('-sn', dest='sn', type=str, required=True, - help='south/north bound of area of interest, format: south/north') + parser.add_argument('-sn', dest='sn', type=float, nargs=2, required=True, + help='south and north bound of area of interest, format: -sn south north') parser.add_argument('-nr', dest='nr', type=int, default=10, help = 'minimum number of acquisitions for same starting ranges. default: 10') @@ -438,7 +438,8 @@ def cmdLineParse(): if __name__ == '__main__': inps = cmdLineParse() - s,n=[float(x) for x in inps.sn.split('/')] + #s,n=[float(x) for x in inps.sn.split('/')] + s,n=inps.sn #group the slices group = get_group(inps.dir) diff --git a/contrib/stack/topsStack/stackSentinel.py b/contrib/stack/topsStack/stackSentinel.py index ef92360e..24639d18 100755 --- a/contrib/stack/topsStack/stackSentinel.py +++ b/contrib/stack/topsStack/stackSentinel.py @@ -14,7 +14,7 @@ import isce import isceobj from isceobj.Sensor.TOPS.Sentinel1 import Sentinel1 -from topsStack.Stack import config, run, sentinelSLC +from topsStack.Stack import config, run, sentinelSLC, ionParamUsr helpstr = """ @@ -607,8 +607,9 @@ def checkCurrentStatusIonosphere(inps): acquisitionDates, stackReferenceDate, secondaryDates, safe_dict = get_dates(inps) pairs_same_starting_ranges, pairs_diff_starting_ranges = selectNeighborPairsIonosphere(safe_dict, inps.num_connections_ion) + dateListIon = getDatesIonosphere(pairs_same_starting_ranges, pairs_diff_starting_ranges) pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update = excludeExistingPairsIonosphere(pairs_same_starting_ranges, pairs_diff_starting_ranges, inps.work_dir) - dateListIon = getDatesIonosphere(pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update) + dateListIonUpdate = getDatesIonosphere(pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update) #report pairs of different swath starting ranges. pdiff = 'ionosphere phase estimation pairs with different swath starting ranges\n' @@ -623,7 +624,7 @@ def checkCurrentStatusIonosphere(inps): with open('pairs_diff_starting_ranges.txt', 'w') as f: f.write(pdiff) - return dateListIon, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict + return dateListIonUpdate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, dateListIon ######################################## @@ -821,12 +822,12 @@ def offsetStack(inps, acquisitionDates, stackReferenceDate, secondaryDates, safe return i -def ionosphereStack(inps, dateListIon, stackReferenceDate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, i): +def ionosphereStack(inps, dateListIonUpdate, dateListIon, stackReferenceDate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, i): i+=1 runObj = run() runObj.configure(inps, 'run_{:02d}_subband_and_resamp'.format(i)) - runObj.subband_and_resamp(dateListIon, stackReferenceDate) + runObj.subband_and_resamp(dateListIonUpdate, stackReferenceDate) runObj.finalize() i+=1 @@ -871,6 +872,38 @@ def ionosphereStack(inps, dateListIon, stackReferenceDate, pairs_same_starting_r runObj.invertIon() runObj.finalize() + + ionParamUsrObj = ionParamUsr(inps.param_ion) + ionParamUsrObj.configure() + if ionParamUsrObj.ION_considerBurstProperties: + i+=1 + runObj = run() + runObj.configure(inps, 'run_{:02d}_filtIonShift'.format(i)) + runObj.filtIonShift(pairs_same_starting_ranges_update + pairs_diff_starting_ranges_update) + runObj.finalize() + + #From now on, need to reprocess all dates/pairs, not only newly added + i+=1 + runObj = run() + runObj.configure(inps, 'run_{:02d}_invertIonShift'.format(i)) + runObj.invertIonShift() + runObj.finalize() + + i+=1 + runObj = run() + runObj.configure(inps, 'run_{:02d}_burstRampIon'.format(i)) + runObj.burstRampIon(dateListIon) + runObj.finalize() + + i+=1 + runObj = run() + runObj.configure(inps, 'run_{:02d}_mergeBurstRampIon'.format(i)) + runObj.mergeBurstRampIon(dateListIon, stackReferenceDate) + runObj.finalize() + + del ionParamUsrObj + + return i @@ -1015,8 +1048,19 @@ def main(iargs=None): elif not os.path.isfile(inps.param_ion): print("Ion parameter file is missing. Ionospheric estimation will not be done.") else: - dateListIon, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict = checkCurrentStatusIonosphere(inps) - i = ionosphereStack(inps, dateListIon, stackReferenceDate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, i) + ############################################# + #check if consider burst properties + ionParamUsrObj = ionParamUsr(inps.param_ion) + ionParamUsrObj.configure() + if ionParamUsrObj.ION_considerBurstProperties: + if inps.coregistration != 'NESD': + raise Exception("Coregistration options '-C', '--coregistration' must be NESD when 'consider burst properties in ionosphere computation' is True") + del ionParamUsrObj + ############################################# + + dateListIonUpdate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, dateListIon = checkCurrentStatusIonosphere(inps) + i = ionosphereStack(inps, dateListIonUpdate, dateListIon, stackReferenceDate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, i) + if __name__ == "__main__": # Main engine