Compute the phase interferometric wavenumber k_z

python/packages/CMakeLists.txt

@@ -39,6 +39,7 @@ set (list_of_exe
      nisar/workflows/gslc.py
      nisar/workflows/gslc_point_target_analysis.py
      nisar/workflows/insar.py
+     nisar/workflows/insar_kz.py
      nisar/workflows/nisar_l0b_dm1_to_science.py
      nisar/workflows/nisar_l0b_dm2_to_dbf.py
      nisar/workflows/noise_estimator.py

python/packages/nisar/workflows/insar.py

@@ -2,10 +2,12 @@
 import time
 
 import journal
+import isce3
 from nisar.workflows import (bandpass_insar, baseline, crossmul, dense_offsets,
                              filter_interferogram, geo2rdr, geocode_insar,
                              h5_prep, ionosphere, offsets_product,
-                             prepare_insar_hdf5, rdr2geo, resample_slc_v2,
+                             prepare_insar_hdf5, rdr2geo,
+                             resample_slc, resample_slc_v2,
                              rubbersheet, solid_earth_tides, split_spectrum,
                              troposphere, unwrap)
 from nisar.workflows.geocode_insar import InputProduct
@@ -36,7 +38,13 @@ def run(cfg: dict, out_paths: dict, run_steps: dict):
         prepare_insar_hdf5.run(cfg)
 
     if run_steps['coarse_resample']:
-        resample_slc_v2.run(cfg, 'coarse')
+        # Resample SLC V2 doesn't have GPU implemented, so run the old version if
+        # GPU was requested.
+        if isce3.core.gpu_check.use_gpu(cfg['worker']['gpu_enabled'],
+                                           cfg['worker']['gpu_id']):
+            resample_slc.run(cfg, 'coarse')
+        else:
+            resample_slc_v2.run(cfg, 'coarse')
 
     if (run_steps['dense_offsets']) and \
             (cfg['processing']['dense_offsets']['enabled']):
@@ -57,7 +65,13 @@ def run(cfg: dict, out_paths: dict, run_steps: dict):
         and cfg['processing']['fine_resample']['enabled']
         and 'RIFG' in out_paths
     ):
-        resample_slc_v2.run(cfg, 'fine')
+        # Resample SLC V2 doesn't have GPU implemented, so run the old version if
+        # GPU was requested.
+        if isce3.core.gpu_check.use_gpu(cfg['worker']['gpu_enabled'],
+                                           cfg['worker']['gpu_id']):
+            resample_slc.run(cfg, 'fine')
+        else:
+            resample_slc_v2.run(cfg, 'fine')
 
     # If fine_resampling is enabled, use fine-coregistered SLC
     # to run crossmul

python/packages/nisar/workflows/insar_kz.py

