Geogrid update for post-processing correction

testGeogrid_ISCE.py

@@ -67,6 +67,12 @@ def cmdLineParse():
             help='Input stable surface mask')
     parser.add_argument('-fo', '--flag_optical', dest='optical_flag', type=bool, required=False, default=0,
             help='flag for reading optical data (e.g. Landsat): use 1 for on and 0 (default) for off')
+    parser.add_argument('-b', '--buffer', dest='buffer', type=bool, required=False, default=0,
+            help='buffer to add to the starting/end range accounting for all passes from the same relative orbit')
+    parser.add_argument('-p', '--parse', dest='parse', action='store_true',
+            default=False, help='Parse the SAFE zip file to get radar image and orbit metadata; no need to run ISCE')
+    parser.add_argument('--orbit-dir', hep='Directory of Sentinel-1 orbit files')
+    parser.add_argument('--aux-dir', hep='Directory of Sentinel-1 aux files')
 
     return parser.parse_args()
 
@@ -113,7 +119,7 @@ def getMergedOrbit(product):
     return orb
 
 
-def loadMetadata(indir):
+def loadMetadata(indir,buffer=0):
     '''
     Input file.
     '''
@@ -135,13 +141,78 @@ def loadMetadata(indir):
     info.prf = 1.0 / frames[0].bursts[0].azimuthTimeInterval
     info.rangePixelSize = frames[0].bursts[0].rangePixelSize
     info.lookSide = -1
+
+    info.startingRange -= buffer * info.rangePixelSize
+    info.farRange += buffer * info.rangePixelSize
+
     info.numberOfLines = int( np.round( (info.sensingStop - info.sensingStart).total_seconds() * info.prf)) + 1
-    info.numberOfSamples = int( np.round( (info.farRange - info.startingRange)/info.rangePixelSize)) + 1
+    info.numberOfSamples = int( np.round( (info.farRange - info.startingRange)/info.rangePixelSize)) + 1  + 2 * buffer
     info.orbit = getMergedOrbit(frames)
 
     return info
 
 
+def get_polarizations(s1_safe: str):
+    '''
+    Return:
+        Tuple[str]
+    '''
+    from typing import Tuple
+    from pathlib import Path
+
+    mapping = {
+        'SH': ('hh',),
+        'SV': ('vv',),
+        'DH': ('hh', 'hv'),
+        'DV': ('vv', 'vh'),
+    }
+    key = Path(s1_safe).name[14:16]
+
+    return mapping[key]
+
+
+def loadParsedata(indir, orbit_dir, aux_dir, buffer=0):
+    '''
+    Input file.
+    '''
+    import os
+    import numpy as np
+    import isce
+    from isceobj.Sensor.TOPS.Sentinel1 import Sentinel1
+
+
+    frames = []
+    for swath in range(1,4):
+        rdr=Sentinel1()
+        rdr.configure()
+#        rdr.safe=['./S1A_IW_SLC__1SDH_20180401T100057_20180401T100124_021272_024972_8CAF.zip']
+        rdr.safe=[indir]
+        rdr.output='reference'
+        rdr.orbitDir=orbit_dir
+        rdr.auxDir=aux_dir
+        rdr.swathNumber=swath
+        rdr.polarization=get_polarizations(indir)[0]
+        rdr.parse()
+        frames.append(rdr.product)
+
+    info = Dummy()
+    info.sensingStart = min([x.sensingStart for x in frames])
+    info.sensingStop = max([x.sensingStop for x in frames])
+    info.startingRange = min([x.startingRange for x in frames])
+    info.farRange = max([x.farRange for x in frames])
+    info.prf = 1.0 / frames[0].bursts[0].azimuthTimeInterval
+    info.rangePixelSize = frames[0].bursts[0].rangePixelSize
+    info.lookSide = -1
+
+    info.startingRange -= buffer * info.rangePixelSize
+    info.farRange += buffer * info.rangePixelSize
+
+    info.numberOfLines = int( np.round( (info.sensingStop - info.sensingStart).total_seconds() * info.prf)) + 1
+    info.numberOfSamples = int( np.round( (info.farRange - info.startingRange)/info.rangePixelSize)) + 1 + 2 * buffer
+    info.orbit = getMergedOrbit(frames)
+
+    return info
+
 def coregisterLoadMetadataOptical(indir_m, indir_s):
     '''
     Input file.
@@ -383,8 +454,12 @@ def main():
         metadata_m, metadata_s = coregisterLoadMetadataOptical(inps.indir_m, inps.indir_s)
         runGeogridOptical(metadata_m, metadata_s, inps.demfile, inps.dhdxfile, inps.dhdyfile, inps.vxfile, inps.vyfile, inps.srxfile, inps.sryfile, inps.csminxfile, inps.csminyfile, inps.csmaxxfile, inps.csmaxyfile, inps.ssmfile)
     else:
-        metadata_m = loadMetadata(inps.indir_m)
-        metadata_s = loadMetadata(inps.indir_s)
+        if inps.parse:
+            metadata_m = loadParsedata(inps.indir_m, inps.orbit_dir, inputs.aux_dir, inps.buffer)
+            metadata_s = loadParsedata(inps.indir_s, inps.orbit_dir, inputs.aux_dir, inps.buffer)
+        else:
+            metadata_m = loadMetadata(inps.indir_m,inps.buffer)
+            metadata_s = loadMetadata(inps.indir_s,inps.buffer)
         runGeogrid(metadata_m, metadata_s, inps.demfile, inps.dhdxfile, inps.dhdyfile, inps.vxfile, inps.vyfile, inps.srxfile, inps.sryfile, inps.csminxfile, inps.csminyfile, inps.csmaxxfile, inps.csmaxyfile, inps.ssmfile)