Merge pull request #770 from jungkyoj/layover_shadow_mask

Add Layover/shadow mask in InSAR workflow

Author: jungkyoJung

python/extensions/pybind_isce3/geometry/rdr2geo.cpp

@@ -178,7 +178,7 @@ void addbinding(py::class_<Topo> & pyRdr2Geo)
             py::arg("extraiter") = 10,
             py::arg("dem_interp_method") = isce3::core::BIQUINTIC_METHOD,
             py::arg("epsg_out") = 4326,
-            py::arg("compute_mask") = false)
+            py::arg("compute_mask") = true)
         .def("topo", py::overload_cast<isce3::io::Raster &, const std::string &>
                 (&Topo::topo),
                 py::arg("dem_raster"),

python/packages/pybind_nisar/workflows/geocode_insar.py

@@ -8,6 +8,7 @@
 
 import h5py
 import journal
+import pathlib
 import pybind_isce3 as isce3
 from pybind_nisar.products.readers import SLC
 from pybind_nisar.workflows import h5_prep
@@ -34,6 +35,7 @@ def run(cfg, runw_hdf5, output_hdf5):
     rg_looks = cfg["processing"]["crossmul"]["range_looks"]
     interp_method = cfg["processing"]["geocode"]["interp_method"]
     gunw_datasets = cfg["processing"]["geocode"]["datasets"]
+    scratch_path = pathlib.Path(cfg['ProductPathGroup']['ScratchPath'])
 
     slc = SLC(hdf5file=ref_hdf5)
 
@@ -67,9 +69,9 @@ def run(cfg, runw_hdf5, output_hdf5):
         for freq in freq_pols.keys():
             pol_list = freq_pols[freq]
 
-            radar_grid = slc.getRadarGrid(freq)
+            radar_grid_slc = slc.getRadarGrid(freq)
             if az_looks > 1 or rg_looks > 1:
-                radar_grid = radar_grid.multilook(az_looks, rg_looks)
+                radar_grid_multilook = radar_grid_slc.multilook(az_looks, rg_looks)
             geo_grid = geogrids[freq]
             geo.geogrid(
                 geo_grid.start_x,
@@ -80,24 +82,39 @@ def run(cfg, runw_hdf5, output_hdf5):
                 geo_grid.length,
                 geo_grid.epsg,
             )
+            src_freq_path = f"/science/LSAR/RUNW/swaths/frequency{freq}/interferogram"
+            dst_freq_path = f"/science/LSAR/GUNW/grids/frequency{freq}/interferogram"
 
             for pol in pol_list:
-                src_group_path = f"/science/LSAR/RUNW/swaths/frequency{freq}/interferogram/{pol}"
-                dst_group_path = f"/science/LSAR/GUNW/grids/frequency{freq}/interferogram/{pol}"
+                src_group_path = f"{src_freq_path}/{pol}"
+                dst_group_path = f"{dst_freq_path}/{pol}"
 
                 # iterate over key: dataset name value: bool flag to perform geocode
                 for dataset_name, geocode_this_dataset in gunw_datasets.items():
                     if not geocode_this_dataset:
                         continue
 
-                    # prepare input raster
-                    input_raster_str = (
-                        f"HDF5:{runw_hdf5}:/{src_group_path}/{dataset_name}"
-                    )
-                    input_raster = isce3.io.Raster(input_raster_str)
+                    if (dataset_name == "layoverShadowMask"):
+                        # prepare input raster
+                        raster_ref = scratch_path / 'rdr2geo' / f'freq{freq}' / 'mask.rdr'
+                        input_raster = isce3.io.Raster(str(raster_ref))
+                       
+                        # access the HDF5 dataset for layover shadow mask
+                        dataset_path = f"{dst_freq_path}/{dataset_name}"
+                        geo.data_interpolator = 'NEAREST'
+                        radar_grid = radar_grid_slc
+                    else:  
+                        # prepare input raster
+                        input_raster_str = (
+                            f"HDF5:{runw_hdf5}:/{src_group_path}/{dataset_name}"
+                        )
+                        input_raster = isce3.io.Raster(input_raster_str)
+
+                        # access the HDF5 dataset for a given frequency and pol
+                        geo.data_interpolator = interp_method
+                        radar_grid = radar_grid_multilook
+                        dataset_path = f"{dst_group_path}/{dataset_name}"
 
-                    # access the HDF5 dataset for a given frequency and pol
-                    dataset_path = f"{dst_group_path}/{dataset_name}"
                     geocoded_dataset = dst_h5[dataset_path]
 
                     # Construct the output ratster directly from HDF5 dataset

python/packages/pybind_nisar/workflows/geocode_insar_runconfig.py

@@ -46,7 +46,7 @@ def yaml_check(self):
             self.cfg['processing']['geocode']['datasets'] = {}
 
         # default to True for datasets not found
-        gunw_datasets = ["connectedComponents", "coherenceMagnitude", "unwrappedPhase"]
+        gunw_datasets = ["connectedComponents", "coherenceMagnitude", "unwrappedPhase", "layoverShadowMask"]
         for gunw_dataset in gunw_datasets:
             if gunw_dataset not in self.cfg['processing']['geocode']:
                 self.cfg['processing']['geocode']['datasets'][gunw_dataset] = True

python/packages/pybind_nisar/workflows/h5_prep.py

@@ -669,6 +669,14 @@ def prep_ds_insar(cfg, dst, dst_h5):
                              'quality',
                              descr=descr, units=" ", long_name='quality')
 
+        # Adding layover-shadow mask
+        if dst in ['GUNW']:
+            descr = f"Layover Shadow mask for frequency{freq} layer, 1 - Radar Shadow. 2 - Radar Layover. 3 - Both"
+            _create_datasets(dst_h5[dst_path_intf], shape, np.byte,
+                             'layoverShadowMask',
+                             descr=descr, units=" ", grids=grids_val,
+                             long_name='layover shadow mask')
+
         # Add datasets in metadata
         dst_cal = os.path.join(common_parent_path,
                                f'{dst}/metadata/calibrationInformation')

python/packages/pybind_nisar/workflows/insar_runconfig.py

@@ -55,7 +55,7 @@ def yaml_check(self):
 
         # default to True for datasets not found
         gunw_datasets = ["connectedComponents", "coherenceMagnitude",
-                         "unwrappedPhase"]
+                         "unwrappedPhase", 'layoverShadowMask']
         for gunw_dataset in gunw_datasets:
             if gunw_dataset not in self.cfg['processing']['geocode']:
                 self.cfg['processing']['geocode']['datasets'][