@@ -0,0 +1,173 @@
+#!/usr/bin/env python3
+import time
+
+import journal
+import tempfile
+from osgeo import gdal
+import numpy as np
+import plant
+
+from nisar.workflows import (bandpass_insar, crossmul, dense_offsets, geo2rdr,
+                             geocode_insar, h5_prep, filter_interferogram,
+                             offsets_product, prepare_insar_hdf5, rdr2geo,
+                             resample_slc, rubbersheet,
+                             split_spectrum, unwrap, ionosphere, baseline,
+                             troposphere, solid_earth_tides)
+
+from nisar.workflows.geocode_insar import InputProduct
+from nisar.workflows.insar_runconfig import InsarRunConfig
+from nisar.workflows.persistence import Persistence
+from nisar.workflows.yaml_argparse import YamlArgparse
+
+
+def _generate_dem_with_offset(dem_file, dem_offset, scratch_dir):
+    '''
+    Generate DEM with offset. It may require:
+
+    >> export GDAL_VRT_ENABLE_PYTHON=YES
+    '''
+
+    dem_with_offset = tempfile.NamedTemporaryFile(
+        dir=scratch_dir, suffix='.vrt').name
+
+    print('building DEM VRT with height offset [m]:', dem_offset)
+    gdal.BuildVRT(dem_with_offset, dem_file)
+
+    pixel_function_str = (
+        '        <PixelFunctionType>offset</PixelFunctionType>\n'
+        '        <PixelFunctionLanguage>Python</PixelFunctionLanguage>\n'
+        f'        <PixelFunctionArguments factor="{dem_offset}"/>\n'
+        '        <PixelFunctionCode><![CDATA[\n'
+        'import numpy as np\n'
+        'def offset(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,\n'
+        ' raster_ysize, buf_radius, gt, **kwargs):\n'
+        '    factor = float(kwargs["factor"])\n'
+        '    out_ar[:] = in_ar[0] + factor]]>\n'
+        '        </PixelFunctionCode>\n')
+
+    with open(dem_with_offset, "r") as f:
+        lines = f.readlines()
+
+    with open(dem_with_offset, "w") as f:
+        for line in lines:
+            if '<VRTRasterBand' in line:
+                line = line.replace('>', ' subClass="VRTDerivedRasterBand">')
+            if '<SimpleSource>' in line or '<ComplexSource>' in line:
+                f.write(pixel_function_str)
+            f.write(line)
+
+    return dem_with_offset
+
+
+def run(cfg: dict, out_paths: dict, run_steps: dict):
+    '''
+    Run INSAR workflow with parameters in cfg dictionary
+    '''
+    info_channel = journal.info("insar.run")
+    info_channel.log("starting INSAR")
+
+    t_all = time.time()
+
+    # if run_steps['bandpass_insar']:
+    #     bandpass_insar.run(cfg)
+    scratch_dir = cfg['product_path_group']['scratch_path']
+
+    # save original parameters
+    product_type_orig = cfg['primary_executable']['product_type']
+    print('product_type_orig:', product_type_orig)
+    # scratch_dir
+    # rdr2geo, geo2rdr, coarse_resample enabled
+
+    # k_z is computed as the InSAR phase difference using dem + 0.5 and
+    # dem - 0.5
+
+    dem_file = cfg['dynamic_ancillary_file_group']['dem_file']
+
+    for dem_offset in [0.5, 0, -0.5]:
+
+        if dem_offset == 0.5:
+            dem_offset_str = '_p0.5'
+        elif dem_offset == -0.5:
+            dem_offset_str = '_m0.5'
+        elif dem_offset == 0:
+            dem_offset_str = '_zero'
+        else:
+            raise NotImplementedError("Invalid option")
+
+        # udpate dem with dem_offset, create VRT and update dem_file ?
+
+        # write new parameters
+        cfg['primary_executable']['product_type'] = 'RIFG'
+
+        dem_with_offset = _generate_dem_with_offset(dem_file, dem_offset,
+                                                    scratch_dir)
+
+        cfg['dynamic_ancillary_file_group']['dem_file'] = dem_with_offset
+
+        # create RFIG writer
+        if run_steps['prepare_insar_hdf5']:
+            prepare_insar_hdf5.run(cfg)
+
+        if run_steps['rdr2geo']:
+            rdr2geo.run(cfg)
+
+        if run_steps['geo2rdr']:
+            geo2rdr.run(cfg)
+
+        if run_steps['coarse_resample']:
+            resample_files_dict = resample_slc.run(cfg, 'coarse', flatten=True,
+                                                   flag_constant_value=True,
+                                                   suffix=dem_offset_str)
+
+            if dem_offset == 0.5:
+                resample_files_dict_p_05 = resample_files_dict
+            elif dem_offset == -0.5:
+                resample_files_dict_m_05 = resample_files_dict
+
+    for (_, pol_dict_p_05), (_, pol_dict_m_05) in \
+            zip(resample_files_dict_p_05.items(),
+                resample_files_dict_m_05.items()):
+
+        for (_, f_p_05), (_, f_m_05) \
+                in zip(pol_dict_p_05.items(), pol_dict_m_05.items()):
+
+            print(f'DEM +0.5 file: {f_p_05}')
+            print(f'DEM -0.5 file: {f_m_05}')
+            data_p_05_band = gdal.Open(f_p_05, gdal.GA_ReadOnly)
+            data_p_05 = data_p_05_band.GetRasterBand(1).ReadAsArray()
+            data_m_05_band = gdal.Open(f_m_05, gdal.GA_ReadOnly)
+            data_m_05 = data_m_05_band.GetRasterBand(1).ReadAsArray()
+            data_diff = np.angle(data_p_05 * np.conj(data_m_05))
+            output_filename = f_p_05.replace('_p0.5', '_kz')
+            plant.save_image(data_diff, output_filename, force=True)
+            print(f'file saved: {output_filename}')
+
+        print('*** resample_files_dict:', resample_files_dict)
+
+    # if run_steps['baseline']:
+    #     baseline.run(cfg, out_paths)
+
+    t_all_elapsed = time.time() - t_all
+    info_channel.log(f"successfully ran INSAR in {t_all_elapsed:.3f} seconds")
+
+
+if __name__ == "__main__":
+    # parse CLI input
+    yaml_parser = YamlArgparse()
+    args = yaml_parser.parse()
+
+    # convert CLI input to run configuration
+    insar_runcfg = InsarRunConfig(args)
+
+    # To allow persistence, a logfile is required. Raise exception
+    # if logfile is None and persistence is requested
+    logfile_path = insar_runcfg.cfg['logging']['path']
+    if (logfile_path is None) and insar_runcfg.args.restart:
+        raise ValueError('InSAR workflow persistence requires to specify a logfile')
+    persist = Persistence(logfile_path, insar_runcfg.args.restart)
+
+    # run InSAR workflow
+    if persist.run:
+        _, out_paths = h5_prep.get_products_and_paths(insar_runcfg.cfg)
+
+        run(insar_runcfg.cfg, out_paths, persist.run_steps)

