Add PHASS unwrapper to InSAR workflow (#777)

* Update defaults and schemas to include both phass and icu unwrappers

* Include unwrapper choice in InSAR workflow

* Include unwrapper preferences in unwrap workflow

* Allocate temporary wrapped phase raster

* Using gdal pixel function to compute phase raster

* Correct rifg variable name to rifg_h5

* Specifying that crossmul_path is required for stand-alone phase unwrapping usage

* Resolve issue with merging develop (misplaced connected component)

* Remove extra space, rename unwrap write attribute function, move code out of if statement

* Improving comments on good_correlation

* Update description of coherence in insar schema/defaults

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

Author: 2 people

python/packages/pybind_nisar/workflows/insar_runconfig.py

@@ -55,14 +55,10 @@ def yaml_check(self):
         else:
             self.cfg['processing']['crossmul']['flatten'] = None
 
-        # set to empty dict and default unwrap values will be used
-        # if phase_unwrap fields not in yaml
-        if 'phase_unwrap' not in self.cfg['processing']:
-            self.cfg['processing']['phase_unwrap'] = {}
-
-        # if phase_unwrap fields not in yaml
-        if self.cfg['processing']['phase_unwrap'] is None:
-            self.cfg['processing']['phase_unwrap'] = {}
+        # Create default unwrap cfg dict depending on unwrapping algorithm
+        algorithm = self.cfg['processing']['phase_unwrap']['algorithm']
+        if algorithm not in self.cfg['processing']['phase_unwrap']:
+            self.cfg['processing']['phase_unwrap'][algorithm]={}
 
         if 'interp_method' not in self.cfg['processing']['geocode']:
             self.cfg['processing']['geocode']['interp_method'] = 'BILINEAR'

python/packages/pybind_nisar/workflows/unwrap.py

@@ -9,7 +9,10 @@
 
 import h5py
 import journal
+import os
 import pybind_isce3 as isce3
+from osgeo import gdal
+import numpy as np
 from pybind_nisar.workflows import h5_prep
 from pybind_nisar.workflows.unwrap_runconfig import UnwrapRunConfig
 from pybind_nisar.workflows.yaml_argparse import YamlArgparse
@@ -22,6 +25,7 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
 
     # pull parameters from dictionary
     freq_pols = cfg['processing']['input_subset']['list_of_frequencies']
+    scratch_path = pathlib.Path(cfg['ProductPathGroup']['ScratchPath'])
     unwrap_args = cfg['processing']['phase_unwrap']
 
     # Create error and info channels
@@ -36,11 +40,21 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
         error_channel.log(err_str)
         raise ValueError(err_str)
 
-    # ICU unwrapping (CPU only)
-    unwrap_obj = isce3.unwrap.ICU()
+    # Instantiate correct unwrap object
+    if unwrap_args['algorithm'] == 'icu':
+        unwrap_obj = isce3.unwrap.ICU()
+    elif unwrap_args['algorithm'] == 'phass':
+        unwrap_obj = isce3.unwrap.Phass()
+    else:
+        err_str = "Invalid unwrapping algorithm"
+        error_channel.log(err_str)
+        raise ValueError(err_str)
+
+    # Depending on unwrapper, set unwrap attributes
     unwrap_obj = set_unwrap_attributes(unwrap_obj, unwrap_args)
     info_channel.log("Running phase unwrapping with: ")
-    info_channel = write_unwrap_attributes(unwrap_obj, info_channel)
+    log_unwrap_attributes(unwrap_obj, info_channel,
+                          unwrap_args['algorithm'])
 
     t_all = time.time()
 
@@ -54,10 +68,9 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
                 src_pol_group_path = f'{src_freq_group_path}/interferogram/{pol}'
                 dst_pol_group_path = f'{dst_freq_group_path}/interferogram/{pol}'
 
-                # Prepare igram input raster
+                # Interferogram filepath
                 igram_path = f'HDF5:{crossmul_path}:/' \
                              f'{src_pol_group_path}/wrappedInterferogram'
-                igram_raster = isce3.io.Raster(igram_path)
 
                 # Prepare correlation input raster
                 corr_path = f'HDF5:{crossmul_path}:/' \
@@ -75,6 +88,54 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
                 conn_comp_dataset = dst_h5[conn_comp_path]
                 conn_comp_raster = isce3.io.Raster(f"IH5:::ID={conn_comp_dataset.id.id}".encode("utf-8"),
                                                    update=True)
+
+                # If unwrapping algorithm is ICU, run it with or without seed
+                if unwrap_args['algorithm'] == 'icu':
+                    # Allocate interferogram as ISCE3 raster
+                    igram_raster = isce3.io.Raster(igram_path)
+                    if 'seed' in unwrap_args['icu']:
+                        unwrap_obj.unwrap(uigram_raster, conn_comp_raster,
+                                          igram_raster,
+                                          corr_raster, seed=unwrap_args['seed'])
+                    else:
+                        unwrap_obj.unwrap(uigram_raster, conn_comp_raster,
+                                          igram_raster,
+                                          corr_raster)
+                else:
+                    # Unwrapping algorithm is PHASS which requires
+                    # the interferometric phase as input raster
+                    unwrap_scratch = scratch_path / f'unwrap/freq{freq}/{pol}'
+                    unwrap_scratch.mkdir(parents=True, exist_ok=True)
+
+                    # Using GDAL pixel function to compute a wrapped phase VRT
+                    ds = gdal.Open(igram_path, gdal.GA_ReadOnly)
+                    vrttmpl = f'''
+                            <VRTDataset rasterXSize="{ds.RasterXSize}" rasterYSize="{ds.RasterYSize}">
+                            <VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
+                            <Description>Phase</Description>
+                            <PixelFunctionType>phase</PixelFunctionType>
+                            <SimpleSource>
+                            <SourceFilename>{igram_path}</SourceFilename>
+                            </SimpleSource>
+                            </VRTRasterBand>
+                            </VRTDataset>'''
+                    ds = None
+                    with open(os.path.join(unwrap_scratch, 'wrapped_phase.vrt'),
+                              'w') as fid:
+                        fid.write(vrttmpl)
+
+                    # Open phase_raster as ISCE3 raster
+                    phase_raster = isce3.io.Raster(os.path.join(unwrap_scratch, 'wrapped_phase.vrt'))
+
+                    # Check if power raster has been allocated
+                    if 'power' in unwrap_args['phass']:
+                        power_raster = isce3.io.Raster(unwrap_args['phass']['power'])
+                        unwrap_obj.unwrap(phase_raster, power_raster,
+                                          corr_raster, uigram_raster, conn_comp_raster)
+                    else:
+                        unwrap_obj.unwrap(phase_raster, corr_raster,
+                                          uigram_raster, conn_comp_raster)
+
                 if 'seed' in unwrap_args:
                     unwrap_obj.unwrap(uigram_raster, conn_comp_raster, igram_raster,
                                       corr_raster, seed=unwrap_args['seed'])
@@ -85,8 +146,7 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
                 # Copy coherence magnitude and culled offsets to RIFG
                 dataset_names = ['coherenceMagnitude', 'alongTrackOffset',
                                  'slantRangeOffset']
-                group_names =['interferogram', 'pixelOffsets',
-                              'pixelOffsets']
+                group_names = ['interferogram', 'pixelOffsets', 'pixelOffsets']
                 for dataset_name, group_name in zip(dataset_names, group_names):
                     dst_path = f'{dst_freq_group_path}/{group_name}/{pol}/{dataset_name}'
                     src_path = f'{src_freq_group_path}/{group_name}/{pol}/{dataset_name}'
@@ -99,29 +159,36 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
     info_channel.log(f"Successfully ran phase unwrapping in {t_all_elapsed:.3f} seconds")
 
 
-def write_unwrap_attributes(unwrap, info):
+def log_unwrap_attributes(unwrap, info, algorithm):
     '''
-    Write unwrap attributes to info channel
+    Write unwrap attributes to info
+    channel depending on unwrapping algorithm
     '''
-    info.log(f"Number of buffer lines: {unwrap.buffer_lines}")
-    info.log(f"Number of overlap lines: {unwrap.overlap_lines}")
-    info.log(f"Use phase gradient neutron: {unwrap.use_phase_grad_neut}")
-    info.log(f"Use intensity neutron: {unwrap.use_intensity_neut}")
-    info.log(f"Phase gradient window size: {unwrap.phase_grad_win_size}")
-    info.log(f"Phase gradient threshold: {unwrap.neut_phase_grad_thr}")
-    info.log(f"Neutron intensity threshold: {unwrap.neut_intensity_thr}")
-    info.log(
-        f"Maximum intesity correlation threshold: {unwrap.neut_correlation_thr}")
-    info.log(f"Number of trees: {unwrap.trees_number}")
-    info.log(f"Maximum branch length: {unwrap.max_branch_length}")
-    info.log(f"Pixel spacing ratio: {unwrap.ratio_dxdy}")
-    info.log(f"Initial correlation threshold: {unwrap.init_corr_thr}")
-    info.log(f"Maximum correlation threshold: {unwrap.max_corr_thr}")
-    info.log(f"Correlation threshold increments: {unwrap.corr_incr_thr}")
-    info.log(f"Minimum tile area fraction: {unwrap.min_cc_area}")
-    info.log(f"Number of bootstraping lines: {unwrap.num_bs_lines}")
-    info.log(f"Minimum overlapping area: {unwrap.min_overlap_area}")
-    info.log(f"Phase variance threshold: {unwrap.phase_var_thr}")
+    info.log(f"Unwrapping algorithm:{algorithm}")
+    if algorithm == 'icu':
+        info.log(f"Number of buffer lines: {unwrap.buffer_lines}")
+        info.log(f"Number of overlap lines: {unwrap.overlap_lines}")
+        info.log(f"Use phase gradient neutron: {unwrap.use_phase_grad_neut}")
+        info.log(f"Use intensity neutron: {unwrap.use_intensity_neut}")
+        info.log(f"Phase gradient window size: {unwrap.phase_grad_win_size}")
+        info.log(f"Phase gradient threshold: {unwrap.neut_phase_grad_thr}")
+        info.log(f"Neutron intensity threshold: {unwrap.neut_intensity_thr}")
+        info.log(f"Maximum intensity correlation "
+                 f"threshold: {unwrap.neut_correlation_thr}")
+        info.log(f"Number of trees: {unwrap.trees_number}")
+        info.log(f"Maximum branch length: {unwrap.max_branch_length}")
+        info.log(f"Pixel spacing ratio: {unwrap.ratio_dxdy}")
+        info.log(f"Initial correlation threshold: {unwrap.init_corr_thr}")
+        info.log(f"Maximum correlation threshold: {unwrap.max_corr_thr}")
+        info.log(f"Correlation threshold increments: {unwrap.corr_incr_thr}")
+        info.log(f"Minimum tile area fraction: {unwrap.min_cc_area}")
+        info.log(f"Number of bootstraping lines: {unwrap.num_bs_lines}")
+        info.log(f"Minimum overlapping area: {unwrap.min_overlap_area}")
+        info.log(f"Phase variance threshold: {unwrap.phase_var_thr}")
+    else:
+        info.log(f"Good correlation: {unwrap.good_correlation}")
+        info.log(
+            f"Minimum size of an unwrapped region: {unwrap.min_pixels_region}")
 
     return info
 
@@ -131,43 +198,60 @@ def set_unwrap_attributes(unwrap, cfg: dict):
     Assign user-defined values in cfg to
     the unwrap object
     '''
+    error_channel = journal.error('unwrap.set_unwrap_attributes')
 
-    if 'buffer_lines' in cfg:
-        unwrap.buffer_lines = cfg['buffer_lines']
-    if 'overlap_lines' in cfg:
-        unwrap.overlap_lines = cfg['overlap_lines']
-    if 'use_phase_gradient_neutron' in cfg:
-        unwrap.use_phase_grad_neut = cfg['use_phase_gradient_neutron']
-    if 'use_intensity_neutron' in cfg:
-        unwrap.use_intensity_neut = cfg['use_intensity_neutron']
-    if 'phase_gradient_window_size' in cfg:
-        unwrap.phase_grad_win_size = cfg['phase_gradient_window_size']
-    if 'neutron_phase_gradient_threshold' in cfg:
-        unwrap.neut_phase_grad_thr = cfg['neutron_phase_gradient_threshold']
-    if 'neutron_intensity_threshold' in cfg:
-        unwrap.neut_intensity_thr = cfg['neutron_intensity_threshold']
-    if 'max_intensity_correlation_threshold' in cfg:
-        unwrap.neut_correlation_thr = cfg['max_intensity_correlation_threshold']
-    if 'trees_number' in cfg:
-        unwrap.trees_number = cfg['trees_number']
-    if 'max_branch_length' in cfg:
-        unwrap.max_branch_length = cfg['max_branch_length']
-    if 'pixel_spacing_ratio' in cfg:
-        unwrap.ratio_dxdy = cfg['pixel_spacing_ratio']
-    if 'initial_correlation_threshold' in cfg:
-        unwrap.init_corr_thr = cfg['initial_correlation_threshold']
-    if 'max_correlation_threshold' in cfg:
-        unwrap.max_corr_thr = cfg['max_correlation_threshold']
+    algorithm = cfg['algorithm']
     if 'correlation_threshold_increments' in cfg:
         unwrap.corr_incr_thr = cfg['correlation_threshold_increments']
-    if 'min_tile_area' in cfg:
-        unwrap.min_cc_area = cfg['min_tile_area']
-    if 'bootstrap_lines' in cfg:
-        unwrap.num_bs_lines = cfg['bootstrap_lines']
-    if 'min_overlap_area' in cfg:
-        unwrap.min_overlap_area = cfg['min_overlap_area']
-    if 'phase_variance_threshold' in cfg:
-        unwrap.phase_var_thr = cfg['phase_variance_threshold']
+
+    if algorithm == 'icu':
+        icu_cfg = cfg['icu']
+        if 'buffer_lines' in icu_cfg:
+            unwrap.buffer_lines = icu_cfg['buffer_lines']
+        if 'overlap_lines' in icu_cfg:
+            unwrap.overlap_lines = icu_cfg['overlap_lines']
+        if 'use_phase_gradient_neutron' in icu_cfg:
+            unwrap.use_phase_grad_neut = icu_cfg['use_phase_gradient_neutron']
+        if 'use_intensity_neutron' in icu_cfg:
+            unwrap.use_intensity_neut = icu_cfg['use_intensity_neutron']
+        if 'phase_gradient_window_size' in icu_cfg:
+            unwrap.phase_grad_win_size = icu_cfg['phase_gradient_window_size']
+        if 'neutron_phase_gradient_threshold' in icu_cfg:
+            unwrap.neut_phase_grad_thr = icu_cfg[
+                'neutron_phase_gradient_threshold']
+        if 'neutron_intensity_threshold' in icu_cfg:
+            unwrap.neut_intensity_thr = icu_cfg['neutron_intensity_threshold']
+        if 'max_intensity_correlation_threshold' in icu_cfg:
+            unwrap.neut_correlation_thr = icu_cfg[
+                'max_intensity_correlation_threshold']
+        if 'trees_number' in icu_cfg:
+            unwrap.trees_number = icu_cfg['trees_number']
+        if 'max_branch_length' in icu_cfg:
+            unwrap.max_branch_length = icu_cfg['max_branch_length']
+        if 'pixel_spacing_ratio' in icu_cfg:
+            unwrap.ratio_dxdy = icu_cfg['pixel_spacing_ratio']
+        if 'initial_correlation_threshold' in icu_cfg:
+            unwrap.init_corr_thr = icu_cfg['initial_correlation_threshold']
+        if 'max_correlation_threshold' in icu_cfg:
+            unwrap.max_corr_thr = icu_cfg['max_correlation_threshold']
+        if 'min_tile_area' in icu_cfg:
+            unwrap.min_cc_area = icu_cfg['min_tile_area']
+        if 'bootstrap_lines' in icu_cfg:
+            unwrap.num_bs_lines = icu_cfg['bootstrap_lines']
+        if 'min_overlap_area' in icu_cfg:
+            unwrap.min_overlap_area = icu_cfg['min_overlap_area']
+        if 'phase_variance_threshold' in icu_cfg:
+            unwrap.phase_var_thr = icu_cfg['phase_variance_threshold']
+    elif algorithm == 'phass':
+        phass_cfg = cfg['phass']
+        if 'good_correlation' in phass_cfg:
+            unwrap.good_correlation = phass_cfg['good_correlation']
+        if 'min_unwrap_area' in phass_cfg:
+            unwrap.min_pixels_region = phass_cfg['min_unwrap_area']
+    else:
+        err_str = "Not a valid unwrapping algorithm"
+        error_channel.log(err_str)
+        raise ValueError(err_str)
 
     return unwrap
 
@@ -187,13 +271,8 @@ def set_unwrap_attributes(unwrap, cfg: dict):
     # Prepare RUNW HDF5
     out_paths = h5_prep.run(unwrap_runconfig.cfg)
 
-    # out_paths['RIFG'] from h5_prep does not have actual interferogram data.
-    # Use RIFG path from CLI or YAML
-    if args.run_config_path is None:
-        rifg = args.crossmul
-    else:
-        rifg = unwrap_runconfig.cfg['processing']['phase_unwrap'][
-            'crossmul_path']
+    # Use RIFG from crossmul_path
+    rifg_h5 = unwrap_runconfig.cfg['processing']['phase_unwrap']['crossmul_path']
 
     # Run phase unwrapping
-    run(unwrap_runconfig.cfg, rifg, out_paths['RUNW'])
+    run(unwrap_runconfig.cfg, rifg_h5, out_paths['RUNW'])

python/packages/pybind_nisar/workflows/unwrap_runconfig.py

@@ -20,16 +20,17 @@ def yaml_check(self):
         '''
         error_channel = journal.error('CrossmulRunConfig.yaml_check')
 
-        if self.cfg['processing']['phase_unwrap'] is None:
-            err_str = "'phase_unwrap' necessary for standalone execution with YAML"
-            error_channel.log(err_str)
-            raise ValueError(err_str)
-
-        if 'crossmul_path' not in self.cfg['processing']['phase_unwrap']:
+        # Check if crossmul_path is provided (needed for stand-alone unwrapping)
+        if self.cfg['processing']['phase_unwrap']['crossmul_path'] is None:
             err_str = "'crossmul_path' file path under `phase_unwrap' required for standalone execution with YAML"
             error_channel.log(err_str)
             raise ValueError(err_str)
 
+        # Allocate, if not present, cfg related to user-unwrapper choice
+        algorithm = self.cfg['processing']['phase_unwrap']['algorithm']
+        if algorithm not in self.cfg['processing']['phase_unwrap']:
+            self.cfg['processing']['phase_unwrap'][algorithm] = {}
+
         # Check if crossmul path is a directory or a file
         crossmul_path = self.cfg['processing']['phase_unwrap']['crossmul_path']
         if not os.path.isfile(crossmul_path):

share/nisar/defaults/insar.yaml

@@ -271,3 +271,15 @@ runconfig:
                 flatten: True
                 oversample: 2
                 rows_per_block: 8192
+
+            phase_unwrap:
+                # Path to HDF5 file or directory containing the input interferogram
+                # and coherence (normalized magnitude of complex correlation) rasters
+                # Not required as we allow the use of intermediate outputs from
+                # the previous InSAR module, which is not user-specified.
+                # Required for running stand-alone phase unwrapping
+                crossmul_path:
+                # Increments to correlation threshold
+                correlation_threshold_increments: 0.1
+                # Unwrapping algorithm to use
+                algorithm: icu