Geocode dense offsets (#812)

* Adddefault margin in default

* Modify insar runconfig to include checks on dense offsets cfg

* Modify radar grid for geocoding dense offsets

* Add checks to dense offsets cfg for stand-alone usage of geocode_insar.py

* Correct bug for stand-alone usage of geocode_insar.py

* Correct mispelled variable in dense_offsets.py

* Add dataset to geocode in insar schema

* isort input and pep8 notation formatting

* Correct ampcor parameter in unit test

* Correct path to layover/shadow dataset

* Correct path to dataset to geocode

* Correct dense offsets search windo for unit test

* Correct bg in margin definition for geocode_insar.py

* Creating dense offsets radar grid from scratch

* Pep8 notation formatting

* Remove radar grid multilook and correct tabs in geocode_insar_runconfig.py

* Correctly allocate the radar grid for interferogram/coherence

* Correcting bug in geocode_insar_runconfig

* Correct bug in insar_runconfig.py

* Update start pixel of offset radar grid

* Remove multiplication for 1.0 in offset prf computation

Co-authored-by: vbrancat <[email protected]>

Author: 2 people

python/packages/pybind_nisar/workflows/dense_offsets.py

@@ -239,7 +239,7 @@ def set_optional_attributes(ampcor_obj, cfg, length, width):
             err_str = "The input gross offset does not match the offset width*offset length"
             error_channel.log(err_str)
             raise RuntimeError(err_str)
-        gross_offset = gross_offset.reshape(window_number, 2)
+        gross_offset = gross_offset.reshape(windows_number, 2)
         gross_azimuth = gross_offset[:, 0]
         gross_range = gross_offset[:, 1]
         ampcor_obj.setVaryingGrossOffset(gross_azimuth, gross_range)

python/packages/pybind_nisar/workflows/geocode_insar.py

@@ -4,17 +4,17 @@
 collection of functions for NISAR geocode workflow
 """
 
+import pathlib
 import time
 
 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
-from pybind_nisar.workflows.h5_prep import add_radar_grid_cubes_to_hdf5
 from pybind_nisar.workflows.geocode_insar_runconfig import \
     GeocodeInsarRunConfig
+from pybind_nisar.workflows.h5_prep import add_radar_grid_cubes_to_hdf5
 from pybind_nisar.workflows.yaml_argparse import YamlArgparse
 
 
@@ -39,6 +39,7 @@ def run(cfg, runw_hdf5, output_hdf5):
     interp_method = cfg["processing"]["geocode"]["interp_method"]
     gunw_datasets = cfg["processing"]["geocode"]["datasets"]
     scratch_path = pathlib.Path(cfg['ProductPathGroup']['ScratchPath'])
+    offset_cfg = cfg["processing"]["dense_offsets"]
 
     slc = SLC(hdf5file=ref_hdf5)
 
@@ -75,7 +76,8 @@ def run(cfg, runw_hdf5, output_hdf5):
 
             radar_grid_slc = slc.getRadarGrid(freq)
             if az_looks > 1 or rg_looks > 1:
-                radar_grid_multilook = radar_grid_slc.multilook(az_looks, rg_looks)
+               radar_grid_mlook = radar_grid_slc.multilook(az_looks, rg_looks)
+
             geo_grid = geogrids[freq]
             geo.geogrid(
                 geo_grid.start_x,
@@ -86,34 +88,90 @@ 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"
+            src_freq_path = f"/science/LSAR/RUNW/swaths/frequency{freq}"
+            dst_freq_path = f"/science/LSAR/GUNW/grids/frequency{freq}"
 
             for pol in pol_list:
-                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
 
-                    if dataset_name in ['alongTrackOffset','slantRangeOffset']:
-                        src_group_path = f"/science/LSAR/RUNW/swaths/frequency{freq}/pixelOffsets/{pol}"
-                        dst_group_path = f"/science/LSAR/GUNW/grids/frequency{freq}/pixelOffsets/{pol}"
+                    # Create radar grid for the offsets (and dataset path)
+                    if dataset_name in ['alongTrackOffset', 'slantRangeOffset']:
+                        src_group_path = f'{src_freq_path}/pixelOffsets/{pol}'
+                        dst_group_path = f'{dst_freq_path}/pixelOffsets/{pol}'
+
+                        # Define margin used during dense offsets execution
+                        margin = max(offset_cfg['margin'],
+                                     offset_cfg['gross_offset_range'],
+                                     offset_cfg['gross_offset_azimuth'])
+
+                        # If not allocated, determine shape of the offsets
+                        if offset_cfg['offset_length'] is None:
+                            length_margin = 2 * margin + 2 * offset_cfg[
+                                'half_search_azimuth'] + \
+                                            offset_cfg['window_azimuth']
+                            offset_cfg['offset_length'] = (radar_grid_slc.length -
+                                                           length_margin) // offset_cfg['skip_azimuth']
+                        if offset_cfg['offset_width'] is None:
+                            width_margin = 2 * margin + 2 * offset_cfg[
+                                'half_search_range'] + \
+                                           offset_cfg['window_range']
+                            offset_cfg['offset_width'] = (radar_grid_slc.width -
+                                                          width_margin) // offset_cfg['skip_azimuth']
+                        # Determine the starting range and sensing start for the offset radar grid
+                        offset_starting_range = radar_grid_slc.starting_range + \
+                                                (offset_cfg['start_pixel_range'] + offset_cfg['window_range']//2)\
+                                                * radar_grid_slc.range_pixel_spacing
+                        offset_sensing_start = radar_grid_slc.sensing_start + \
+                                               (offset_cfg['start_pixel_azimuth'] + offset_cfg['window_azimuth']//2)\
+                                               / radar_grid_slc.prf
+                        # Range spacing for offsets
+                        offset_range_spacing = radar_grid_slc.range_pixel_spacing * offset_cfg['skip_range']
+                        offset_prf = radar_grid_slc.prf / offset_cfg['skip_azimuth']
+
+                        # Create offset radar grid
+                        radar_grid = isce3.product.RadarGridParameters(offset_sensing_start,
+                                                                       radar_grid_slc.wavelength,
+                                                                       offset_prf,
+                                                                       offset_starting_range,
+                                                                       offset_range_spacing,
+                                                                       radar_grid_slc.lookside,
+                                                                       offset_cfg['offset_length'],
+                                                                       offset_cfg['offset_width'],
+                                                                       radar_grid_slc.ref_epoch)
+                        # 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
+                        dataset_path = f"{dst_group_path}/{dataset_name}"
 
                     # prepare input raster
-
-                    if (dataset_name == "layoverShadowMask"):
+                    elif (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}"
+                        dataset_path = f"{dst_freq_path}/interferogram/{dataset_name}"
                         geo.data_interpolator = 'NEAREST'
                         radar_grid = radar_grid_slc
-                    else:  
+                    else:
+                        # Assign correct radar grid
+                        if az_looks > 1 or rg_looks > 1:
+                            radar_grid = radar_grid_mlook
+                        else:
+                            radar_grid = radar_grid_slc
+
+                        # Prepare input path
+                        src_group_path = f'{src_freq_path}/interferogram/{pol}'
+                        dst_group_path = f'{dst_freq_path}/interferogram/{pol}'
+
                         # prepare input raster
                         input_raster_str = (
                             f"HDF5:{runw_hdf5}:/{src_group_path}/{dataset_name}"
@@ -122,10 +180,8 @@ def run(cfg, runw_hdf5, output_hdf5):
 
                         # 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}"
 
-
                     geocoded_dataset = dst_h5[dataset_path]
 
                     # Construct the output ratster directly from HDF5 dataset
@@ -165,17 +221,16 @@ def run(cfg, runw_hdf5, output_hdf5):
             computer cubes values outside radar-grid boundaries
             '''
             native_doppler.bounds_error = False
-            add_radar_grid_cubes_to_hdf5(dst_h5, cube_group_name, 
-                                         cube_geogrid, radar_grid_cubes_heights, 
-                                         radar_grid, orbit, native_doppler, 
-                                         grid_zero_doppler, threshold_geo2rdr, 
+            add_radar_grid_cubes_to_hdf5(dst_h5, cube_group_name,
+                                         cube_geogrid, radar_grid_cubes_heights,
+                                         radar_grid, orbit, native_doppler,
+                                         grid_zero_doppler, threshold_geo2rdr,
                                          iteration_geo2rdr)
 
     t_all_elapsed = time.time() - t_all
     info_channel.log(f"Successfully ran geocode in {t_all_elapsed:.3f} seconds")
 
 
-
 if __name__ == "__main__":
     """
     run geocode from command line
@@ -190,10 +245,10 @@ def run(cfg, runw_hdf5, output_hdf5):
 
     # prepare RIFG HDF5
     out_paths = h5_prep.run(geocode_insar_runconfig.cfg)
-    if args.run_config_path is None:
-        out_paths['RUNW'] = args.runw_h5
-    else:
-        out_paths['RUNW'] = geocode_insar_runconfig.cfg['processing']['geocode']['runw_path']
+    runw_path = geocode_insar_runconfig.cfg['processing']['geocode'][
+        'runw_path']
+    if runw_path is not None:
+        out_paths['RUNW'] = runw_path
 
     # Run geocode
     run(geocode_insar_runconfig.cfg, out_paths["RUNW"], out_paths["GUNW"])

python/packages/pybind_nisar/workflows/geocode_insar_runconfig.py

@@ -49,7 +49,7 @@ def yaml_check(self):
         gunw_datasets = ["connectedComponents", "coherenceMagnitude", "unwrappedPhase",
                          "alongTrackOffset", "slantRangeOffset", "layoverShadowMask"]
         for gunw_dataset in gunw_datasets:
-            if gunw_dataset not in self.cfg['processing']['geocode']:
+            if gunw_dataset not in self.cfg['processing']['geocode']['datasets']:
                 self.cfg['processing']['geocode']['datasets'][gunw_dataset] = True
 
         if self.cfg['processing']['dem_margin'] is None:
@@ -74,3 +74,22 @@ def yaml_check(self):
             err_str = f"range looks = {rg_looks} not an odd integer."
             error_channel.log(err_str)
             raise ValueError(err_str)
+
+        # To geocode the offsets we need the offset field shape and
+        # the start pixel in range and azimuth. Note, margin and gross_offsets
+        # are allocated as defaults in share/nisar/defaults/insar.yaml
+        geocode_azimuth_offset = self.cfg['processing']['geocode']['datasets'][
+                'alongTrackOffset']
+        geocode_range_offset = self.cfg['processing']['geocode']['datasets'][
+                'slantRangeOffset']
+        if geocode_azimuth_offset or geocode_range_offset:
+            offset_cfg = self.cfg['processing']['dense_offsets']
+            margin = max(offset_cfg['margin'],
+                         offset_cfg['gross_offset_range'],
+                         offset_cfg['gross_offset_azimuth'])
+            if offset_cfg['start_pixel_range'] is None:
+               offset_cfg['start_pixel_range'] = margin + offset_cfg[
+                          'half_search_range']
+            if offset_cfg['start_pixel_azimuth'] is None:
+               offset_cfg['start_pixel_azimuth'] = margin + offset_cfg[
+                          'half_search_azimuth']

python/packages/pybind_nisar/workflows/insar_runconfig.py

@@ -1,3 +1,5 @@
+import journal
+import numpy as np
 from pybind_nisar.workflows.geo2rdr_runconfig import Geo2rdrRunConfig
 
 
@@ -13,6 +15,7 @@ def yaml_check(self):
         Check submodule paths from YAML
         '''
         scratch_path = self.cfg['ProductPathGroup']['ScratchPath']
+        error_channel = journal.error('InsarRunConfig.yaml_check')
 
         # If dense_offsets is disabled and rubbersheet is enabled
         # throw an exception and do not run the workflow
@@ -73,6 +76,24 @@ def yaml_check(self):
                          'layoverShadowMask']
 
         for gunw_dataset in gunw_datasets:
-            if gunw_dataset not in self.cfg['processing']['geocode']:
+            if gunw_dataset not in self.cfg['processing']['geocode']['datasets']:
                 self.cfg['processing']['geocode']['datasets'][
                     gunw_dataset] = True
+
+        # To geocode the offsets we need the offset field shape and
+        # the start pixel in range and azimuth. Note, margin and gross_offsets
+        # are allocated as defaults in share/nisar/defaults/insar.yaml
+        geocode_azimuth_offset = self.cfg['processing'][
+            'geocode']['datasets']['alongTrackOffset']
+        geocode_range_offset = self.cfg['processing'][
+            'geocode']['datasets']['slantRangeOffset']
+        if geocode_azimuth_offset or geocode_range_offset:
+            offset_cfg = self.cfg['processing']['dense_offsets']
+            margin = max(offset_cfg['margin'], offset_cfg['gross_offset_range'],
+                         offset_cfg['gross_offset_azimuth'])
+            if offset_cfg['start_pixel_range'] is None:
+                offset_cfg['start_pixel_range'] = margin + offset_cfg[
+                    'half_search_range']
+            if offset_cfg['start_pixel_azimuth'] is None:
+                offset_cfg['start_pixel_azimuth'] = margin + offset_cfg[
+                    'half_search_azimuth']

share/nisar/defaults/insar.yaml

@@ -171,7 +171,7 @@ runconfig:
                 # Number of lines to skip in reference image for next offset estimate
                 skip_azimuth: 32
                 # Margin around image edges to avoid for offset computation
-                margin:
+                margin: 0
                 # Number of columns of output offset (covariance, snr) file
                 offset_width:
                 # Number of lines of output offset (covariance, snr) file

share/nisar/schemas/insar.yaml

@@ -405,7 +405,9 @@ gunw_datasets:
     connectedComponents: bool(required=False)
     coherenceMagnitude: bool(required=False)
     unwrappedPhase: bool(required=False)
-
+    layoverShadowMask: bool(required=False)
+    alongTrackOffset: bool(required=False)
+    slantRangeOffset: bool(required=False)
 
 qa_options:
     # Enabled file format validation

tests/data/insar_test.yaml

@@ -117,8 +117,8 @@ runconfig:
                 enabled: True
                 window_range: 64
                 window_azimuth: 64
-                half_search_range: 20
-                half_search_azimuth: 20
+                half_search_range: 10
+                half_search_azimuth: 10
                 skip_range: 8
                 skip_azimuth: 8
                 correlation_surface_zoom: 8