diff --git a/src/compass/s1_geocode_metadata.py b/src/compass/s1_geocode_metadata.py index 73775374..ccec8fba 100755 --- a/src/compass/s1_geocode_metadata.py +++ b/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) diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py index 61926ecb..769d2bf0 100755 --- a/src/compass/s1_geocode_slc.py +++ b/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)] diff --git a/src/compass/utils/browse_image.py b/src/compass/utils/browse_image.py index ba0ef807..e05324aa 100644 --- a/src/compass/utils/browse_image.py +++ b/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 diff --git a/src/compass/utils/fill_value.py b/src/compass/utils/fill_value.py new file mode 100644 index 00000000..d31c20e9 --- /dev/null +++ b/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) diff --git a/src/compass/utils/h5_helpers.py b/src/compass/utils/h5_helpers.py index cb0abe17..96eb719e 100644 --- a/src/compass/utils/h5_helpers.py +++ b/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,22 +80,26 @@ 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 @@ -102,6 +107,8 @@ def init_geocoded_dataset(grid_group, dataset_name, geo_grid, dtype, 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