prep_hyp3: update to be compatible with current hyp3 products (insarl…

…ab#573)

- removes hard coding of EARTH_RADIUS and HEIGHT
- enables use of HyP3 DEM and incidence angle products
- reorganizes main() and add_hyp3_metadata() functions
- merge add_sentinel1_metadata() into add_hyp3_metadata()
- updates help information
- update hyp3 related docs

docs/README.md

@@ -8,7 +8,7 @@
 
 ## MintPy ##
 
-The Miami INsar Time-series software in PYthon (MintPy as /mɪnt paɪ/) is an open-source package for Interferometric Synthetic Aperture Radar (InSAR) time series analysis. It reads the stack of interferograms (coregistered and unwrapped) in [ISCE](https://github.com/isce-framework/isce2), [ARIA](https://github.com/aria-tools/ARIA-tools), [FRInGE](https://github.com/isce-framework/fringe), [GMTSAR](https://github.com/gmtsar/gmtsar), [SNAP](http://step.esa.int/), [GAMMA](https://www.gamma-rs.ch/no_cache/software.html), [HyP3](https://hyp3-docs.asf.alaska.edu/) or ROI_PAC format, and produces three dimensional (2D in space and 1D in time) ground surface displacement in line-of-sight direction. It includes a routine time series analysis (`smallbaselineApp.py`) and some independent toolbox.
+The Miami INsar Time-series software in PYthon (MintPy as /mɪnt paɪ/) is an open-source package for Interferometric Synthetic Aperture Radar (InSAR) time series analysis. It reads the stack of interferograms (coregistered and unwrapped) in [ISCE](https://github.com/isce-framework/isce2), [ARIA](https://github.com/aria-tools/ARIA-tools), [FRInGE](https://github.com/isce-framework/fringe), [HyP3](https://hyp3-docs.asf.alaska.edu/), [GMTSAR](https://github.com/gmtsar/gmtsar), [SNAP](http://step.esa.int/), [GAMMA](https://www.gamma-rs.ch/no_cache/software.html) or ROI_PAC format, and produces three dimensional (2D in space and 1D in time) ground surface displacement in line-of-sight direction. It includes a routine time series analysis (`smallbaselineApp.py`) and some independent toolbox.
 
 This package was called PySAR before version 1.1.1. For version 1.1.2 and onward, we use MintPy instead.
 

docs/dir_structure.md

@@ -378,6 +378,43 @@ mintpy.load.azAngleFile      = $DATA_DIR/SanFranSenDT42/azimuthAngle/*.vrt
 mintpy.load.waterMaskFile    = $DATA_DIR/SanFranSenDT42/mask/watermask.msk
 ```
 
+### [ASF HyP3](https://hyp3-docs.asf.alaska.edu/)
+
+**WARNING: The current incidence angle file offered by HyP3 is not compatible with MintPy and should not be used!**
+1. Request and download GUNW products using [hyp3_sdk](https://nbviewer.jupyter.org/github/ASFHyP3/hyp3-sdk/blob/main/docs/sdk_example.ipynb).
+2. For at least one GUNW product, download the accompanying DEM.
+3. Clip DEM and all interferograms to the same area using hyp3lib/[cutGeotiffs.py](https://github.com/ASFHyP3/hyp3-lib/blob/develop/hyp3lib/cutGeotiffs.py) script.
+
+```
+$DATA_DIR/TongariroSenA
+├── hyp3
+|   ├── S1BA_20161229T070618_20170116T070658_VVP018_INT80_G_ueF_8108
+│   │   ├── S1BA_20161229T070618_20170116T070658_VVP018_INT80_G_ueF_8108_unw_phase_clip.tif
+│   │   ├── S1BA_20161229T070618_20170116T070658_VVP018_INT80_G_ueF_8108_corr_clip.tif
+│   │   ├── S1BA_20161229T070618_20170116T070658_VVP018_INT80_G_ueF_8108_dem_clip.tif
+│   │   ├── S1BA_20161229T070618_20170116T070658_VVP018_INT80_G_ueF_8108.txt
+│   │   └── ...
+│   ├── S1AB_20170116T070658_20170122T070616_VVP006_INT80_G_ueF_0209
+│   │   ├── S1AB_20170116T070658_20170122T070616_VVP006_INT80_G_ueF_0209_unw_phase_clip.tif
+│   │   ├── S1AB_20170116T070658_20170122T070616_VVP006_INT80_G_ueF_0209_corr_clip.tif
+│   │   ├── S1AB_20170116T070658_20170122T070616_VVP006_INT80_G_ueF_0209.txt
+│   │   └── ...
+│   └── ...
+└── mintpy
+    └── TongariroSenA.txt
+```
+
+The corresponding template options for `load_data`:
+
+```cfg
+mintpy.load.processor        = hyp3
+##---------interferogram datasets:
+mintpy.load.unwFile          = $DATA_DIR/TongariroSenA/hyp3/*/*unw_phase_clip.tif
+mintpy.load.corFile          = $DATA_DIR/TongariroSenA/hyp3/*/*corr_clip.tif
+##---------geometry datasets:
+mintpy.load.demFile          = $DATA_DIR/TongariroSenA/hyp3/*/*dem_clip.tif
+```
+
 ### [GMTSAR](https://github.com/gmtsar/gmtsar) ###
 
 Below is a recipe to prepare a stack of interferograms from Sentinel-1:
@@ -504,41 +541,6 @@ mintpy.load.corFile          = $DATA_DIR/WCapeSenAT29/interferograms/*/*/coh_*.i
 mintpy.load.demFile          = $DATA_DIR/WCapeSenAT29/dem_tc.data/dem*.img
 ```
 
-### [ASF HyP3](https://hyp3-docs.asf.alaska.edu/)
-
-1. Request and download GUNW products using [hyp3_sdk](https://nbviewer.jupyter.org/github/ASFHyP3/hyp3-sdk/blob/main/docs/sdk_example.ipynb).
-2. Download the corresponding DEM used in processing using [hyp3lib](https://github.com/ASFHyP3/hyp3-lib)/[getDEMfor.getDemFile()](https://github.com/ASFHyP3/hyp3-lib/blob/develop/hyp3lib/getDemFor.py#L16) function.
-3. Paste HyP3 interferogram metadata file (e.g. S1BB_20170510T070618_20170522T070619_VVP012_INT80_G_ueF_FF85.txt) into the same directory as your dem and give it the same name as your dem (e.g. dem.txt)
-4. Clip DEM and all interferograms to the same area using hyp3lib/[cutGeotiffs.py](https://github.com/ASFHyP3/hyp3-lib/blob/develop/hyp3lib/cutGeotiffs.py) script.
-
-```
-$DATA_DIR/TongariroSen
-├── DEM
-│   ├── dem.tif
-│   ├── dem_clip.tif
-│   └── dem_clip.txt
-├── S1BB_20170510T070618_20170522T070619_VVP012_INT80_G_ueF_FF85
-│   ├── S1BB_20170510T070618_20170522T070619_VVP012_INT80_G_ueF_FF85_unw_phase_clip.tif
-│   ├── S1BB_20170510T070618_20170522T070619_VVP012_INT80_G_ueF_FF85_corr_clip.tif
-│   ├── S1BB_20170510T070618_20170522T070619_VVP012_INT80_G_ueF_FF85.txt
-│   └── ...
-├── S1BB_20170428T070618_20170522T070619_VVP024_INT80_G_ueF_0CE0
-├── ...
-└── mintpy
-    └── TongariroSen.txt
-```
-
-The corresponding template options for `load_data`:
-
-```cfg
-mintpy.load.processor        = hyp3
-##---------interferogram datasets:
-mintpy.load.unwFile          = $DATA_DIR/TongariroSen/*/*unw_phase_clip.tif
-mintpy.load.corFile          = $DATA_DIR/TongariroSen/*/*corr_clip.tif
-##---------geometry datasets:
-mintpy.load.demFile          = $DATA_DIR/TongariroSen/dem_clip.tif
-```
-
 ### ROI_PAC (rsmas version) ###
 
 ```

mintpy/defaults/smallbaselineApp.cfg

@@ -25,7 +25,7 @@ mintpy.compute.config    = auto #[none / slurm / pbs / lsf ], auto for none (sam
 ## no   - save   0% disk usage, fast [default]
 ## lzf  - save ~57% disk usage, relative slow
 ## gzip - save ~62% disk usage, very slow [not recommend]
-mintpy.load.processor      = auto  #[isce, aria, gmtsar, snap, gamma, hyp3, roipac], auto for isce
+mintpy.load.processor      = auto  #[isce, aria, hyp3, gmtsar, snap, gamma, roipac], auto for isce
 mintpy.load.autoPath       = auto  #[yes / no], auto for no, use pre-defined auto path
 mintpy.load.updateMode     = auto  #[yes / no], auto for yes, skip re-loading if HDF5 files are complete
 mintpy.load.compression    = auto  #[gzip / lzf / no], auto for no.

mintpy/load_data.py

@@ -32,7 +32,7 @@
 
 
 #################################################################
-PROCESSOR_LIST = ['isce', 'aria', 'gmtsar', 'snap', 'gamma', 'hyp3', 'roipac']
+PROCESSOR_LIST = ['isce', 'aria', 'hyp3', 'gmtsar', 'snap', 'gamma', 'roipac']
 
 datasetName2templateKey = {
     'unwrapPhase'     : 'mintpy.load.unwFile',
@@ -505,7 +505,7 @@ def read_inps_dict2geometry_dict_object(iDict):
         # for processors with lookup table in geo-coordinates, remove latitude/longitude
         iDict['ds_name2key'].pop('latitude')
         iDict['ds_name2key'].pop('longitude')
-    elif iDict['processor'] in ['aria', 'gmtsar', 'snap', 'hyp3']:
+    elif iDict['processor'] in ['aria', 'gmtsar', 'hyp3', 'snap']:
         # for processors with geocoded products support only, do nothing for now.
         # check again when adding products support in radar-coordiantes
         pass

mintpy/prep_hyp3.py

@@ -15,17 +15,15 @@
 
 
 SPEED_OF_LIGHT = 299792458  # m/s
-EARTH_RADIUS = 6371e3       # Earth radius in meters
 
 
+#########################################################################
 EXAMPLE = """example:
-  prep_hyp3.py  interferograms/*/*unw_phase.tif
-  prep_hyp3.py  interferograms/*/*corr.tif
   prep_hyp3.py  interferograms/*/*unw_phase_clip.tif
   prep_hyp3.py  interferograms/*/*corr_clip.tif
+  prep_hyp3.py  interferograms/*/*dem_clip.tif
+  prep_hyp3.py  interferograms/*/*inc_map_clip.tif
   prep_hyp3.py  interferograms/*/*clip.tif
-  prep_hyp3.py  dem.tif
-  prep_hyp3.py  dem_clip.tif
 """
 
 DESCRIPTION = """
@@ -34,8 +32,9 @@
   2) S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr.tif
   3) S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2.txt
 
-  A DEM is also needed e.g.:
-  1) dem.tif
+  A DEM filename is needed and a incidence angle filename is recommended  e.g.:
+  1) S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_dem.tif
+  2) S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_inc_map.tif
 
   This script will read these files, read the geospatial metadata from GDAL,
   find the corresponding HyP3 metadata file (for interferograms and coherence),
@@ -46,32 +45,33 @@
 
   Before loading:
       For each interferogram, 3 files are needed:
-          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_unw_phase.tif
-          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr.tif
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_unw_phase_clip.tif
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr_clip.tif
           S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2.txt
-      For the geometry file 1 file is needed:
-          dem.tif (a DEM with any name)
+      For the geometry file 2 file are recommended:
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_dem_clip.tif     (required)
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_inc_map_clip.tif (optional but recommended)
 
   After running prep_hyp3.py:
       For each interferogram:
-          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_unw_phase.tif
-          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_unw_phase.tif.rsc
-          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr.tif
-          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr.tif.rsc
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_unw_phase_clip.tif
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_unw_phase_clip.tif.rsc
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr_clip.tif
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr_clip.tif.rsc
           S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2.txt
-      For the input DEM geometry file:
-          dem.tif
-          dem.tif.rsc
+      For the input geometry files:
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_dem_clip.tif
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_dem_clip.tif.rsc
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_inc_map_clip.tif
+          S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_inc_map_clip.tif.rsc
 
   Notes:
     HyP3 currently only supports generation of Sentinel-1 interferograms, so
     some Sentinel-1 metadata is hard-coded. If HyP3 adds processing of interferograms
     from other satellites, changes will be needed. 
 """
 
-#########################################################################
 
-# CMD parsing
 def create_parser():
     parser = argparse.ArgumentParser(description='Prepare attributes file for HyP3 InSAR product.\n'+
                                      DESCRIPTION,
@@ -90,82 +90,40 @@ def cmd_line_parse(iargs=None):
 
 
 #########################################################################
-# extract data from HyP3 interferogram metadata
 def add_hyp3_metadata(fname,meta,is_ifg=True):
     '''Read/extract attribute data from HyP3 metadata file and add to metadata dictionary
     Inputs:
-        *unw_phase.tif or *corr.tif file name, *dem.tif e.g.
-            S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_unw_phase.tif
-            S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr.tif
-            dem.tif
+        *unw_phase.tif, *corr.tif file name, *dem.tif, *inc_map.tif, e.g.
+            S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_unw_phase_clip.tif
+            S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_corr_clip.tif
+            S1AA_20161223T070700_20170116T070658_VVP024_INT80_G_ueF_74C2_dem_clip.tif
         Metadata dictionary (meta)
     Output:
         Metadata dictionary (meta)
     '''
 
-    if is_ifg:
-        # determine interferogram pair info and hyp3 metadata file name
-        sat, date1_string, date2_string, pol, res, soft, proc, ids, *_ = os.path.basename(fname).split('_')
-
-        job_id = '_'.join([sat, date1_string, date2_string, pol, res, soft, proc, ids])
-        directory = os.path.dirname(fname)
-        meta_file = f'{os.path.join(directory,job_id)}.txt'
-
-        date1 = datetime.strptime(date1_string,'%Y%m%dT%H%M%S')
-        date2 = datetime.strptime(date2_string,'%Y%m%dT%H%M%S')
-        date_avg = date1 + (date2 - date1) / 2
-        date_avg_seconds = (date_avg - date_avg.replace(hour=0, minute=0, second=0, microsecond=0)).total_seconds()
-
-        # open and read hyp3 metadata
-        hyp3_meta = {}
-        with open(meta_file, 'r') as f:
-            for line in f:
-                key, value = line.strip().split(': ')
-                hyp3_meta[key] = value
-
-        # add relevant data to meta object
-        meta['CENTER_LINE_UTC'] = date_avg_seconds
-        meta['DATE12'] = f'{date1.strftime("%y%m%d")}-{date2.strftime("%y%m%d")}'
-        meta['HEADING'] = hyp3_meta['Heading']
-        meta['ALOOKS'] = hyp3_meta['Azimuth looks']
-        meta['RLOOKS'] = hyp3_meta['Range looks']
-        meta['P_BASELINE_TOP_HDR'] = hyp3_meta['Baseline']
-        meta['P_BASELINE_BOTTOM_HDR'] = hyp3_meta['Baseline']
-
-    else:
-        meta_file = os.path.splitext(fname)[0]+'.txt'
-
-        # open and read hyp3 metadata
-        hyp3_meta = {}
-        with open(meta_file, 'r') as f:
-            for line in f:
-                key, value = line.strip().split(': ')
-                hyp3_meta[key] = value
-
-        # add relevant data to meta object
-        meta['HEADING'] = hyp3_meta['Heading']
-
-    return(meta)
-
-
-# add sentinel-1 metadata and remove unnessecary metadata
-def add_sentinel1_metadata(meta):
-    '''Add Sentinel-1 attribute data to metadata dictionary
-    Inputs:
-        Metadata dictionary (meta)
-    Output:
-        Metadata dictionary (meta)
-    '''
+    # determine interferogram pair info and hyp3 metadata file name
+    sat, date1_string, date2_string, pol, res, soft, proc, ids, *_ = os.path.basename(fname).split('_')
+    job_id = '_'.join([sat, date1_string, date2_string, pol, res, soft, proc, ids])
+    directory = os.path.dirname(fname)
+    meta_file = f'{os.path.join(directory,job_id)}.txt'
 
-    # note: HyP3 currently only supports Sentinel-1 data, so Sentinel-1
-    #       configuration is hard-coded.
+    # open and read hyp3 metadata
+    hyp3_meta = {}
+    with open(meta_file, 'r') as f:
+        for line in f:
+            key, value = line.strip().split(': ')
+            hyp3_meta[key] = value
 
+    # add universal hyp3 metadata
     meta['PROCESSOR'] = 'hyp3'
-    meta['PLATFORM'] = 'Sen'
-    meta['ANTENNA_SIDE'] = -1
-    meta['WAVELENGTH'] = SPEED_OF_LIGHT / sensor.SEN['carrier_frequency']
-    meta['HEIGHT'] = 693000.0 # nominal altitude of Sentinel1 orbit
-    meta['EARTH_RADIUS'] = EARTH_RADIUS
+    meta['HEADING'] = hyp3_meta['Heading']
+    meta['ALOOKS'] = hyp3_meta['Azimuth looks']
+    meta['RLOOKS'] = hyp3_meta['Range looks']
+    meta['P_BASELINE_TOP_HDR'] = hyp3_meta['Baseline']
+    meta['P_BASELINE_BOTTOM_HDR'] = hyp3_meta['Baseline']
+    meta['EARTH_RADIUS'] = hyp3_meta['Earth radius at nadir']
+    meta['HEIGHT'] = hyp3_meta['Spacecraft height']
 
     # add LAT/LON_REF1/2/3/4 based on whether satellite ascending or descending
     meta['ORBIT_DIRECTION'] = 'ASCENDING' if float(meta['HEADING']) > -90 else 'DESCENDING'
@@ -174,7 +132,6 @@ def add_sentinel1_metadata(meta):
     S = N + float(meta['Y_STEP']) * int(meta['LENGTH'])
     E = W + float(meta['X_STEP']) * int(meta['WIDTH'])
 
-    # NOTE by ZY, 27-03-2021: this should be converted from meters to degrees for pyaps to use it properly.
     if meta['ORBIT_DIRECTION'] == 'ASCENDING':
         meta['LAT_REF1'] = str(S)
         meta['LAT_REF2'] = str(S)
@@ -194,25 +151,37 @@ def add_sentinel1_metadata(meta):
         meta['LON_REF3'] = str(E)
         meta['LON_REF4'] = str(W)
 
+    # note: HyP3 currently only supports Sentinel-1 data, so Sentinel-1
+    #       configuration is hard-coded.
+    meta['PLATFORM'] = 'Sen'
+    meta['ANTENNA_SIDE'] = -1
+    meta['WAVELENGTH'] = SPEED_OF_LIGHT / sensor.SEN['carrier_frequency']
+
+    # add metadata that is only relevant to interferogram files
+    if is_ifg:
+        date1 = datetime.strptime(date1_string,'%Y%m%dT%H%M%S')
+        date2 = datetime.strptime(date2_string,'%Y%m%dT%H%M%S')
+        date_avg = date1 + (date2 - date1) / 2
+        date_avg_seconds = (date_avg - date_avg.replace(hour=0, minute=0, second=0, microsecond=0)).total_seconds()
+
+        meta['CENTER_LINE_UTC'] = date_avg_seconds
+        meta['DATE12'] = f'{date1.strftime("%y%m%d")}-{date2.strftime("%y%m%d")}'
+
     return(meta)
 
 
 #########################################################################
-
 def main(iargs=None):
     # read in arguments
     inps = cmd_line_parse(iargs)
 
     # for each filename, generate metadata rsc file
     for fname in inps.file:
+        is_ifg = any([x in fname for x in ['unw_phase','corr']])
         meta = readfile.read_gdal_vrt(fname)
+        meta = add_hyp3_metadata(fname, meta, is_ifg=is_ifg)
 
-        if any([x in fname for x in ['unw_phase','corr','inc_map']]):
-            meta = add_hyp3_metadata(fname, meta, is_ifg=True)
-        else:
-            meta = add_hyp3_metadata(fname, meta, is_ifg=False)
-
-        meta = add_sentinel1_metadata(meta)
+        # write
         rsc_file = fname+'.rsc'
         writefile.write_roipac_rsc(meta, out_file=rsc_file)