python/packages/nisar/workflows/resample_slc.py

@@ -7,17 +7,54 @@
 import os
 import pathlib
 import time
+import tempfile
 
 import journal
 import isce3
 from osgeo import gdal
 from nisar.products.readers import SLC
-from nisar.workflows.helpers import copy_raster
+from nisar.workflows.helpers import copy_raster, complex_raster_path_from_h5
 from nisar.workflows.resample_slc_runconfig import ResampleSlcRunConfig
 from nisar.workflows.yaml_argparse import YamlArgparse
 
 
-def run(cfg, resample_type):
+def build_empty_vrt(filename, length, width, fill_value, dtype='CFloat32'):
+    """Build an empty VRT file, i.e, not pointing to any rasters,
+       with given input dimensions (length and width), data type, and
+       fill value.
+       Parameters
+       ----------
+       filename: str
+              VRT file name
+       length: int
+              VRT data length
+       width: int
+              VRT data width
+       fill_value: scalar
+              VRT data fill value
+       dtype: str
+              VRT data type
+    """
+
+    vrt_contents = f'<VRTDataset rasterXSize="{width}"'
+    vrt_contents += f' rasterYSize="{length}"> \n'
+
+    vrt_contents += (
+        f'  <VRTRasterBand dataType="{dtype}" band="1"> \n'
+        f'    <NoDataValue>{fill_value}</NoDataValue> \n'
+        f'    <HideNoDataValue>{fill_value}</HideNoDataValue> \n'
+        f'  </VRTRasterBand> \n'
+        f'</VRTDataset> \n')
+
+    with open(filename, 'w') as out:
+        out.write(vrt_contents)
+
+    if os.path.isfile(filename):
+        print('file saved:', filename)
+
+
+def run(cfg, resample_type, flatten=False, flag_constant_value=False,
+        suffix=''):
     '''
     run resample_slc
     '''
@@ -30,6 +67,10 @@ def run(cfg, resample_type):
 
     # Get SLC parameters
     slc = SLC(hdf5file=input_hdf5)
