prep_nisar: update for changes in NISAR GUNW products

src/mintpy/cli/prep_nisar.py

@@ -41,12 +41,12 @@ def create_parser(subparsers=None):
     )
 
     parser.add_argument(
-        "-o",
-        "--out-dir",
-        dest="out_dir",
+        "-m",
+        "--mask",
+        dest="mask_file",
         type=str,
-        default="./mintpy",
-        help="output directory (default: %(default)s).",
+        default=None,
+        help="path to the mask (default: %(default)s).",
     )
 
     parser.add_argument(

src/mintpy/load_data.py

@@ -634,9 +634,13 @@ def prepare_metadata(iDict):
     elif processor == 'nisar':
         dem_file = iDict['mintpy.load.demFile']
         gunw_files = iDict['mintpy.load.unwFile']
+        water_mask = iDict['mintpy.load.waterMaskFile']
 
         # run prep_*.py
-        iargs = ['-i', gunw_files, '-d', dem_file, '-o', '../mintpy']
+        iargs = ['-i', gunw_files, '-d', dem_file]
+
+        if os.path.exists(water_mask):
+            iargs = iargs + ['--mask', water_mask]
 
         if iDict['mintpy.subset.yx']:
             warnings.warn('Subset in Y/X is not implemented for NISAR. \n'

src/mintpy/prep_nisar.py

@@ -27,8 +27,8 @@
     'unw'              : f"{DATASET_ROOT_UNW}/POL/unwrappedPhase",
     'cor'              : f"{DATASET_ROOT_UNW}/POL/coherenceMagnitude",
     'connComp'         : f"{DATASET_ROOT_UNW}/POL/connectedComponents",
-    'layoverShadowMask': f"{DATASET_ROOT_UNW}/layoverShadowMask",
-    'waterMask'        : f"{DATASET_ROOT_UNW}/waterMask",
+    'mask'             : f"{DATASET_ROOT_UNW}/mask",
+    # 'waterMask'        : f"{DATASET_ROOT_UNW}/mask",
     'epsg'             : f"{DATASET_ROOT_UNW}/projection",
     'xSpacing'         : f"{DATASET_ROOT_UNW}/xCoordinateSpacing",
     'ySpacing'         : f"{DATASET_ROOT_UNW}/yCoordinateSpacing",
@@ -61,16 +61,16 @@ def load_nisar(inps):
     print(f"Found {len(input_files)} unwrapped files")
 
     if inps.subset_lat:
-        bbox = (inps.subset_lat[0], inps.subset_lon[0], inps.subset_lat[1], inps.subset_lon[1])
+        bbox = (inps.subset_lon[0], inps.subset_lat[0], inps.subset_lon[1], inps.subset_lat[1])
     else:
-        bbox=None
+        bbox = None
 
     # extract metadata
     metadata, bounds = extract_metadata(input_files, bbox)
 
     # output filename
-    stack_file = os.path.join(inps.out_dir, "inputs/ifgramStack.h5")
-    geometry_file = os.path.join(inps.out_dir, "inputs/geometryGeo.h5")
+    stack_file = os.path.join("./inputs/ifgramStack.h5")
+    geometry_file = os.path.join("./inputs/geometryGeo.h5")
 
     # get date info: date12_list
     date12_list = _get_date_pairs(input_files)
@@ -80,7 +80,8 @@ def load_nisar(inps):
         metaFile=input_files[0],
         bbox=bounds,
         metadata=metadata,
-        demFile=inps.dem_file
+        demFile=inps.dem_file,
+        maskFile=inps.mask_file
         )
 
     prepare_stack(
@@ -136,7 +137,15 @@ def extract_metadata(input_files, bbox=None, polarization='HH'):
     meta["Y_FIRST"] = y_origin - pixel_height//2
     meta["X_STEP"] = pixel_width
     meta["Y_STEP"] = pixel_height
-    meta["X_UNIT"] = meta["Y_UNIT"] = "meters"
+
+    if meta["EPSG"] == 4326:
+        meta["X_UNIT"] = meta["Y_UNIT"] = 'degree'
+    else:
+        meta["X_UNIT"] = meta["Y_UNIT"] = "meters"
+        if str(meta["EPSG"]).startswith('326'):
+             meta["UTM_ZONE"] = str(meta["EPSG"])[3:] + 'N'
+        else:
+             meta["UTM_ZONE"] = str(meta["EPSG"])[3:] + 'S'
     meta["EARTH_RADIUS"] = EARTH_RADIUS
 
     # NISAR Altitude
@@ -155,7 +164,6 @@ def extract_metadata(input_files, bbox=None, polarization='HH'):
         utm_bbox = None
     bounds = common_raster_bound(input_files, utm_bbox)
     meta['bbox'] = ",".join([str(b) for b in bounds])
-
     col1, row1, col2, row2 = get_rows_cols(xcoord, ycoord, bounds)
     length = row2 - row1
     width = col2 - col1
@@ -225,8 +233,8 @@ def read_subset(inp_file, bbox, geometry=False):
 
         if geometry:
             # Set all values to 1 temporarily because water mask is zero
-            dataset['water_mask'] = ds[DATASETS['waterMask']][row1:row2, col1:col2] * 0 + 1
-            dataset['layover_shadow_mask'] = ds[DATASETS['layoverShadowMask']][row1:row2, col1:col2]
+            # dataset['water_mask'] = ds[DATASETS['waterMask']][row1:row2, col1:col2] * 0 + 1
+            dataset['mask'] = ds[DATASETS['mask']][row1:row2, col1:col2]
             dataset['xybbox'] = (col1, row1, col2, row2)
         else:
             dataset['unw_data'] = ds[DATASETS['unw']][row1:row2, col1:col2]
@@ -237,13 +245,16 @@ def read_subset(inp_file, bbox, geometry=False):
     return dataset
 
 
-def read_and_interpolate_geometry(gunw_file, dem_file, xybbox):
-    """Read the DEM, change projection and interpolate to data grid, interpolate slant range and incidence angle to data grid"""
-    dataset = gdal.Open(dem_file, gdal.GA_ReadOnly)
-    geotransform = dataset.GetGeoTransform()
-    proj = gdal.osr.SpatialReference(wkt=dataset.GetProjection())
-    src_epsg = int(proj.GetAttrValue('AUTHORITY', 1))
-    raster_array = dataset.ReadAsArray()
+def read_and_interpolate_geometry(gunw_file, dem_file, xybbox, mask_file=None):
+    """Read the DEM and mask, change projection and interpolate to data grid,
+       interpolate slant range and incidence angle to data grid"""
+    dem_dataset = gdal.Open(dem_file, gdal.GA_ReadOnly)
+    geotransform = dem_dataset.GetGeoTransform()
+    proj = gdal.osr.SpatialReference(wkt=dem_dataset.GetProjection())
+    dem_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1))
+    dem_raster_array = dem_dataset.ReadAsArray()
+    del dem_dataset
+
     rdr_coords = {}
 
     with h5py.File(gunw_file, 'r') as ds:
