Correctly populate product fill values in geocoded datasets

src/compass/s1_geocode_metadata.py

@@ -64,7 +64,7 @@ def _fix_layover_shadow_mask(static_layers_dict, h5_root, geo_grid,
         del h5_root[layover_shadow_path]
         desc = 'Layover shadow mask. 0=no layover, no shadow; 1=shadow; 2=layover; 3=shadow and layover.'
         _ = init_geocoded_dataset(h5_root[DATA_PATH], dst_ds_name, geo_grid,
-                                  dtype=None,
+                                  dtype=np.int8,
                                   description=np.string_(desc),
                                   data=temp_arr, output_cfg=output_params)
 
@@ -473,10 +473,8 @@ def geocode_noise_luts(geo_burst_h5, burst, cfg,
     for _, bursts in bursts_grouping_generator(cfg.bursts):
         burst = bursts[0]
 
-        # Generate required static layers
-        if cfg.rdr2geo_params.enabled:
-            s1_rdr2geo.run(cfg, save_in_scratch=True)
+        # Generate required un-geocoded static layers
+        s1_rdr2geo.run(cfg, save_in_scratch=True)
 
-            # Geocode static layers if needed
-            if cfg.rdr2geo_params.geocode_metadata_layers:
-                run(cfg, burst, fetch_from_scratch=True)
+        # Geocode static layers generated in step above
+        run(cfg, burst, fetch_from_scratch=True)

src/compass/s1_geocode_slc.py

@@ -129,7 +129,7 @@ def run(cfg: GeoRunConfig):
             ctype = h5py.h5t.py_create(np.complex64)
             ctype.commit(geo_burst_h5['/'].id, np.string_('complex64'))
 
-            grid_group = geo_burst_h5.require_group(DATA_PATH)
+            data_group = geo_burst_h5.require_group(DATA_PATH)
             check_eap = is_eap_correction_necessary(burst.ipf_version)
 
             # Initialize source/radar and destination/geo dataset into lists
@@ -167,7 +167,7 @@ def run(cfg: GeoRunConfig):
                 rdr_data_blks.append(rdr_dataset.ReadAsArray())
 
                 # Prepare output dataset of current polarization in HDF5
-                geo_ds = init_geocoded_dataset(grid_group, pol, geo_grid,
+                geo_ds = init_geocoded_dataset(data_group, pol, geo_grid,
                                                'complex64',
                                                f'{pol} geocoded CSLC image',
                                                output_cfg=cfg.output_params)
@@ -193,7 +193,7 @@ def run(cfg: GeoRunConfig):
             ((carrier_phase_data_blk, carrier_phase_ds),
              (flatten_phase_data_blk, flatten_phase_ds)) = \
             [(np.full(out_shape, np.nan).astype(np.float64),
-                  init_geocoded_dataset(grid_group, ds_name, geo_grid,
+                  init_geocoded_dataset(data_group, ds_name, geo_grid,
                                         np.float64, ds_desc,
                                         output_cfg=cfg.output_params))
                  for ds_name, ds_desc in zip(phase_names, phase_descrs)]

src/compass/utils/browse_image.py

@@ -206,14 +206,14 @@ def make_browse_image(filename, path_h5, bursts, complex_to_real='amplitude', pe
     derived_netcdf_to_grid = f'{derived_ds_str}:NETCDF:{path_h5}:/{DATA_PATH}'
 
     with h5py.File(path_h5, 'r', swmr=True) as h5_obj:
-        grid_group = h5_obj[DATA_PATH]
+        data_group = h5_obj[DATA_PATH]
 
         for b in bursts:
             # get polarization to extract geocoded raster
             pol = b.polarization
 
             # compute browse shape
-            full_shape = grid_group[pol].shape
+            full_shape = data_group[pol].shape
             browse_h, browse_w = _scale_to_max_pixel_dimension(full_shape)
 
             # create in memory GDAL raster for GSLC as real value array

src/compass/utils/fill_value.py

@@ -0,0 +1,93 @@
+'''
+Class and function for helping set and determine dataset fill values
+'''
+from dataclasses import dataclass
+
+import journal
+import numpy as np
+
+
+@dataclass(frozen=True)
+class FillValues:
+    '''
+    Dataclass of fill values for float, complex, and int types with default
+    values.
+    '''
+    # Catch all fill value for all float types (float32, float64, ...)
+    float_fill: np.single = np.nan
+
+    # Catch all fill value for all complex types (complex64, complex128, ...)
+    complex_fill: np.csingle = np.nan * (0 + 1j)
+
+    # Catch all fill value for all int types (int8, byte8, ...)
+    # Currently hard coded for int8/layover_shadow
+    int_fill: np.intc = 127
+
+    @classmethod
+    def from_user_defined_value(cls, user_defined_value):
+        '''
+        Create and return a FillValues class object populated with all default
+        values populated to a single user defined value.
+
+        Parameters
+        ----------
+        user_defined_value: float
+            User defined value to be assigned to default value of all types
+
+        Returns
+        -------
+        FillValues
+            FillValues object with all default values set to user defined value
+        '''
+        return cls(np.single(user_defined_value),
+                   np.csingle(user_defined_value),
+                   np.intc(user_defined_value))
+
+
+def determine_fill_value(dtype, usr_fill_val=None):
+    '''
+    Helper function to determine COMPASS specific fill values based on h5py
+    Dataset type (dtype)
+
+    Parameters
+    ----------
+    dtype: type
+        Given numeric type whose corresponding fill value of same type is to be
+        determined
+    usr_fill_val: float
+        User specified non-default dataset fill value
+
+    Returns
+    -------
+        Fill value of type dtype. An exception is raised if no appropriate
+        value is found.
+    '''
+    if usr_fill_val is None:
+        fill_values = FillValues()
+    else:
+        # Assign user provided non-default value to all fill values in
+        # FillValues object with correct typing. Logic below will return on
+        # accordingly.
+        fill_values = FillValues.from_user_defined_value(usr_fill_val)
+
+    # Check if float type and return float fill
+    float_types = [np.double, np.single, np.float32, np.float64, 'float32',
+                   'float64']
+    if dtype in float_types:
+        return fill_values.float_fill
+
+    # Check if complex type and return complex fill
+    complex_types = [np.complex128, 'complex64', np.complex64, 'complex32']
+    if dtype in complex_types:
+        return fill_values.complex_fill
+
+    # Check if int type and return int fill
+    int_types = [np.byte, np.int8]
+    if dtype in int_types:
+        return fill_values.int_fill
+
+    # No appropriate fill value found above. Raise exception.
+    err_str = f'Unexpected COMPASS geocoded dataset type: {dtype}'
+    error_channel = journal.error('fill_value.determine_fill_value')
+    error_channel.log(err_str)
+    raise ValueError(err_str)

src/compass/utils/h5_helpers.py

@@ -14,6 +14,7 @@
 import shapely
 
 import compass
+from compass.utils.fill_value import determine_fill_value
 
 
 TIME_STR_FMT = '%Y-%m-%d %H:%M:%S.%f'
@@ -79,29 +80,35 @@ def add_dataset_and_attrs(group, meta_item):
         val_ds.attrs[key] = _as_np_string_if_needed(val)
 
 
-def init_geocoded_dataset(grid_group, dataset_name, geo_grid, dtype,
-                          description, data=None, output_cfg=None):
+def init_geocoded_dataset(data_group, dataset_name, geo_grid, dtype,
+                          description, data=None, output_cfg=None,
+                          fill_val=None):
     '''
     Create and allocate dataset for isce.geocode.geocode_slc to write to that
-    is CF-compliant
+    is CF-compliant. If data parameter not provided, then an appropriate fill
+    value is found and used to fill dataset.
 
     Parameters
     ----------
-    grid_group: h5py.Group
+    data_group: h5py.Group
         h5py group where geocoded dataset will be created in
     dataset_name: str
         Name of dataset to be created
     geo_grid: isce3.product.GeoGridParameters
         Geogrid of output
-    dtype: str
-        Data type of dataset to be geocoded
+    dtype: Union(str, type)
+        Data type of dataset to be geocoded to be passed to require_dataset.
+        require_dataset can take string values e.g. "float32" or types e.g.
+        numpy.float32
     description: str
         Description of dataset to be geocoded
     data: np.ndarray
         Array to set dataset raster with
     output_cfg: dict
         Optional dict containing output options in runconfig to apply to
         created datasets
+    fill_val: float
+        Optional value to fill an empty dataset
 
     Returns
     -------
@@ -120,15 +127,28 @@ def init_geocoded_dataset(grid_group, dataset_name, geo_grid, dtype,
         output_kwargs['compression_opts'] = output_cfg.compression_level
         output_kwargs['shuffle'] = output_cfg.shuffle
 
+    # Shape of dataset is defined by the geo grid
     shape = (geo_grid.length, geo_grid.width)
+
+    # Determine fill value of dataset to either correctly init empty dataset
+    # and/or populate dataset attribute
+    _fill_val = determine_fill_value(dtype, fill_val)
+
+    # If data is None, create dataset to specified parameters and fill with
+    # specified fill value. If data is not None, create a dataset with
+    # provided data.
     if data is None:
-        cslc_ds = grid_group.require_dataset(dataset_name, dtype=dtype,
-                                             shape=shape, **output_kwargs)
+        # Create a dataset with shape and a fill value from above
+        cslc_ds = data_group.require_dataset(dataset_name, dtype=dtype,
+                                             shape=shape, fillvalue=_fill_val,
+                                             **output_kwargs)
     else:
-        cslc_ds = grid_group.create_dataset(dataset_name, data=data,
+        # Create a dataset with provided data
+        cslc_ds = data_group.create_dataset(dataset_name, data=data,
                                             **output_kwargs)
 
     cslc_ds.attrs['description'] = description
+    cslc_ds.attrs['fill_value'] = _fill_val
 
     # Compute x scale
     dx = geo_grid.spacing_x
@@ -144,9 +164,9 @@ def init_geocoded_dataset(grid_group, dataset_name, geo_grid, dtype,
 
     # following copied and pasted (and slightly modified) from:
     # https://github-fn.jpl.nasa.gov/isce-3/isce/wiki/CF-Conventions-and-Map-Projections
-    x_ds = grid_group.require_dataset('x_coordinates', dtype='float64',
+    x_ds = data_group.require_dataset('x_coordinates', dtype='float64',
                                       data=x_vect, shape=x_vect.shape)
-    y_ds = grid_group.require_dataset('y_coordinates', dtype='float64',
+    y_ds = data_group.require_dataset('y_coordinates', dtype='float64',
                                       data=y_vect, shape=y_vect.shape)
 
     # Mapping of dimension scales to datasets is not done automatically in HDF5
@@ -160,6 +180,7 @@ def init_geocoded_dataset(grid_group, dataset_name, geo_grid, dtype,
     # Associate grid mapping with data - projection created later
     cslc_ds.attrs['grid_mapping'] = np.string_("projection")
 
+    # Build list of metadata to be inserted to accompany dataset
     grid_meta_items = [
         Meta('x_spacing', geo_grid.spacing_x,
              'Spacing of the geographical grid along X-direction',
@@ -169,14 +190,14 @@ def init_geocoded_dataset(grid_group, dataset_name, geo_grid, dtype,
              {'units': 'meters'})
     ]
     for meta_item in grid_meta_items:
-        add_dataset_and_attrs(grid_group, meta_item)
+        add_dataset_and_attrs(data_group, meta_item)
 
     # Set up osr for wkt
     srs = osr.SpatialReference()
     srs.ImportFromEPSG(geo_grid.epsg)
 
     #Create a new single int dataset for projections
-    projection_ds = grid_group.require_dataset('projection', (), dtype='i')
+    projection_ds = data_group.require_dataset('projection', (), dtype='i')
     projection_ds[()] = geo_grid.epsg
 
     # WGS84 ellipsoid