+    # If flattening is enabled, load the reference SLC
+    if flatten:
+        reference_hdf5 = cfg['input_file_group']['reference_rslc_file']
+        ref_slc = SLC(hdf5file=reference_hdf5)
 
     info_channel = journal.info('resample_slc.run')
     info_channel.log('starting resampling SLC')
@@ -45,7 +86,12 @@ def run(cfg, resample_type):
 
     t_all = time.time()
 
+    resample_files_dict = {}
+
     for freq in freq_pols.keys():
+
+        resample_files_dict[freq] = {}
+
         # Get frequency specific parameters
         radar_grid = slc.getRadarGrid(frequency=freq)
         native_doppler = slc.getDopplerCentroid(frequency=freq)
@@ -79,7 +125,14 @@ def run(cfg, resample_type):
         else:
             Resamp = isce3.image.ResampSlc
 
-        resamp_obj = Resamp(radar_grid, native_doppler)
+        # If flattening is enabled, add the reference SLC radar grid to
+        # the call to the constructor of the ResampSlc module
+        if flatten:
+            ref_radar_grid = ref_slc.getRadarGrid(frequency=freq)
+            resamp_obj = Resamp(radar_grid, native_doppler,
+                                ref_rdr_grid=ref_radar_grid)
+        else:
+            resamp_obj = Resamp(radar_grid, native_doppler)
 
         # If lines per tile is > 0, assign it to resamp_obj
         if resamp_args['lines_per_tile']:
@@ -92,19 +145,46 @@ def run(cfg, resample_type):
             # Create directory for each polarization
             out_dir = resample_slc_scratch_path / pol
             out_dir.mkdir(parents=True, exist_ok=True)
-            out_path = out_dir / 'coregistered_secondary.slc'
+
+            out_path = out_dir / f'coregistered_secondary{suffix}.slc'
+
+            # If necessary, perform complex32 to complex64 conversion on input
+            input_as_c32_path = str(out_dir/'secondary.slc')
+            input_raster_path, _ = complex_raster_path_from_h5(
+                slc, freq, pol, input_hdf5, resamp_args['lines_per_tile'],
+                input_as_c32_path)
+
+            input_raster = isce3.io.Raster(input_raster_path)
 
             # Dump secondary RSLC on disk
-            raster_path = str(out_dir/'secondary.slc')
-            copy_raster(input_hdf5, freq, pol, 1000,
-                        raster_path, file_type='ENVI')
-            input_raster = isce3.io.Raster(raster_path)
+            # raster_path = str(out_dir/'secondary.slc')
+            # copy_raster(input_hdf5, freq, pol, 1000,
+            #            raster_path, file_type='ENVI')
+            # input_raster = isce3.io.Raster(raster_path)
+
+            # Dump secondary RSLC on disk
+            # raster_path = str(out_dir/'secondary.slc')
+            # copy_raster(input_hdf5, freq, pol, 1000,
+            #            raster_path, file_type='ENVI')
+            # input_raster = isce3.io.Raster(raster_path)
+
+            if flag_constant_value:
+                print('*** constant value!')
+                fill_value = 1
+                input_raster_vrt = tempfile.NamedTemporaryFile(
+                    dir=scratch_path, suffix='.vrt').name
+                build_empty_vrt(input_raster_vrt, input_raster.length,
+                                input_raster.width, fill_value)
+                input_raster = isce3.io.Raster(input_raster_vrt)
 
             # Create output raster
             resamp_slc = isce3.io.Raster(str(out_path), rg_off.width,
                                          rg_off.length, rg_off.num_bands,
                                          gdal.GDT_CFloat32, 'ENVI')
-            resamp_obj.resamp(input_raster, resamp_slc, rg_off, az_off)
+            resamp_obj.resamp(input_raster, resamp_slc, rg_off, az_off,
+                              flatten=flatten)
+
+            resample_files_dict[freq][pol] = str(out_path)
 
     t_all_elapsed = time.time() - t_all
     info_channel.log(f"successfully ran resample in {t_all_elapsed:.3f} seconds")