@@ -263,20 +274,47 @@ def read_and_interpolate_geometry(gunw_file, dem_file, xybbox):
     # dem_subset_array = np.zeros((subset_rows, subset_cols), dtype=raster_array.dtype)
     Y_2d, X_2d = np.meshgrid(ycoord, xcoord, indexing='ij')
 
-    if not src_epsg == dst_epsg:
-        coord_transform = Transformer.from_crs(dst_epsg, src_epsg, always_xy=True)
+    if not dem_src_epsg == dst_epsg:
+        coord_transform = Transformer.from_crs(dst_epsg, dem_src_epsg, always_xy=True)
         x_dem, y_dem = coord_transform.transform(X_2d.flatten(), Y_2d.flatten())
     else:
         x_dem, y_dem = X_2d.flatten(), Y_2d.flatten()
 
     cols = ((y_dem - geotransform[3]) / geotransform[5]).astype(int)
     rows = ((x_dem - geotransform[0]) / geotransform[1]).astype(int)
 
-    dem_subset_array = raster_array[cols.reshape(subset_rows, subset_cols), rows.reshape(subset_rows, subset_cols)]
+    dem_subset_array = dem_raster_array[cols.reshape(subset_rows, subset_cols), rows.reshape(subset_rows, subset_cols)]
 
     slant_range, incidence_angle = interpolate_geometry(X_2d, Y_2d, dem_subset_array, rdr_coords)
 
