update gmtsar support for the new ifg stack structure (#1226)

This PR updates support for GMTSAR interferogram stack as prepared for the EarthScope short course 2024 by Xiaohua Xu.

+ utils.readfile: add GMTSAR_SENSOR_ID2NAME, and call it in _attribute_gmtsar2roipac() to get PLATFORM

+ utils.ptime: add yyyyddd2yyyymmdd() to convert from GMTSAR date format into YYYYMMDD format

+ prep_gmtsar: update to the new structure as prepared by Xiaohua Xu for the 2024 EarthScope short course. Major changes include: 
   - use supermaster.PRM for common metadata extraction;
   - use baseline_table.dat to read the bperp TS directly;
   - support inc/az angles files

+ docs/dir_structure: update to the latest structure.

docs/dir_structure.md

@@ -459,40 +459,43 @@ Below is a recipe to prepare a stack of interferograms from Sentinel-1:
 + [https://topex.ucsd.edu/gmtsar/tar/sentinel_time_series_2.pdf](https://topex.ucsd.edu/gmtsar/tar/sentinel_time_series_2.pdf)
 
 ```
-$DATA_DIR/StHelensEnvDT156
+$DATA_DIR/SanFranBaySenD42
+├── baseline_table.dat
+├── supermaster.PRM
 ├── geometry
-│   └── topo_ll.grd
+|   ├── azimuth_angle.grd
+|   ├── dem.grd
+│   └── incidence_angle.grd
 ├── interferograms
 │   ├── 2004114_2004324
-│   │   ├── 20040423.LED
-│   │   ├── 20040423.PRM
-│   │   ├── 20041119.LED
-│   │   ├── 20041119.PRM
-│   │   ├── baseline.txt   #generated by SAT_baseline
-│   │   ├── corr.grd
 │   │   ├── corr_ll.grd
-│   │   ├── unwrap.grd
 │   │   └── unwrap_ll.grd
 │   ├── 2004114_2004359
 │   └── ...
 └── mintpy
-    └── StHelensEnvDT156.txt
+    └── SanFranBaySenD42.txt
 ```
 
 The corresponding template options for `load_data`:
 
 ```cfg
 ## manually specify the following attributes since they are missing from gmtsar products
+ALOOKS          = 8          #[int], number of looks in the azimuth direction
+RLOOKS          = 32         #[int], number of looks in the range direction
 HEADING         = -168.0     #[float], satellite heading angle, measured from the north in clockwise as positive
                              # One could open the *.kml file in Google Earth and measure it manually
 ORBIT_DIRECTION = DESCENDING #[ASCENDING, DESCENDING]
 
-mintpy.load.processor   = gmtsar
+mintpy.load.processor    = gmtsar
+mintpy.load.metaFile     = $DATA_DIR/SanFranBaySenD42/supermaster.PRM
+mintpy.load.baselineDir  = $DATA_DIR/SanFranBaySenD42/baseline_table.dat
 ##---------interferogram datasets:
-mintpy.load.unwFile     = $DATA_DIR/StHelensEnvDT156/interferograms/*/unwrap_ll.grd
-mintpy.load.corFile     = $DATA_DIR/StHelensEnvDT156/interferograms/*/corr_ll.grd
+mintpy.load.unwFile      = $DATA_DIR/SanFranBaySenD42/interferograms/*/unwrap_ll*.grd
+mintpy.load.corFile      = $DATA_DIR/SanFranBaySenD42/interferograms/*/corr_ll*.grd
 ##---------geometry datasets:
-mintpy.load.demFile     = $DATA_DIR/StHelensEnvDT156/geometry/topo_ll.grd
+mintpy.load.demFile      = $DATA_DIR/SanFranBaySenD42/geometry/dem*.grd
+mintpy.load.incAngleFile = $DATA_DIR/SanFranBaySenD42/geometry/incidence_angle.grd
+mintpy.load.azAngleFile  = $DATA_DIR/SanFranBaySenD42/geometry/azimuth_angle.grd
 ```
 
 ### Gamma ###

pyproject.toml

@@ -110,7 +110,7 @@ include-package-data = true
 dependencies = { file = ["requirements.txt"] }
 readme = { file = ["docs/README.md"], content-type = "text/markdown" }
 
-# extra requirements: `pip install pysolid[test]` or `pip install .[test]`
+# extra requirements: `pip install mintpy[test]` or `pip install .[test]`
 [tool.setuptools.dynamic.optional-dependencies.test]
 file = ["tests/requirements.txt"]
 

src/mintpy/load_data.py

@@ -382,7 +382,7 @@ def read_inps_dict2ifgram_stack_dict_object(iDict, ds_name2template_key):
     pairsDict = {}
     for i, dsPath0 in enumerate(dsPathDict[dsName0]):
         # date string used in the file/dir path
-        # YYYYDDD       for gmtsar [modern Julian date]
+        # YYYYDDD       for gmtsar [day of the year - 1]
         # YYYYMMDDTHHMM for uavsar
         # YYYYMMDD      for all the others
         date6s = readfile.read_attribute(dsPath0)['DATE12'].replace('_','-').split('-')

src/mintpy/prep_gmtsar.py

@@ -19,6 +19,26 @@
 
 
 #########################################################################
+def get_ifg_file(ifg_dir, coord_type='geo', fbases=['corr', 'phase', 'phasefilt', 'unwrap']):
+    ifg_file = None
+
+    # search interferogram data file
+    for fbase in fbases:
+        if coord_type == 'geo':
+            fnames = glob.glob(os.path.join(ifg_dir, f'{fbase}_ll*.grd'))
+        else:
+            fnames = glob.glob(os.path.join(ifg_dir, f'{fbase}.grd'))
+
+        if len(fnames) > 0:
+            ifg_file = fnames[0]
+            break
+
+    if not ifg_file:
+        print(f'WARNING: No {coord_type}-coord files found in {ifg_dir} with suffix: {fbases}')
+
+    return ifg_file
+
+
 def get_prm_files(ifg_dir):
     """Get the date1/2.prm files in the interferogram directory."""
     prm_files = sorted(glob.glob(os.path.join(ifg_dir, '*.PRM')))
@@ -31,12 +51,10 @@ def get_prm_files(ifg_dir):
 def get_multilook_number(ifg_dir, fbases=['corr', 'phase', 'phasefilt', 'unwrap']):
     """Get the number of multilook in range and azimuth direction."""
     # grab an arbitrary file in radar-coordiantes
-    rdr_files = [os.path.join(ifg_dir, f'{i}.grd') for i in fbases]
-    if len(rdr_files) == 0:
-        raise ValueError(f'No radar-coord files found in {ifg_dir} with suffix: {fbases}')
+    rdr_file = get_ifg_file(ifg_dir, coord_type='radar')
 
     # read step info into multilook number
-    ds = gdal.Open(rdr_files[0], gdal.GA_ReadOnly)
+    ds = gdal.Open(rdr_file, gdal.GA_ReadOnly)
     transform = ds.GetGeoTransform()
     rlooks = int(abs(transform[1]))
     alooks = int(abs(transform[5]))
@@ -46,13 +64,12 @@ def get_multilook_number(ifg_dir, fbases=['corr', 'phase', 'phasefilt', 'unwrap'
 def get_lalo_ref(ifg_dir, prm_dict, fbases=['corr', 'phase', 'phasefilt', 'unwrap']):
     """Get the LAT/LON_REF1/2/3/4 from *_ll.grd file."""
     # grab an arbitrary file in geo-coordiantes
-    geo_files = [os.path.join(ifg_dir, f'{i}_ll.grd') for i in fbases]
-    if len(geo_files) == 0:
-        print(f'WARNING: No geo-coord files found in {ifg_dir} with suffix: {fbases}')
+    geo_file = get_ifg_file(ifg_dir, coord_type='geo')
+    if not geo_file:
         return prm_dict
 
     # read corners lat/lon info
-    ds = gdal.Open(geo_files[0], gdal.GA_ReadOnly)
+    ds = gdal.Open(geo_file, gdal.GA_ReadOnly)
     transform = ds.GetGeoTransform()
     x_step = abs(transform[1])
     y_step = abs(transform[5]) * -1.
@@ -88,13 +105,13 @@ def get_lalo_ref(ifg_dir, prm_dict, fbases=['corr', 'phase', 'phasefilt', 'unwra
 
 def get_slant_range_distance(ifg_dir, prm_dict, fbases=['corr', 'phase', 'phasefilt', 'unwrap']):
     """Get a constant slant range distance in the image center, for dataset in geo-coord."""
-    # grab an arbitrary file in radar-coordiantes
-    rdr_files = [os.path.join(ifg_dir, f'{i}.grd') for i in fbases]
-    if len(rdr_files) == 0:
+    # grab an arbitrary file in geo/radar-coordiantes
+    geo_file = get_ifg_file(ifg_dir, coord_type='geo')
+    if not geo_file:
         raise ValueError(f'No radar-coord files found in {ifg_dir} with suffix: {fbases}')
 
     # read width from rdr_file
-    ds = gdal.Open(rdr_files[0], gdal.GA_ReadOnly)
+    ds = gdal.Open(geo_file, gdal.GA_ReadOnly)
     width = ds.RasterXSize
 
     near_range = float(prm_dict['STARTING_RANGE'])
@@ -105,8 +122,30 @@ def get_slant_range_distance(ifg_dir, prm_dict, fbases=['corr', 'phase', 'phasef
     return prm_dict
 
 
+def read_baseline_table(fname):
+    """Read GMTSAR baseline table file.
+
+    Example: baseline_table.dat file
+    #  file_ID         yyyyddd.fraction day_cnt    b_para           b_perp
+    S1_20141231_ALL_F2 2014364.5885345591 364  -58.814087517826  -16.031404204594
+    S1_20150301_ALL_F2 2015059.5885222249 424  -33.273566547687  -17.102628888348
+    S1_20150325_ALL_F2 2015083.5885250464 448 -104.664278131966 -142.776284597889
+
+    Parameters: fname - str, path to the baseline table file
+    Returns:    bdict - dict, perpendicular baseline dictionary
+    """
+    print(f'read baseline time-series from file: {fname}')
+    fc = np.loadtxt(fname, dtype=str, usecols=(1,4))
+    date_list = ptime.yyyyddd2yyyymmdd([x[:7] for x in fc[:,0]])
+    bperp_list = fc[:,1].astype(np.float32)
+    bdict = {}
+    for date_str, bperp in zip(date_list, bperp_list):
+        bdict[date_str] = bperp
+    return bdict
+
+
 #########################################################################
-def extract_gmtsar_metadata(unw_file, template_file, rsc_file=None, update_mode=True):
+def extract_gmtsar_metadata(meta_file, unw_file, template_file, rsc_file=None, update_mode=True):
     """Extract metadata from GMTSAR interferogram stack."""
 
     # update_mode: check existing rsc_file
@@ -116,22 +155,22 @@ def extract_gmtsar_metadata(unw_file, template_file, rsc_file=None, update_mode=
     ifg_dir = os.path.dirname(unw_file)
 
     # 1. read *.PRM file
-    prm_file = get_prm_files(ifg_dir)[0]
-    meta = readfile.read_gmtsar_prm(prm_file)
+    #prm_file = get_prm_files(ifg_dir)[0]
+    meta = readfile.read_gmtsar_prm(meta_file)
     meta['PROCESSOR'] = 'gmtsar'
 
     # 2. read template file: HEADING, ORBIT_DIRECTION
     template = readfile.read_template(template_file)
-    for key in ['HEADING', 'ORBIT_DIRECTION']:
+    for key in ['HEADING', 'ORBIT_DIRECTION', 'ALOOKS', 'RLOOKS']:
         if key in template.keys():
             meta[key] = template[key].lower()
         else:
             raise ValueError('Attribute {} is missing! Please manually specify it in the template file.')
 
     # 3. grab A/RLOOKS from radar-coord data file
-    meta['ALOOKS'], meta['RLOOKS'] = get_multilook_number(ifg_dir)
-    meta['AZIMUTH_PIXEL_SIZE'] *= meta['ALOOKS']
-    meta['RANGE_PIXEL_SIZE'] *= meta['RLOOKS']
+    #meta['ALOOKS'], meta['RLOOKS'] = get_multilook_number(ifg_dir)
+    meta['AZIMUTH_PIXEL_SIZE'] *= int(meta['ALOOKS'])
+    meta['RANGE_PIXEL_SIZE'] *= int(meta['RLOOKS'])
 
     # 4. grab LAT/LON_REF1/2/3/4 from geo-coord data file
     meta = get_lalo_ref(ifg_dir, meta)
@@ -202,7 +241,7 @@ def prepare_geometry(geom_files, meta, update_mode=True):
     return
 
 
-def prepare_stack(unw_files, meta, update_mode=True):
+def prepare_stack(unw_files, meta, bperp_dict, update_mode=True):
     """Prepare .rsc file for all unwrapped interferogram files."""
     num_file = len(unw_files)
     if num_file == 0:
@@ -222,22 +261,15 @@ def prepare_stack(unw_files, meta, update_mode=True):
         ifg_meta.update(readfile.read_gdal_vrt(unw_file))
 
         # add DATE12
-        prm_files = get_prm_files(ifg_dir)
-        date1, date2 = (os.path.splitext(os.path.basename(i))[0] for i in prm_files)
+        date1, date2 = os.path.basename(ifg_dir).replace('-', '_').split('_')
+        date1 = ptime.yyyyddd2yyyymmdd(date1)
+        date2 = ptime.yyyyddd2yyyymmdd(date2)
         ifg_meta['DATE12'] = f'{ptime.yymmdd(date1)}-{ptime.yymmdd(date2)}'
 
         # and P_BASELINE_TOP/BOTTOM_HDR
-        baseline_file = os.path.join(ifg_dir, 'baseline.txt')
-        if os.path.isfile(baseline_file):
-            bDict = readfile.read_template(baseline_file)
-            ifg_meta['P_BASELINE_TOP_HDR'] = bDict['B_perpendicular']
-            ifg_meta['P_BASELINE_BOTTOM_HDR'] = bDict['B_perpendicular']
-        else:
-            ifg_meta['P_BASELINE_TOP_HDR'] = '0'
-            ifg_meta['P_BASELINE_BOTTOM_HDR'] = '0'
-            msg = f'WARNING: No baseline file found in: {baseline_file}. '
-            msg += 'Set P_BASELINE* to 0 and continue.'
-            print(msg)
+        bperp = bperp_dict[date2] - bperp_dict[date1]
+        ifg_meta['P_BASELINE_TOP_HDR'] = bperp
+        ifg_meta['P_BASELINE_BOTTOM_HDR'] = bperp
 
         # write .rsc file
         rsc_file = unw_file+'.rsc'
@@ -256,22 +288,30 @@ def prepare_stack(unw_files, meta, update_mode=True):
 def prep_gmtsar(inps):
     # read file path from template file
     template = readfile.read_template(inps.template_file)
+    inps.meta_file = glob.glob(template['mintpy.load.metaFile'])[0]
+    inps.bperp_file = glob.glob(template['mintpy.load.baselineDir'])[0]
     inps.unw_files = sorted(glob.glob(template['mintpy.load.unwFile']))
     inps.cor_files = sorted(glob.glob(template['mintpy.load.corFile']))
     inps.dem_file = glob.glob(template['mintpy.load.demFile'])[0]
+    inps.inc_angle_file = glob.glob(template['mintpy.load.incAngleFile'])[0]
+    inps.az_angle_file = glob.glob(template['mintpy.load.azAngleFile'])[0]
 
     # extract common metadata
-    rsc_file = os.path.join(inps.mintpy_dir, 'inputs/data.rsc')
+    rsc_file = os.path.join(os.path.dirname(inps.meta_file), 'data.rsc')
     meta = extract_gmtsar_metadata(
+        meta_file=inps.meta_file,
         unw_file=inps.unw_files[0],
         template_file=inps.template_file,
         rsc_file=rsc_file,
         update_mode=inps.update_mode,
     )
 
+    # read baseline
+    bperp_dict = read_baseline_table(inps.bperp_file)
+
     # prepare metadata for geometry files
     prepare_geometry(
-        geom_files=[inps.dem_file],
+        geom_files=[inps.dem_file, inps.inc_angle_file, inps.az_angle_file],
         meta=meta,
         update_mode=inps.update_mode,
     )
@@ -280,6 +320,7 @@ def prep_gmtsar(inps):
     prepare_stack(
         unw_files=inps.unw_files,
         meta=meta,
+        bperp_dict=bperp_dict,
         update_mode=inps.update_mode,
     )
 

src/mintpy/utils/ptime.py

@@ -236,6 +236,29 @@ def yymmdd2yyyymmdd(date):
     return date
 
 
+def yyyyddd2yyyymmdd(date_in):
+    """Convert GMTSAR folder name in YYYYDDD into YYYYMMDD.
+    Parameters: date_in  - str/list, GMTSAR date format in YYYYDDD, where DDD is the day of year - 1
+    Returns:    date_out - str/list, date in YYYYMMDD format
+    """
+    if isinstance(date_in, str):
+        year, doy = date_in[:4], date_in[-3:]
+        dt_obj = dt.datetime(int(year), 1, 1) + dt.timedelta(days=int(doy))
+        date_out = dt_obj.strftime('%Y%m%d')
+
+    elif isinstance(date_in, list):
+        date_out = []
+        for date_str in date_in:
+            year, doy = date_str[:4], date_str[-3:]
+            dt_obj = dt.datetime(int(year), 1, 1) + dt.timedelta(days=int(doy))
+            date_out.append(dt_obj.strftime('%Y%m%d'))
+
+    else:
+        return None
+
+    return date_out
+
+
 def yy2yyyy(year):
     """Convert year str from YY to YYYY format"""
     if year[0] == '9':

src/mintpy/utils/readfile.py

@@ -207,6 +207,21 @@
     'nan'   : np.nan,
 }
 
+GMTSAR_SENSOR_ID2NAME = {
+    1 : 'ers',
+    2 : 'ers',
+    3 : 'rs1',
+    4 : 'env',
+    5 : 'alos',
+    6 : 'tsx',
+    7 : 'csk',
+    8 : 'csk',
+    9 : 'rs2',
+    10: 'sen',
+    11: 'gf3',
+    12: 'lt1',
+}
+
 
 #########################################################################
 def numpy_to_gdal_dtype(np_dtype: DTypeLike) -> int:
@@ -1816,6 +1831,9 @@ def _attribute_gmtsar2roipac(prm_dict_in):
     t_center = dt_center - int(dt_center)
     prm_dict['CENTER_LINE_UTC'] = str(t_center * 24. * 60. * 60.)
 
+    # SC_identity -> PLATFORM
+    prm_dict['PLATFORM'] = GMTSAR_SENSOR_ID2NAME[int(prm_dict['SC_identity'])]
+
     return prm_dict