ISO metadata for InSAR product (#935)

* Add ISO metadata to RIFG layers

* Add ISO metadata to RUNW layers

* Generate stats in RUNW instead of copying

* Create external compute_stats function

* Use external compute_stats inside workflow

* Reusing the compute_stats function in SLC writer

Author: Virginia Brancato

python/packages/nisar/products/writers/SLC.py

@@ -11,6 +11,8 @@
 from nisar.types import complex32
 from nisar.products.readers.Raw import Raw
 from nisar.workflows.h5_prep import add_geolocation_grid_cubes_to_hdf5
+from nisar.workflows.compute_stats import compute_stats_complex_data
+
 
 log = logging.getLogger("SLCWriter")
 
@@ -378,17 +380,7 @@ def compute_stats(self, frequency, pol):
         '''
         h5_ds = self.swath(frequency)[pol]
         raster = isce3.io.Raster(f"IH5:::ID={h5_ds.id.id}".encode("utf-8"))
-        stats_obj = isce3.math.compute_raster_stats_real_imag(raster)[0]
-        h5_ds.attrs.create('min_real_value', data=stats_obj.min_real)
-        h5_ds.attrs.create('mean_real_value', data=stats_obj.mean_real)
-        h5_ds.attrs.create('max_real_value', data=stats_obj.max_real)
-        h5_ds.attrs.create('sample_stddev_real', 
-                           data=stats_obj.sample_stddev_real)
-        h5_ds.attrs.create('min_imag_value', data=stats_obj.min_imag)
-        h5_ds.attrs.create('mean_imag_value', data=stats_obj.mean_imag)
-        h5_ds.attrs.create('max_imag_value', data=stats_obj.max_imag)
-        h5_ds.attrs.create('sample_stddev_imag',
-                           data=stats_obj.sample_stddev_imag)
+        compute_stats_complex_data(raster, h5_ds)
 
 
     def add_calibration_section(self, frequency, pol, 

python/packages/nisar/workflows/compute_stats.py

@@ -0,0 +1,44 @@
+import isce3
+from osgeo import gdal
+
+
+def compute_stats_complex_data(raster, h5_ds):
+    """
+    Compute raster statistics for complex datasets
+
+    raster: isce3.io.Raster
+        ISCE3 Raster object
+    h5_ds: h5py.File
+        h5py file
+    """
+    stats_obj = isce3.math.compute_raster_stats_real_imag(raster)[0]
+    h5_ds.attrs.create('min_real_value', data=stats_obj.min_real)
+    h5_ds.attrs.create('mean_real_value', data=stats_obj.mean_real)
+    h5_ds.attrs.create('max_real_value', data=stats_obj.max_real)
+    h5_ds.attrs.create('sample_stddev_real',
+                       data=stats_obj.sample_stddev_real)
+    h5_ds.attrs.create('min_imag_value', data=stats_obj.min_imag)
+    h5_ds.attrs.create('mean_imag_value', data=stats_obj.mean_imag)
+    h5_ds.attrs.create('max_imag_value', data=stats_obj.max_imag)
+    h5_ds.attrs.create('sample_stddev_imag',
+                       data=stats_obj.sample_stddev_imag)
+
+
+def compute_stats_real_data(raster, h5_ds):
+    """
+    Compute raster statistics for real datasets
+
+    raster: isce3.io.Raster
+       ISCE3 Raster object
+    h5_ds: h5py.File
+        h5py file
+    """
+    if raster.datatype() == gdal.GDT_Float64:
+        stats_obj = isce3.math.compute_raster_stats_float64(raster)[0]
+    else:
+        stats_obj = isce3.math.compute_raster_stats_float32(raster)[0]
+    h5_ds.attrs.create('min_value', data=stats_obj.min)
+    h5_ds.attrs.create('mean_value', data=stats_obj.mean)
+    h5_ds.attrs.create('max_value', data=stats_obj.max)
+    h5_ds.attrs.create('sample_stddev',
+                       data=stats_obj.sample_stddev)

python/packages/nisar/workflows/crossmul.py

@@ -9,10 +9,13 @@
 import h5py
 import journal
 import isce3
+
+from osgeo import gdal
 from nisar.products.readers import SLC
 from nisar.workflows import h5_prep
 from nisar.workflows.crossmul_runconfig import CrossmulRunConfig
 from nisar.workflows.yaml_argparse import YamlArgparse
+from nisar.workflows.compute_stats import compute_stats_real_data
 
 
 def run(cfg: dict, output_hdf5: str = None, resample_type='coarse'):
@@ -137,19 +140,45 @@ def run(cfg: dict, output_hdf5: str = None, resample_type='coarse'):
                         crossmul.crossmul(ref_slc_raster, sec_slc_raster,
                                           igram_raster, coherence_raster)
 
+                    # Allocate raster statistics for coherence
+                    compute_stats_real_data(coherence_raster, coherence_dataset)
+
                     del coherence_raster
                 else:
                     # no coherence without multilook
                     crossmul.crossmul(ref_slc_raster, sec_slc_raster,
                                       igram_raster)
-
                 del igram_raster
 
+                # Allocate stats for rubbersheet offsets
+                stats_offsets(dst_h5, freq, pol)
+
     t_all_elapsed = time.time() - t_all
     info_channel.log(
         f"successfully ran crossmul in {t_all_elapsed:.3f} seconds")
 
 
+def stats_offsets(h5_ds, freq, pol):
+    """
+    Allocate statistics for dense offsets
+    h5_ds: h5py.File
+       h5py File
+    freq: str
+       Frequency to process (A or B)
+    pol: str
+       Polarization to process (HH, HV, VH, VV)
+    """
+    path = f'science/LSAR/RIFG/swaths/frequency{freq}/pixelOffsets/{pol}/'
+    offset_layer = ['slantRangeOffset', 'alongTrackOffset']
+
+    for layer in offset_layer:
+        offset_path = f'{path}/{layer}'
+        offset_dataset = h5_ds[offset_path]
+        offset_raster = isce3.io.Raster(
+            f"IH5:::ID={offset_dataset.id.id}".encode("utf-8"))
+        compute_stats_real_data(offset_raster, offset_dataset)
+
+
 if __name__ == "__main__":
     '''
     run crossmul from command line

python/packages/nisar/workflows/geocode_insar.py

@@ -19,6 +19,7 @@
 from nisar.workflows.geocode_insar_runconfig import \
     GeocodeInsarRunConfig
 from nisar.workflows.yaml_argparse import YamlArgparse
+from nisar.workflows.compute_stats import compute_stats_real_data
 
 
 def run(cfg, runw_hdf5, output_hdf5):
@@ -357,7 +358,7 @@ def cpu_run(cfg, runw_hdf5, output_hdf5):
                     if (dataset_name != "layoverShadowMask"):
                         # Layover/shadow masks dont't have min/max/mean/stddev
                         # stats attributes
-                        _compute_stats(geocoded_raster, geocoded_dataset)
+                        compute_stats_real_data(geocoded_raster, geocoded_dataset)
 
             # spec for NISAR GUNW does not require freq B so skip radar cube
             if freq.upper() == 'B':
@@ -565,8 +566,7 @@ def gpu_run(cfg, runw_hdf5, output_hdf5):
                         if (dataset_name != "layoverShadowMask"):
                             # Layover/shadow masks dont't have min/max/mean/stddev
                             # stats attributes
-                            _compute_stats(geocoded_raster, geocoded_dataset)
-
+                            compute_stats_real_data(geocoded_raster, geocoded_dataset)
 
                 if gunw_datasets['layoverShadowMask']:
                     skip_layover_shadow = True
@@ -581,18 +581,6 @@ def gpu_run(cfg, runw_hdf5, output_hdf5):
     info_channel.log(f"Successfully ran geocode in {t_all_elapsed:.3f} seconds")
 
 
-def _compute_stats(raster, h5_ds):
-    if raster.datatype() == gdal.GDT_Float64:
-        stats_obj = isce3.math.compute_raster_stats_float64(raster)[0]
-    else:
-        stats_obj = isce3.math.compute_raster_stats_float32(raster)[0]
-    h5_ds.attrs.create('min_value', data=stats_obj.min)
-    h5_ds.attrs.create('mean_value', data=stats_obj.mean)
-    h5_ds.attrs.create('max_value', data=stats_obj.max)
-    h5_ds.attrs.create('sample_stddev',
-                       data=stats_obj.sample_stddev)
-
-
 if __name__ == "__main__":
     """
     run geocode from command line

python/packages/nisar/workflows/gslc.py

@@ -12,6 +12,7 @@
 from nisar.workflows.h5_prep import add_radar_grid_cubes_to_hdf5
 from nisar.workflows.yaml_argparse import YamlArgparse
 from nisar.workflows.gslc_runconfig import GSLCRunConfig
+from nisar.workflows.compute_stats import compute_stats_complex_data
 
 
 def run(cfg):
@@ -89,7 +90,7 @@ def run(cfg):
 
                 # output_raster_ref = f'HDF5:{output_hdf5}:/{dataset_path}'
                 gslc_raster = isce3.io.Raster(f"IH5:::ID={gslc_dataset.id.id}".encode("utf-8"))
-                _compute_stats(gslc_raster, gslc_dataset)
+                compute_stats_complex_data(gslc_raster, gslc_dataset)
 
                 t_pol_elapsed = time.time() - t_pol
                 info_channel.log(f'polarization {polarization} ran in {t_pol_elapsed:.3f} seconds')
@@ -123,20 +124,6 @@ def run(cfg):
     info_channel.log(f"successfully ran geocode SLC in {t_all_elapsed:.3f} seconds")
 
 
-def _compute_stats(raster, h5_ds):
-    stats_obj = isce3.math.compute_raster_stats_real_imag(raster)[0]
-    h5_ds.attrs.create('min_real_value', data=stats_obj.min_real)
-    h5_ds.attrs.create('mean_real_value', data=stats_obj.mean_real)
-    h5_ds.attrs.create('max_real_value', data=stats_obj.max_real)
-    h5_ds.attrs.create('sample_stddev_real', 
-                       data=stats_obj.sample_stddev_real)
-    h5_ds.attrs.create('min_imag_value', data=stats_obj.min_imag)
-    h5_ds.attrs.create('mean_imag_value', data=stats_obj.mean_imag)
-    h5_ds.attrs.create('max_imag_value', data=stats_obj.max_imag)
-    h5_ds.attrs.create('sample_stddev_imag',
-                       data=stats_obj.sample_stddev_imag)
-
-
 if __name__ == "__main__":
     yaml_parser = YamlArgparse()
     args = yaml_parser.parse()

python/packages/nisar/workflows/unwrap.py

@@ -17,6 +17,8 @@
 from nisar.workflows import h5_prep
 from nisar.workflows.unwrap_runconfig import UnwrapRunConfig
 from nisar.workflows.yaml_argparse import YamlArgparse
+from nisar.workflows.compute_stats import compute_stats_real_data
+
 
 def run(cfg: dict, input_hdf5: str, output_hdf5: str):
     '''
@@ -143,7 +145,10 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
                     unwrap_obj.unwrap(uigram_raster, conn_comp_raster, igram_raster,
                                       corr_raster)
 
-                # Copy coherence magnitude and culled offsets to RIFG
+                # Compute stats for unwrapped phase
+                compute_stats_real_data(uigram_raster, uigram_dataset)
+
+                # Copy coherence magnitude and culled offsets from RIFG
                 dataset_names = ['coherenceMagnitude', 'alongTrackOffset',
                                  'slantRangeOffset']
                 group_names = ['interferogram', 'pixelOffsets', 'pixelOffsets']
@@ -152,6 +157,13 @@ def run(cfg: dict, input_hdf5: str, output_hdf5: str):
                     src_path = f'{src_freq_group_path}/{group_name}/{pol}/{dataset_name}'
                     dst_h5[dst_path][:, :] = src_h5[src_path][()]
 
+                    dst_dataset = dst_h5[dst_path]
+                    dst_raster = isce3.io.Raster(
+                        f"IH5:::ID={dst_dataset.id.id}".encode("utf-8"),
+                        update=True)
+
+                    compute_stats_real_data(dst_raster, dst_dataset)
+
                 del uigram_raster
                 del conn_comp_raster