-    return dem_subset_array, slant_range, incidence_angle
+
+    if mask_file in ['auto', 'None', None]:
+        print('*** No mask was found ***')
+        mask_subset_array = np.ones(dem_subset_array.shape, dtype='byte')
+    else:
+        try:
+            mask_dataset = gdal.Open(mask_file, gdal.GA_ReadOnly)
+            geotransform = mask_dataset.GetGeoTransform()
+            proj = gdal.osr.SpatialReference(wkt=mask_dataset.GetProjection())
+            mask_src_epsg = int(proj.GetAttrValue('AUTHORITY', 1))
+            mask_raster_array = mask_dataset.ReadAsArray()
+            del mask_dataset
+
+            if not mask_src_epsg == dst_epsg:
+                coord_transform = Transformer.from_crs(dst_epsg, mask_src_epsg, always_xy=True)
+                x_mask, y_mask = coord_transform.transform(X_2d.flatten(), Y_2d.flatten())
+            else:
+                x_mask, y_mask = X_2d.flatten(), Y_2d.flatten()
+
+            cols = ((y_mask - geotransform[3]) / geotransform[5]).astype(int)
+            rows = ((x_mask - geotransform[0]) / geotransform[1]).astype(int)
+            mask_subset_array = mask_raster_array[cols.reshape(subset_rows, subset_cols),
+                                                rows.reshape(subset_rows, subset_cols)]
+        except:
+            raise OSError('*** Mask is not gdal readable ***')
+
+
+    return dem_subset_array, slant_range, incidence_angle, mask_subset_array
 
 
 def interpolate_geometry(X_2d, Y_2d, dem, rdr_coords):
@@ -315,7 +353,8 @@ def prepare_geometry(
         metaFile,
         metadata,
         bbox,
-        demFile
+        demFile,
+        maskFile
 ):
     """Prepare the geometry file."""
     print("-" * 50)
@@ -326,17 +365,20 @@ def prepare_geometry(
 
     # Read waterMask, LayoverShadowMask, xybbox:
     geo_ds = read_subset(metaFile, bbox, geometry=True)
-    dem_subset_array, slant_range, incidence_angle = read_and_interpolate_geometry(metaFile, demFile, geo_ds['xybbox'])
+    dem_subset_array, slant_range, incidence_angle, mask = read_and_interpolate_geometry(metaFile, demFile,
+                                                                                   geo_ds['xybbox'], mask_file=maskFile)
 
     length, width = dem_subset_array.shape
 
     ds_name_dict = {
         "height": [np.float32, (length, width), dem_subset_array],
         "incidenceAngle": [np.float32, (length, width), incidence_angle],
         "slantRangeDistance": [np.float32, (length, width), slant_range],
-        "shadowMask": [np.bool_, (length, width), geo_ds['layover_shadow_mask']],
-        "waterMask": [np.bool_, (length, width), geo_ds['water_mask']],
+        "shadowMask": [np.bool_, (length, width), geo_ds['mask']],
+        #"waterMask": [np.bool_, (length, width), geo_ds['water_mask']],
     }
+    if maskFile:
+        ds_name_dict['waterMask'] = [np.bool_, (length, width), mask]
 
     # initiate HDF5 file
     meta["FILE_TYPE"] = "geometry"

src/mintpy/tropo_pyaps3.py

@@ -419,7 +419,13 @@ def dload_grib_files(grib_files, tropo_model='ERA5', snwe=None):
 
     # Download grib file using PyAPS
     if len(date_list2dload) > 0:
-        hour = re.findall(r'\d{8}[-_]\d{2}', os.path.basename(grib_files2dload[0]))[0].replace('-', '_').split('_')[1]
+        # This was the default for a file pattern of "ERA5_N30_N50_W130_W110_20200403_14.grb"
+        pattern = re.findall(r'\d{8}[-_]\d{2}', os.path.basename(grib_files2dload[0]))
+        if len(pattern) == 0:
+            # This is an exception for when the file has a naming format like: "ERA5_N30_N40_W130_W110_20080125T000000_06.grb"
+            pattern = re.findall(r'\d{8}T\d{6}[-_]\d{2}', os.path.basename(grib_files2dload[0]))
+        hour = pattern[0].replace('-', '_').split('_')[1]
+        #hour = re.findall(r'\d{8}[-_]\d{2}', os.path.basename(grib_files2dload[0]))[0].replace('-', '_').split('_')[1]
         grib_dir = os.path.dirname(grib_files2dload[0])
 
         # Check for non-empty account info in PyAPS config file