From b5e9341fdb3a1cd9a624a24410b1b05a215c48da Mon Sep 17 00:00:00 2001 From: Liang Yu Date: Mon, 5 Dec 2022 15:34:08 -0800 Subject: [PATCH] HDF5 output (#62) * untested geocode SLC to hdf5; no metadata * write metadata to hdf5 save yaml contents to text to save to metadata fix bugs from previous commit * add raster axis dataset * save topo layers to hdf5 * fix syntax error * remove commented code hdf5 to h5 file extension add docstring * more descriptive function name * add attributes to datasets * snake to camel case better description * description fixes * output directory and file paths generated in runconfig pep8 fixes more attr fixes * avoid gdal + IH5 + io.Raster crash * fix indentation error * cfg paths to rdr workflow more camel to snake * adopt burst ID class remove redundant attributes from backscatter remove np.string_ reorg metadata to reflect organization in specs * uncommenting commented out block of code * validation with HDF5 more metadata grouping to remove solo floaters * adjustments to runconfig formation for new schema bbox isnt used any more burst db now required burst_id is a list, not a string * spacing * removed JSON metadata added h5_helpers for hdf5 methods for dataset init and metadata writing * common time format string add EAP and more input file to metadata * Update src/compass/utils/h5_helpers.py correctly attach x+y ds scales Co-authored-by: Scott Staniewicz * equal spacing in y Co-authored-by: Scott Staniewicz * fix paths correctly place CF attr * fix indentation * add type at h5 root extension back to h5 Co-authored-by: Scott Staniewicz Co-authored-by: Scott Staniewicz --- src/compass/s1_geo2rdr.py | 13 +- src/compass/s1_geocode_metadata.py | 99 +++-- src/compass/s1_geocode_slc.py | 108 +++--- src/compass/s1_geocode_stack.py | 30 +- src/compass/s1_rdr2geo.py | 13 +- src/compass/s1_resample.py | 19 +- src/compass/utils/geo_grid.py | 2 +- src/compass/utils/geo_metadata.py | 278 -------------- src/compass/utils/geo_runconfig.py | 34 +- src/compass/utils/h5_helpers.py | 448 ++++++++++++++++++++++ src/compass/utils/range_split_spectrum.py | 2 +- src/compass/utils/raster_polygon.py | 15 +- src/compass/utils/runconfig.py | 84 +++- src/compass/utils/validate_cslc.py | 133 +++++-- 14 files changed, 766 insertions(+), 512 deletions(-) delete mode 100644 src/compass/utils/geo_metadata.py create mode 100644 src/compass/utils/h5_helpers.py diff --git a/src/compass/s1_geo2rdr.py b/src/compass/s1_geo2rdr.py index dd0e3190..63c50f73 100755 --- a/src/compass/s1_geo2rdr.py +++ b/src/compass/s1_geo2rdr.py @@ -58,13 +58,10 @@ def run(cfg: dict): # Run geo2rdr once per burst ID + date pair for burst in cfg.bursts: # Extract date string and create directory - burst_id = burst.burst_id + burst_id = str(burst.burst_id) date_str = burst.sensing_start.strftime("%Y%m%d") id_date = (burst_id, date_str) - - # Create top output path - top_output_path = f'{cfg.scratch_path}/{burst_id}' - os.makedirs(top_output_path, exist_ok=True) + out_paths = cfg.output_paths[id_date] # This ensures running geo2rdr only once; avoiding running for the different polarizations of the same burst_id if id_date in id_dates_processed: @@ -75,10 +72,6 @@ def run(cfg: dict): ref_burst_path = cfg.reference_radar_info.path topo_raster = isce3.io.Raster(f'{ref_burst_path}/topo.vrt') - # Create date directory - burst_output_path = f'{top_output_path}/{date_str}' - os.makedirs(burst_output_path, exist_ok=True) - # Get radar grid and orbit rdr_grid = burst.as_isce3_radargrid() orbit = burst.orbit @@ -90,7 +83,7 @@ def run(cfg: dict): blocksize) # Execute geo2rdr - geo2rdr_obj.geo2rdr(topo_raster, burst_output_path) + geo2rdr_obj.geo2rdr(topo_raster, out_paths.output_directory) dt = str(timedelta(seconds=time.time() - t_start)).split(".")[0] info_channel.log(f"{module_name} burst successfully ran in {dt} (hr:min:sec)") diff --git a/src/compass/s1_geocode_metadata.py b/src/compass/s1_geocode_metadata.py index 171fe05a..c07f424d 100755 --- a/src/compass/s1_geocode_metadata.py +++ b/src/compass/s1_geocode_metadata.py @@ -6,8 +6,10 @@ import os import time +import h5py import isce3 import journal +import numpy as np from osgeo import gdal from compass.utils.runconfig import RunConfig @@ -15,9 +17,9 @@ from compass.utils.yaml_argparse import YamlArgparse -def run(cfg, fetch_from_scratch=False): +def run(cfg, burst, fetch_from_scratch=False): ''' - Geocode metadata layers + Geocode metadata layers in single HDF5 Parameters ---------- @@ -25,7 +27,6 @@ def run(cfg, fetch_from_scratch=False): Dictionary with user runconfig options fetch_from_scratch: bool If True grabs metadata layers from scratch dir - ''' module_name = get_module_name(__file__) info_channel = journal.info(f"{module_name}.run") @@ -43,29 +44,24 @@ def run(cfg, fetch_from_scratch=False): threshold = cfg.geo2rdr_params.threshold iters = cfg.geo2rdr_params.numiter blocksize = cfg.geo2rdr_params.lines_per_block - output_epsg = cfg.geocoding_params.output_epsg output_format = cfg.geocoding_params.output_format # process one burst only - burst = cfg.bursts[0] date_str = burst.sensing_start.strftime("%Y%m%d") - burst_id = burst.burst_id + burst_id = str(burst.burst_id) geo_grid = cfg.geogrids[burst_id] - - os.makedirs(cfg.output_dir, exist_ok=True) - - scratch_path = f'{cfg.scratch_path}/{burst_id}/{date_str}' - os.makedirs(scratch_path, exist_ok=True) + output_epsg = geo_grid.epsg radar_grid = burst.as_isce3_radargrid() orbit = burst.orbit - # Initialize input/output path - input_path = f'{cfg.product_path}/{burst_id}/{date_str}' - output_path = input_path + # Initialize input/output paths + burst_id_date_key = (burst_id, date_str) + out_paths = cfg.output_paths[burst_id_date_key] + + input_path = out_paths.output_directory if fetch_from_scratch: - input_path = f'{cfg.scratch_path}/{burst_id}/{date_str}' - os.makedirs(output_path, exist_ok=True) + input_path = out_paths.scratch_directory # Initialize geocode object geo = isce3.geocode.GeocodeFloat32() @@ -79,46 +75,47 @@ def run(cfg, fetch_from_scratch=False): geo_grid.spacing_x, geo_grid.spacing_y, geo_grid.width, geo_grid.length, geo_grid.epsg) + # Geocode list of products + geotransform = [geo_grid.start_x, geo_grid.spacing_x, 0, + geo_grid.start_y, 0, geo_grid.spacing_y] + # Get the metadata layers to compute meta_layers = {'x': cfg.rdr2geo_params.compute_longitude, 'y': cfg.rdr2geo_params.compute_latitude, 'z': cfg.rdr2geo_params.compute_height, 'incidence': cfg.rdr2geo_params.compute_incidence_angle, - 'localIncidence': cfg.rdr2geo_params.compute_local_incidence_angle, - 'heading': cfg.rdr2geo_params.compute_azimuth_angle} - input_rasters = [ - isce3.io.Raster(f'{input_path}/{fname}.rdr') - for fname, enabled in meta_layers.items() if enabled] - output_rasters = [ - isce3.io.Raster(f'{output_path}/{fname}.geo', - geo_grid.width, geo_grid.length, 1, gdal.GDT_Float32, - output_format) - for fname, enabled in meta_layers.items() if enabled] - - # Geocode list of products - geotransform = [geo_grid.start_x, geo_grid.spacing_x, 0, - geo_grid.start_y, 0, geo_grid.spacing_y] - for input_raster, output_raster in zip(input_rasters, output_rasters): - geo.geocode(radar_grid=radar_grid, input_raster=input_raster, - output_raster=output_raster, dem_raster=dem_raster, - output_mode=isce3.geocode.GeocodeOutputMode.INTERP) - output_raster.set_geotransform(geotransform) - output_raster.set_epsg(output_epsg) - del output_raster - - # Geocode layover shadow separately - if cfg.rdr2geo_params.compute_layover_shadow_mask: - input_raster = isce3.io.Raster(f'{input_path}/layoverShadowMask.rdr') - output_raster = isce3.io.Raster(f'{output_path}/layoverShadowMask.geo', - geo_grid.width, geo_grid.length, 1, - gdal.GDT_Byte, output_format) - geo.data_interpolator = 'NEAREST' - geo.geocode(radar_grid=radar_grid, input_raster=input_raster, - output_raster=output_raster, dem_raster=dem_raster, - output_mode=isce3.geocode.GeocodeOutputMode.INTERP) - output_raster.set_geotransform(geotransform) - output_raster.set_epsg(output_epsg) - del output_raster + 'local_incidence': cfg.rdr2geo_params.compute_local_incidence_angle, + 'heading': cfg.rdr2geo_params.compute_azimuth_angle, + 'layover_shadow_mask': cfg.rdr2geo_params.compute_layover_shadow_mask} + + out_h5 = f'{out_paths.output_directory}/topo.h5' + shape = (geo_grid.length, geo_grid.width) + with h5py.File(out_h5, 'w') as topo_h5: + for layer_name, enabled in meta_layers.items(): + if not enabled: + continue + dtype = np.single + # layoverShadowMask is last option, no need to change data type + # and interpolator afterwards + if layer_name == 'layover_shadow_mask': + geo.data_interpolator = 'NEAREST' + dtype = np.byte + + topo_ds = topo_h5.create_dataset(layer_name, dtype=dtype, + shape=shape) + topo_ds.attrs['description'] = np.string_(layer_name) + output_raster = isce3.io.Raster(f"IH5:::ID={topo_ds.id.id}".encode("utf-8"), + update=True) + + input_raster = isce3.io.Raster(f'{input_path}/{layer_name}.rdr') + + geo.geocode(radar_grid=radar_grid, input_raster=input_raster, + output_raster=output_raster, dem_raster=dem_raster, + output_mode=isce3.geocode.GeocodeOutputMode.INTERP) + output_raster.set_geotransform(geotransform) + output_raster.set_epsg(output_epsg) + del input_raster + del output_raster dt = str(timedelta(seconds=time.time() - t_start)).split(".")[0] info_channel.log( diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py index 853d06b9..270247a1 100755 --- a/src/compass/s1_geocode_slc.py +++ b/src/compass/s1_geocode_slc.py @@ -3,25 +3,24 @@ '''wrapper for geocoded SLC''' from datetime import timedelta -import os import time +import h5py import isce3 import journal import numpy as np -from osgeo import gdal +from s1reader.s1_reader import is_eap_correction_necessary from compass import s1_rdr2geo from compass import s1_geocode_metadata - from compass.utils.elevation_antenna_pattern import apply_eap_correction -from compass.utils.geo_metadata import GeoCslcMetadata from compass.utils.geo_runconfig import GeoRunConfig +from compass.utils.h5_helpers import (init_geocoded_dataset, + geo_burst_metadata_to_hdf5) from compass.utils.helpers import get_module_name from compass.utils.lut import compute_geocoding_correction_luts from compass.utils.range_split_spectrum import range_split_spectrum from compass.utils.yaml_argparse import YamlArgparse -from s1reader.s1_reader import is_eap_correction_necessary def run(cfg: GeoRunConfig): ''' @@ -47,7 +46,6 @@ def run(cfg: GeoRunConfig): blocksize = cfg.geo2rdr_params.lines_per_block flatten = cfg.geocoding_params.flatten - # process one burst only for burst in cfg.bursts: # Reinitialize the dem raster per burst to prevent raster artifacts # caused by modification in geocodeSlc @@ -59,21 +57,13 @@ def run(cfg: GeoRunConfig): date_str = burst.sensing_start.strftime("%Y%m%d") burst_id = str(burst.burst_id) pol = burst.polarization - id_pol = f"{burst_id}_{pol}" geo_grid = cfg.geogrids[burst_id] - # Create top output path - burst_output_path = f'{cfg.product_path}/{burst_id}/{date_str}' - os.makedirs(burst_output_path, exist_ok=True) - - scratch_path = f'{cfg.scratch_path}/{burst_id}/{date_str}' - os.makedirs(scratch_path, exist_ok=True) - - # Get range and azimuth LUTs rg_lut, az_lut = compute_geocoding_correction_luts(burst, rg_step=cfg.lut_params.range_spacing, az_step=cfg.lut_params.azimuth_spacing) + radar_grid = burst.as_isce3_radargrid() native_doppler = burst.doppler.lut2d orbit = burst.orbit @@ -85,13 +75,20 @@ def run(cfg: GeoRunConfig): if cfg.rdr2geo_params.enabled: s1_rdr2geo.run(cfg, save_in_scratch=True) if cfg.rdr2geo_params.geocode_metadata_layers: - s1_geocode_metadata.run(cfg, fetch_from_scratch=True) + s1_geocode_metadata.run(cfg, burst, fetch_from_scratch=True) + + # Get output paths for current burst + burst_id_date_key = (burst_id, date_str) + out_paths = cfg.output_paths[burst_id_date_key] + + # Create scratch as needed + scratch_path = out_paths.scratch_directory # Load the input burst SLC - temp_slc_path = f'{scratch_path}/{id_pol}_temp.vrt' + temp_slc_path = f'{scratch_path}/{out_paths.file_name_stem}_temp.vrt' burst.slc_to_vrt_file(temp_slc_path) - # Check if EAP correction is necessary + # Apply EAP correction if necessary check_eap = is_eap_correction_necessary(burst.ipf_version) if check_eap.phase_correction: temp_slc_path_corrected = temp_slc_path.replace('_temp.vrt', @@ -100,6 +97,7 @@ def run(cfg: GeoRunConfig): temp_slc_path, temp_slc_path_corrected, check_eap) + # Replace the input burst if the correction is applied temp_slc_path = temp_slc_path_corrected @@ -113,43 +111,55 @@ def run(cfg: GeoRunConfig): else: rdr_burst_raster = isce3.io.Raster(temp_slc_path) - # Generate output geocoded burst raster - geo_burst_raster = isce3.io.Raster( - f'{burst_output_path}/{id_pol}.slc', - geo_grid.width, geo_grid.length, - rdr_burst_raster.num_bands, gdal.GDT_CFloat32, - cfg.geocoding_params.output_format) - # Extract burst boundaries b_bounds = np.s_[burst.first_valid_line:burst.last_valid_line, - burst.first_valid_sample:burst.last_valid_sample] + burst.first_valid_sample:burst.last_valid_sample] # Create sliced radar grid representing valid region of the burst sliced_radar_grid = burst.as_isce3_radargrid()[b_bounds] - # Geocode - isce3.geocode.geocode_slc(geo_burst_raster, rdr_burst_raster, - dem_raster, - radar_grid, sliced_radar_grid, - geo_grid, orbit, - native_doppler, - image_grid_doppler, ellipsoid, threshold, - iters, blocksize, flatten, - azimuth_carrier=az_carrier_poly2d) - - # Set geo transformation - geotransform = [geo_grid.start_x, geo_grid.spacing_x, 0, - geo_grid.start_y, 0, geo_grid.spacing_y] - geo_burst_raster.set_geotransform(geotransform) - geo_burst_raster.set_epsg(geo_grid.epsg) - del geo_burst_raster - del dem_raster # modified in geocodeSlc - - # Save burst metadata - metadata = GeoCslcMetadata.from_georunconfig(cfg, burst_id) - json_path = f'{burst_output_path}/{id_pol}.json' - with open(json_path, 'w') as f_json: - metadata.to_file(f_json, 'json') + output_hdf5 = out_paths.hdf5_path + with h5py.File(output_hdf5, 'w') as geo_burst_h5: + geo_burst_h5.attrs['Conventions'] = "CF-1.8" + + # add type to root for GDAL recognition of datasets + ctype = h5py.h5t.py_create(np.complex64) + ctype.commit(geo_burst_h5['/'].id, np.string_('complex64')) + + backscatter_group = geo_burst_h5.require_group('complex_backscatter') + init_geocoded_dataset(backscatter_group, pol, geo_grid, + 'complex64', f'{pol} geocoded SLC image') + + # access the HDF5 dataset for a given frequency and polarization + dataset_path = f'/complex_backscatter/{pol}' + gslc_dataset = geo_burst_h5[dataset_path] + + # Construct the output raster directly from HDF5 dataset + geo_burst_raster = isce3.io.Raster(f"IH5:::ID={gslc_dataset.id.id}".encode("utf-8"), update=True) + + # Geocode + isce3.geocode.geocode_slc(geo_burst_raster, rdr_burst_raster, + dem_raster, + radar_grid, sliced_radar_grid, + geo_grid, orbit, + native_doppler, + image_grid_doppler, ellipsoid, threshold, + iters, blocksize, flatten, + azimuth_carrier=az_carrier_poly2d) + + # Set geo transformation + geotransform = [geo_grid.start_x, geo_grid.spacing_x, 0, + geo_grid.start_y, 0, geo_grid.spacing_y] + geo_burst_raster.set_geotransform(geotransform) + geo_burst_raster.set_epsg(epsg) + del geo_burst_raster + del dem_raster # modified in geocodeSlc + + # Save burst metadata with new h5py File instance because io.Raster things + with h5py.File(output_hdf5, 'a') as geo_burst_h5: + geo_burst_metadata_to_hdf5(geo_burst_h5, burst, geo_grid, cfg) + geo_burst_h5['metadata/runconfig'] = cfg.yaml_string + dt = str(timedelta(seconds=time.time() - t_start)).split(".")[0] info_channel.log(f"{module_name} burst successfully ran in {dt} (hr:min:sec)") diff --git a/src/compass/s1_geocode_stack.py b/src/compass/s1_geocode_stack.py index 1627634d..08ee86d5 100755 --- a/src/compass/s1_geocode_stack.py +++ b/src/compass/s1_geocode_stack.py @@ -122,17 +122,10 @@ def generate_burst_map(zip_files, orbit_dir, output_epsg=None, bbox=None, if epsg is None: # Flag for skipping burst continue - burst_map['burst_id'].append(burst.burst_id) + burst_map['burst_id'].append(str(burst.burst_id)) # keep the burst object so we don't have to re-parse burst_map['burst'].append(burst) - left, bottom, right, top = bbox_utm - burst_map['x_top_left'].append(left) - burst_map['y_top_left'].append(top) - burst_map['x_bottom_right'].append(right) - burst_map['y_bottom_right'].append(bottom) - - burst_map['epsg'].append(epsg) burst_map['date'].append(burst.sensing_start.strftime("%Y%m%d")) # Save the file paths for creating the runconfig burst_map['orbit_path'].append(orbit_path) @@ -151,7 +144,7 @@ def _get_burst_epsg_and_bbox(burst, output_epsg, bbox, bbox_epsg, burst_db_file) if output_epsg is None: if os.path.exists(burst_db_file): epsg, _ = helpers.burst_bbox_from_db( - burst.burst_id, burst_db_file + str(burst.burst_id), burst_db_file ) else: # Fallback: ust the burst center UTM zone @@ -172,7 +165,7 @@ def _get_burst_epsg_and_bbox(burst, output_epsg, bbox, bbox_epsg, burst_db_file) return None, None else: epsg_db, bbox_utm = helpers.burst_bbox_from_db( - burst.burst_id, burst_db_file + str(burst.burst_id), burst_db_file ) if epsg_db != epsg: bbox_utm = helpers.transform_bbox( @@ -232,7 +225,8 @@ def get_common_burst_ids(data): def create_runconfig(burst_map_row, dem_file, work_dir, flatten, enable_rss, - low_band, high_band, pol, x_spac, y_spac, enable_metadata): + low_band, high_band, pol, x_spac, y_spac, enable_metadata, + burst_db_file): """ Create runconfig to process geocoded bursts @@ -282,8 +276,9 @@ def create_runconfig(burst_map_row, dem_file, work_dir, flatten, enable_rss, burst = burst_map_row.burst inputs['safe_file_path'] = [burst_map_row.zip_file] inputs['orbit_file_path'] = [burst_map_row.orbit_path] - inputs['burst_id'] = burst.burst_id + inputs['burst_id'] = [str(burst.burst_id)] groups['dynamic_ancillary_file_group']['dem_file'] = dem_file + groups['dynamic_ancillary_file_group']['burst_database_file'] = burst_db_file # Product path product['product_path'] = work_dir @@ -296,14 +291,6 @@ def create_runconfig(burst_map_row, dem_file, work_dir, flatten, enable_rss, geocode['x_posting'] = x_spac geocode['y_posting'] = y_spac - geocode['top_left']['x'] = burst_map_row.x_top_left - geocode['top_left']['y'] = burst_map_row.y_top_left - geocode['bottom_right']['x'] = burst_map_row.x_bottom_right - geocode['bottom_right']['y'] = burst_map_row.y_bottom_right - # geocode['x_snap'] = None - # geocode['y_snap'] = None - geocode['output_epsg'] = burst_map_row.epsg - # Range split spectrum rss['enabled'] = enable_rss rss['low_band_bandwidth'] = low_band @@ -315,7 +302,7 @@ def create_runconfig(burst_map_row, dem_file, work_dir, flatten, enable_rss, date_str = burst.sensing_start.strftime("%Y%m%d") os.makedirs(f'{work_dir}/runconfigs', exist_ok=True) - runconfig_path = f'{work_dir}/runconfigs/geo_runconfig_{date_str}_{burst.burst_id}.yaml' + runconfig_path = f'{work_dir}/runconfigs/geo_runconfig_{date_str}_{str(burst.burst_id)}.yaml' with open(runconfig_path, 'w') as yaml_file: yaml.dump(yaml_cfg, yaml_file, default_flow_style=False) return runconfig_path @@ -484,6 +471,7 @@ def run(slc_dir, dem_file, burst_id, start_date=None, end_date=None, exclude_dat x_spac, y_spac, do_metadata, + burst_db_file=burst_db_file, ) date_str = row.burst.sensing_start.strftime("%Y%m%d") runfile_name = f'{run_dir}/run_{date_str}_{row.burst.burst_id}.sh' diff --git a/src/compass/s1_rdr2geo.py b/src/compass/s1_rdr2geo.py index d609e6bd..298f593e 100755 --- a/src/compass/s1_rdr2geo.py +++ b/src/compass/s1_rdr2geo.py @@ -60,14 +60,15 @@ def run(cfg, save_in_scratch=False): for burst in cfg.bursts: # extract date string and create directory date_str = burst.sensing_start.strftime("%Y%m%d") - burst_id = burst.burst_id + burst_id = str(burst.burst_id) pol = burst.polarization # init output directory in product_path - output_path = f'{cfg.product_path}/{burst_id}/{date_str}' + burst_id_date_key = (burst_id, date_str) + out_paths = cfg.output_paths[burst_id_date_key] + output_path = out_paths.output_directory if save_in_scratch: - output_path = f'{cfg.scratch_path}/{burst_id}/{date_str}' - os.makedirs(output_path, exist_ok=True) + output_path = out_paths.scratch_directory # save SLC to ENVI for all bursts # run rdr2geo for only 1 burst avoid redundancy @@ -101,12 +102,12 @@ def run(cfg, save_in_scratch=False): topo_output = {'x': (rdr2geo_cfg.compute_longitude, gdal.GDT_Float64), 'y': (rdr2geo_cfg.compute_latitude, gdal.GDT_Float64), 'z': (rdr2geo_cfg.compute_height, gdal.GDT_Float64), - 'layoverShadowMask': ( + 'layover_shadow_mask': ( cfg.rdr2geo_params.compute_layover_shadow_mask, gdal.GDT_Byte), 'incidence': (rdr2geo_cfg.compute_incidence_angle, gdal.GDT_Float32), - 'localIncidence': ( + 'local_incidence': ( rdr2geo_cfg.compute_local_incidence_angle, gdal.GDT_Float32), 'heading': ( diff --git a/src/compass/s1_resample.py b/src/compass/s1_resample.py index b31ad29e..1cb5a25b 100755 --- a/src/compass/s1_resample.py +++ b/src/compass/s1_resample.py @@ -46,21 +46,17 @@ def run(cfg: dict): # Process all bursts for burst in cfg.bursts: - # get burst ID of current burst - burst_id = burst.burst_id + # get burst ID and date string of current burst + burst_id = str(burst.burst_id) + date_str = burst.sensing_start.strftime("%Y%m%d") # Create top output path - top_output_path = f'{cfg.product_path}/{burst_id}' - os.makedirs(top_output_path, exist_ok=True) + burst_id_date_key = (burst_id, date_str) + out_paths = cfg.output_paths[burst_id_date_key] # Get reference burst radar grid ref_rdr_grid = cfg.reference_radar_info.grid - # Extract date string and create directory - date_str = burst.sensing_start.strftime("%Y%m%d") - burst_output_path = f'{top_output_path}/{date_str}' - os.makedirs(burst_output_path, exist_ok=True) - # Extract polarization pol = burst.polarization @@ -76,8 +72,7 @@ def run(cfg: dict): resamp_obj.lines_per_tile = blocksize # Get range and azimuth offsets - offset_path = f'{cfg.scratch_path}/' \ - f'{burst_id}/{date_str}' + offset_path = out_paths.scratch_directory rg_off_raster = isce3.io.Raster(f'{offset_path}/range.off') az_off_raster = isce3.io.Raster(f'{offset_path}/azimuth.off') @@ -87,7 +82,7 @@ def run(cfg: dict): original_raster = isce3.io.Raster(sec_burst_path) # Prepare resamled SLC as raster object - coreg_burst_path = f'{burst_output_path}/{burst_id}_{date_str}_{pol}.slc' + coreg_burst_path = f'{out_paths.output_directory}/{out_paths.file_name_stem}.slc' resampled_raster = isce3.io.Raster(coreg_burst_path, rg_off_raster.width, rg_off_raster.length, diff --git a/src/compass/utils/geo_grid.py b/src/compass/utils/geo_grid.py index 098d6657..f7629bab 100644 --- a/src/compass/utils/geo_grid.py +++ b/src/compass/utils/geo_grid.py @@ -387,7 +387,7 @@ def generate_geogrids(bursts, geo_dict, dem): geo_grids = {} for burst in bursts: - burst_id = burst.burst_id + burst_id = str(burst.burst_id) if burst_id in geo_grids: continue diff --git a/src/compass/utils/geo_metadata.py b/src/compass/utils/geo_metadata.py deleted file mode 100644 index 4febe510..00000000 --- a/src/compass/utils/geo_metadata.py +++ /dev/null @@ -1,278 +0,0 @@ -from dataclasses import dataclass -from datetime import datetime -import json -from types import SimpleNamespace - -import isce3 -from isce3.core import LUT2d, Poly1d, Orbit -from isce3.product import GeoGridParameters -import numpy as np -from ruamel.yaml import YAML -from shapely.geometry import Point, Polygon - -from compass.utils.geo_runconfig import GeoRunConfig -from compass.utils.raster_polygon import get_boundary_polygon -from compass.utils.wrap_namespace import wrap_namespace, unwrap_to_dict - -from s1reader.s1_reader import is_eap_correction_necessary - -def _poly1d_from_dict(poly1d_dict) -> Poly1d: - return Poly1d(poly1d_dict['coeffs'], poly1d_dict['mean'], - poly1d_dict['std']) - - -def _lut2d_from_dict(lut2d_dict) -> LUT2d: - lut2d_shape = (lut2d_dict['length'], lut2d_dict['width']) - lut2d_data = np.array(lut2d_dict['data']).reshape(lut2d_shape) - return LUT2d(lut2d_dict['x_start'], lut2d_dict['y_start'], - lut2d_dict['x_spacing'], lut2d_dict['y_spacing'], - lut2d_data) - - -def _orbit_from_dict(orbit_dict) -> Orbit: - ref_epoch = isce3.core.DateTime(orbit_dict['ref_epoch']) - - # build state vector - dt = float(orbit_dict['time']['spacing']) - t0 = ref_epoch + isce3.core.TimeDelta(float(orbit_dict['time']['first'])) - n_pts = int(orbit_dict['time']['size']) - orbit_sv = [[]] * n_pts - for i in range(n_pts): - t = t0 + isce3.core.TimeDelta(i * dt) - pos = [float(orbit_dict[f'position_{xyz}'][i]) for xyz in 'xyz'] - vel = [float(orbit_dict[f'velocity_{xyz}'][i]) for xyz in 'xyz'] - orbit_sv[i] = isce3.core.StateVector(t, pos, vel) - - return Orbit(orbit_sv, ref_epoch) - - -@dataclass(frozen=True) -class GeoCslcMetadata(): - # subset of burst class attributes - sensing_start: datetime - sensing_stop: datetime - radar_center_frequency: float - wavelength: float - azimuth_steer_rate: float - azimuth_time_interval: float - slant_range_time: float - starting_range: float - range_sampling_rate: float - range_pixel_spacing: float - azimuth_fm_rate: Poly1d - doppler: Poly1d - range_bandwidth: float - polarization: str # {VV, VH, HH, HV} - burst_id: str # t{track_number}_{burst_number}_iw{1,2,3} - platform_id: str # S1{A,B} - center: Point # {center lon, center lat} in degrees - border: Polygon # list of lon, lat coordinate tuples (in degrees) representing burst border - orbit: isce3.core.Orbit - orbit_direction: str - # VRT params - tiff_path: str # path to measurement tiff in SAFE/zip - i_burst: int - # window parameters - range_window_type: str - range_window_coefficient: float - - runconfig: SimpleNamespace - geogrid: GeoGridParameters - nodata: str - input_data_ipf_version: str - isce3_version: str - - # correction applied - eap_correction_by_opera: bool - - @classmethod - def from_georunconfig(cls, cfg: GeoRunConfig, burst_id: str): - '''Create GeoBurstMetadata class from GeoRunConfig object - - Parameter: - --------- - cfg : GeoRunConfig - GeoRunConfig containing geocoded burst metadata - burst_id : str - ID of burst to create metadata object for - ''' - burst = None - for b in cfg.bursts: - if b.burst_id == burst_id: - burst = b - - if burst is None: - err_str = f'{burst_id} not found in cfg.bursts' - raise ValueError(err_str) - - geogrid = cfg.geogrids[burst_id] - - # get boundary from geocoded raster - date_str = burst.sensing_start.strftime("%Y%m%d") - pol = burst.polarization - burst_output_path = f'{cfg.product_path}/{burst_id}/{date_str}' - geo_raster_path = f'{burst_output_path}/{burst_id}_{pol}.slc' - geo_boundary = get_boundary_polygon(geo_raster_path, np.nan) - center = geo_boundary.centroid - - # place holders - nodata_val = '?' - ipf_ver = '?' - isce3_ver = '?' - - # correction applied - check_eap = is_eap_correction_necessary(burst.ipf_version) - eap_correction_applied = True if check_eap.phase_correction else False - - return cls(burst.sensing_start, burst.sensing_stop, - burst.radar_center_frequency, burst.wavelength, - burst.azimuth_steer_rate, burst.azimuth_time_interval, - burst.slant_range_time, burst.starting_range, - burst.range_sampling_rate, burst.range_pixel_spacing, - burst.azimuth_fm_rate, burst.doppler.poly1d, - burst.range_bandwidth, burst.polarization, burst_id, - burst.platform_id, center, geo_boundary, burst.orbit, - burst.orbit_direction, burst.tiff_path, burst.i_burst, - burst.range_window_type, burst.range_window_coefficient, - cfg.groups, geogrid, nodata_val, ipf_ver, isce3_ver, - eap_correction_applied) - - - @classmethod - def from_file(cls, file_path: str, fmt: str): - '''Create GeoBurstMetadata class from json file - - Parameter: - --------- - file_path: str - File containing geocoded burst metadata - ''' - if fmt == 'yaml': - yaml = YAML(typ='safe') - load = yaml.load - elif fmt == 'json': - load = json.load - else: - raise ValueError(f'{fmt} unsupported. Only "json" or "yaml" supported') - - with open(file_path, 'r') as fid: - meta_dict = load(fid) - - datetime_fmt = "%Y-%m-%d %H:%M:%S.%f" - sensing_start = datetime.strptime(meta_dict['sensing_start'], - datetime_fmt) - sensing_stop = datetime.strptime(meta_dict['sensing_stop'], - datetime_fmt) - - azimuth_fm_rate = _poly1d_from_dict(meta_dict['azimuth_fm_rate']) - - dopp_poly1d = _poly1d_from_dict(meta_dict['doppler']) - - orbit = _orbit_from_dict(meta_dict['orbit']) - - # init geo_runconfig - cfg = wrap_namespace(meta_dict['runconfig']) - - # init geogrid - grid_dict = meta_dict['geogrid'] - geogrid = GeoGridParameters(grid_dict['start_x'], grid_dict['start_y'], - grid_dict['spacing_x'], - grid_dict['spacing_y'], - grid_dict['length'], grid_dict['width'], - grid_dict['epsg']) - - # get boundary from geocoded raster - product_path = cfg.product_path_group.product_path - date_str = sensing_start.strftime("%Y%m%d") - burst_id = meta_dict['burst_id'] - pol = meta_dict['polarization'] - output_dir = f'{product_path}/{burst_id}/{date_str}' - file_stem = f'geo_{burst_id}_{pol}' - geo_raster_path = f'{output_dir}/{file_stem}' - geo_boundary = get_boundary_polygon(geo_raster_path, np.nan) - center = geo_boundary.centroid - - return cls(sensing_start, sensing_stop, - meta_dict['radar_center_frequency'], - meta_dict['wavelength'], meta_dict['azimuth_steer_rate'], - meta_dict['azimuth_time_interval'], - meta_dict['slant_range_time'], meta_dict['starting_range'], - meta_dict['range_sampling_rate'], - meta_dict['range_pixel_spacing'], azimuth_fm_rate, - dopp_poly1d, meta_dict['range_bandwidth'], pol, - meta_dict['burst_id'], meta_dict['platform_id'], - center, geo_boundary, orbit, meta_dict['orbit_direction'], - meta_dict['tiff_path'], meta_dict['i_burst'], - meta_dict['range_window_type'], - meta_dict['range_window_coefficient'], cfg, geogrid, - meta_dict['nodata'], meta_dict['input_data_ipf_version'], - meta_dict['isce3_version'], - meta_dict['opera_corrected_antenna_pattern']) - - def as_dict(self): - ''' Convert self to dict for write to YAML/JSON - ''' - self_as_dict = {} - for key, val in self.__dict__.items(): - if key in ['border', 'center', 'sensing_start', 'sensing_stop']: - val = str(val) - elif isinstance(val, np.float64): - val = float(val) - elif key in ['azimuth_fm_rate', 'doppler']: - temp = {} - temp['order'] = val.order - temp['mean'] = val.mean - temp['std'] = val.std - temp['coeffs'] = val.coeffs - val = temp - elif key == 'orbit': - temp = {} - temp['ref_epoch'] = str(val.reference_epoch) - temp['time'] = {} - temp['time']['first'] = val.time.first - temp['time']['spacing'] = val.time.spacing - temp['time']['last'] = val.time.last - temp['time']['size'] = val.time.size - temp['position_x'] = val.position[:,0].tolist() - temp['position_y'] = val.position[:,1].tolist() - temp['position_z'] = val.position[:,2].tolist() - temp['velocity_x'] = val.velocity[:,0].tolist() - temp['velocity_y'] = val.velocity[:,1].tolist() - temp['velocity_z'] = val.velocity[:,2].tolist() - val = temp - elif key == 'runconfig': - val = unwrap_to_dict(val) - elif key == 'geogrid': - temp = {} - temp['start_x'] = val.start_x - temp['start_y'] = val.start_y - temp['spacing_x'] = val.spacing_x - temp['spacing_y'] = val.spacing_y - temp['length'] = val.length - temp['width'] = val.width - temp['epsg'] = val.epsg - val = temp - - self_as_dict[key] = val - - return self_as_dict - - def to_file(self, dst, fmt:str): - '''Write self to file - - Parameter: - --------- - dst: file pointer - File object to write metadata to - fmt: ['yaml', 'json'] - Format of output - ''' - self_as_dict = self.as_dict() - - if fmt == 'yaml': - yaml = YAML(typ='safe') - yaml.dump(self_as_dict, dst) - elif fmt == 'json': - json.dump(self_as_dict, dst, indent=4) - else: - raise ValueError(f'{fmt} unsupported. Only "json" or "yaml" supported') diff --git a/src/compass/utils/geo_runconfig.py b/src/compass/utils/geo_runconfig.py index 410ac5a4..74ef8cee 100644 --- a/src/compass/utils/geo_runconfig.py +++ b/src/compass/utils/geo_runconfig.py @@ -10,6 +10,7 @@ from compass.utils.geo_grid import (generate_geogrids_from_db, generate_geogrids, geogrid_as_dict) from compass.utils.runconfig import ( + create_output_paths, runconfig_to_bursts, load_validate_yaml, RunConfig) @@ -80,8 +81,14 @@ def load_from_yaml(cls, yaml_path: str, workflow_name: str) -> GeoRunConfig: # Empty reference dict for base runconfig class constructor empty_ref_dict = {} + with open(yaml_path, 'r') as f_yaml: + entire_yaml = f_yaml.read() + + # Get scratch and output paths + output_paths = create_output_paths(sns, bursts) + return cls(cfg['runconfig']['name'], sns, bursts, empty_ref_dict, - geogrids) + entire_yaml, output_paths, geogrids) @property def geocoding_params(self) -> dict: @@ -91,18 +98,6 @@ def geocoding_params(self) -> dict: def rdr2geo_params(self) -> dict: return self.groups.processing.rdr2geo - @property - def burst_id(self) -> str: - return self.bursts[0].burst_id - - @property - def sensing_start(self): - return self.bursts[0].sensing_start - - @property - def polarization(self) -> str: - return self.bursts[0].polarization - @property def split_spectrum_params(self) -> dict: return self.groups.processing.range_split_spectrum @@ -111,19 +106,6 @@ def split_spectrum_params(self) -> dict: def lut_params(self) -> dict: return self.groups.processing.correction_luts - @property - def output_dir(self) -> str: - date_str = self.sensing_start.strftime("%Y%m%d") - burst_id = self.burst_id - return f'{super().product_path}/{burst_id}/{date_str}' - - @property - def file_stem(self) -> str: - burst_id = self.burst_id - pol = self.polarization - return f'geo_{burst_id}_{pol}' - - def as_dict(self): ''' Convert self to dict for write to YAML/JSON diff --git a/src/compass/utils/h5_helpers.py b/src/compass/utils/h5_helpers.py new file mode 100644 index 00000000..4896bf54 --- /dev/null +++ b/src/compass/utils/h5_helpers.py @@ -0,0 +1,448 @@ +from datetime import datetime + +import isce3 +import numpy as np +from osgeo import osr + +import compass +from compass.utils.raster_polygon import get_boundary_polygon + +def init_geocoded_dataset(geocoded_group, dataset_name, geo_grid, dtype, + description): + ''' + Create and allocate dataset for isce.geocode.geocode_slc to write to that + is NC compliant + + geocoded_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 + description: str + Description of dataset to be geocoded + ''' + shape = (geo_grid.length, geo_grid.width) + geocoded_ds = geocoded_group.require_dataset(dataset_name, dtype=dtype, + shape=shape) + + geocoded_ds.attrs['description'] = description + + # Compute x scale + dx = geo_grid.spacing_x + x0 = geo_grid.start_x + 0.5 * dx + xf = x0 + (geo_grid.width - 1) * dx + x_vect = np.linspace(x0, xf, geo_grid.width, dtype=np.float64) + + # Compute y scale + dy = geo_grid.spacing_y + y0 = geo_grid.start_y + 0.5 * dy + yf = y0 + (geo_grid.length - 1) * dy + y_vect = np.linspace(y0, yf, geo_grid.length, dtype=np.float64) + + # 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 = geocoded_group.require_dataset('x', dtype='float64', data=x_vect, + shape=x_vect.shape) + y_ds = geocoded_group.require_dataset('y', dtype='float64', data=y_vect, + shape=y_vect.shape) + + # Mapping of dimension scales to datasets is not done automatically in HDF5 + # We should label appropriate arrays as scales and attach them to datasets + # explicitly as show below. + x_ds.make_scale() + geocoded_ds.dims[1].attach_scale(x_ds) + y_ds.make_scale() + geocoded_ds.dims[0].attach_scale(y_ds) + + # Associate grid mapping with data - projection created later + geocoded_ds.attrs['grid_mapping'] = np.string_("projection") + + # Set up osr for wkt + srs = osr.SpatialReference() + srs.ImportFromEPSG(geo_grid.epsg) + + #Create a new single int dataset for projections + projection_ds = geocoded_group.require_dataset('projection', (), dtype='i') + projection_ds[()] = geo_grid.epsg + + # WGS84 ellipsoid + projection_ds.attrs['semi_major_axis'] = 6378137.0 + projection_ds.attrs['inverse_flattening'] = 298.257223563 + projection_ds.attrs['ellipsoid'] = np.string_("WGS84") + + # Additional fields + projection_ds.attrs['epsg_code'] = geo_grid.epsg + + # CF 1.7+ requires this attribute to be named "crs_wkt" + # spatial_ref is old GDAL way. Using that for testing only. + # For NISAR replace with "crs_wkt" + projection_ds.attrs['spatial_ref'] = np.string_(srs.ExportToWkt()) + + # Here we have handcoded the attributes for the different cases + # Recommended method is to use pyproj.CRS.to_cf() as shown above + # To get complete set of attributes. + + # Geodetic latitude / longitude + if geo_grid.epsg == 4326: + #Set up grid mapping + projection_ds.attrs['grid_mapping_name'] = np.string_('latitude_longitude') + projection_ds.attrs['longitude_of_prime_meridian'] = 0.0 + + #Setup units for x and y + x_ds.attrs['standard_name'] = np.string_("longitude") + x_ds.attrs['units'] = np.string_("degrees_east") + + y_ds.attrs['standard_name'] = np.string_("latitude") + y_ds.attrs['units'] = np.string_("degrees_north") + + # UTM zones + elif (geo_grid.epsg > 32600 and geo_grid.epsg < 32661) or \ + (geo_grid.epsg > 32700 and geo_grid.epsg < 32761): + #Set up grid mapping + projection_ds.attrs['grid_mapping_name'] = np.string_('universal_transverse_mercator') + projection_ds.attrs['utm_zone_number'] = geo_grid.epsg % 100 + + #Setup units for x and y + x_ds.attrs['standard_name'] = np.string_("projection_x_coordinate") + x_ds.attrs['long_name'] = np.string_("x coordinate of projection") + x_ds.attrs['units'] = np.string_("m") + + y_ds.attrs['standard_name'] = np.string_("projection_y_coordinate") + y_ds.attrs['long_name'] = np.string_("y coordinate of projection") + y_ds.attrs['units'] = np.string_("m") + + # Polar Stereo North + elif geo_grid.epsg == 3413: + #Set up grid mapping + projection_ds.attrs['grid_mapping_name'] = np.string_("polar_stereographic") + projection_ds.attrs['latitude_of_projection_origin'] = 90.0 + projection_ds.attrs['standard_parallel'] = 70.0 + projection_ds.attrs['straight_vertical_longitude_from_pole'] = -45.0 + projection_ds.attrs['false_easting'] = 0.0 + projection_ds.attrs['false_northing'] = 0.0 + + #Setup units for x and y + x_ds.attrs['standard_name'] = np.string_("projection_x_coordinate") + x_ds.attrs['long_name'] = np.string_("x coordinate of projection") + x_ds.attrs['units'] = np.string_("m") + + y_ds.attrs['standard_name'] = np.string_("projection_y_coordinate") + y_ds.attrs['long_name'] = np.string_("y coordinate of projection") + y_ds.attrs['units'] = np.string_("m") + + # Polar Stereo south + elif geo_grid.epsg == 3031: + #Set up grid mapping + projection_ds.attrs['grid_mapping_name'] = np.string_("polar_stereographic") + projection_ds.attrs['latitude_of_projection_origin'] = -90.0 + projection_ds.attrs['standard_parallel'] = -71.0 + projection_ds.attrs['straight_vertical_longitude_from_pole'] = 0.0 + projection_ds.attrs['false_easting'] = 0.0 + projection_ds.attrs['false_northing'] = 0.0 + + #Setup units for x and y + x_ds.attrs['standard_name'] = np.string_("projection_x_coordinate") + x_ds.attrs['long_name'] = np.string_("x coordinate of projection") + x_ds.attrs['units'] = np.string_("m") + + y_ds.attrs['standard_name'] = np.string_("projection_y_coordinate") + y_ds.attrs['long_name'] = np.string_("y coordinate of projection") + y_ds.attrs['units'] = np.string_("m") + + # EASE 2 for soil moisture L3 + elif geo_grid.epsg == 6933: + #Set up grid mapping + projection_ds.attrs['grid_mapping_name'] = np.string_("lambert_cylindrical_equal_area") + projection_ds.attrs['longitude_of_central_meridian'] = 0.0 + projection_ds.attrs['standard_parallel'] = 30.0 + projection_ds.attrs['false_easting'] = 0.0 + projection_ds.attrs['false_northing'] = 0.0 + + #Setup units for x and y + x_ds.attrs['standard_name'] = np.string_("projection_x_coordinate") + x_ds.attrs['long_name'] = np.string_("x coordinate of projection") + x_ds.attrs['units'] = np.string_("m") + + y_ds.attrs['standard_name'] = np.string_("projection_y_coordinate") + y_ds.attrs['long_name'] = np.string_("y coordinate of projection") + y_ds.attrs['units'] = np.string_("m") + + # Europe Equal Area for Deformation map (to be implemented in isce3) + elif geo_grid.epsg == 3035: + #Set up grid mapping + projection_ds.attrs['grid_mapping_name'] = np.string_("lambert_azimuthal_equal_area") + projection_ds.attrs['longitude_of_projection_origin']= 10.0 + projection_ds.attrs['latitude_of_projection_origin'] = 52.0 + projection_ds.attrs['standard_parallel'] = -71.0 + projection_ds.attrs['straight_vertical_longitude_from_pole'] = 0.0 + projection_ds.attrs['false_easting'] = 4321000.0 + projection_ds.attrs['false_northing'] = 3210000.0 + + #Setup units for x and y + x_ds.attrs['standard_name'] = np.string_("projection_x_coordinate") + x_ds.attrs['long_name'] = np.string_("x coordinate of projection") + x_ds.attrs['units'] = np.string_("m") + + y_ds.attrs['standard_name'] = np.string_("projection_y_coordinate") + y_ds.attrs['long_name'] = np.string_("y coordinate of projection") + y_ds.attrs['units'] = np.string_("m") + + else: + raise NotImplementedError('Waiting for implementation / Not supported in ISCE3') + + +def geo_burst_metadata_to_hdf5(dst_h5, burst, geogrid, cfg): + '''Write burst metadata to HDF5 + + Parameter: + --------- + dst_h5: h5py File + HDF5 file meta data will be written to + burst: Sentinel1BurstSlc + Burst + ''' + time_str_fmt = 'time_str_fmt' + metadata_group = dst_h5.require_group('metadata') + + def add_dataset_and_attrs(group, name, value, attr_dict): + '''Write isce3.core.Poly1d properties to hdf5 + + Parameters + ---------- + group: h5py.Group + h5py Group to store poly1d parameters in + name: str + Name of dataset to add + value: object + Value to be added + attr_dict: dict[str: object] + Dict with attribute name as key and some object as value + ''' + if name in group: + del group[name] + + group[name] = value + val_ds = group[name] + for key, val in attr_dict.items(): + val_ds.attrs[key] = val + + # product identification and processing information + id_proc_group = metadata_group.require_group('identifcation_and_processing') + add_dataset_and_attrs(id_proc_group, 'product_id', 'L2_CSLC_S1', {}) + add_dataset_and_attrs(id_proc_group, 'product_version', '?', {}) + add_dataset_and_attrs(id_proc_group, 'software_version', + compass.__version__, + {'description': 'COMPASS version used to generate the L2_CSLC_S1 product'}) + add_dataset_and_attrs(id_proc_group, 'isce3_version', + isce3.__version__, + {'description': 'ISCE3 version used to generate the L2_CSLC_S1 product'}) + add_dataset_and_attrs(id_proc_group, 'project', 'OPERA', {}) + add_dataset_and_attrs(id_proc_group, 'product_level', '2', {}) + add_dataset_and_attrs(id_proc_group, 'product_type', 'CSLC_S1', {}) + add_dataset_and_attrs(id_proc_group, 'processing_datetime', + datetime.now().strftime(time_str_fmt), + {'description': 'L2_CSLC_S1 product processing date and time', + 'format': 'YYYY-MM-DDD HH:MM:SS'}) + add_dataset_and_attrs(id_proc_group, 'spacecraft_name', + burst.platform_id, + {'description': 'Name of Sensor platform (e.g., S1-A/B)'}) + + # burst metadata + s1ab_group = metadata_group.require_group('s1ab_burst_metadata') + add_dataset_and_attrs(s1ab_group, 'sensing_start', + burst.sensing_start.strftime('time_str_fmt.%f'), + {'description': 'Sensing start time of the burst', + 'format': 'YYYY-MM-DD HH:MM:SS.6f'}) + add_dataset_and_attrs(s1ab_group, 'sensing_stop', + burst.sensing_stop.strftime('time_str_fmt.%f'), + {'description':'Sensing stop time of the burst', + 'format': 'YYYY-MM-DD HH:MM:SS.6f'}) + add_dataset_and_attrs(s1ab_group, 'radar_center_frequency', + burst.radar_center_frequency, + {'description':'Radar center frequency', + 'units':'Hz'}) + add_dataset_and_attrs(s1ab_group, 'wavelength', burst.wavelength, + {'description':'Wavelength of the transmitted signal', + 'units':'meters'}) + add_dataset_and_attrs(s1ab_group, 'azimuth_steer_rate', + burst.azimuth_steer_rate, + {'description':'Azimuth steering rate of IW and EW modes', + 'units':'degrees/second'}) + # TODO add input width and length + add_dataset_and_attrs(s1ab_group, 'azimuth_time_interval', + burst.azimuth_time_interval, + {'description':'Time spacing between azimuth lines of the burst', + 'units':'seconds'}) + add_dataset_and_attrs(s1ab_group, 'starting_range', + burst.starting_range, + {'description':'Slant range of the first sample of the input burst', + 'units':'meters'}) + # TODO do we need this? It's not in the specs. + # TODO add far_range? + add_dataset_and_attrs(s1ab_group, 'slant_range_time', + burst.slant_range_time, + {'description':'two-way slant range time of Doppler centroid frequency estimate', + 'units':'seconds'}) + add_dataset_and_attrs(s1ab_group, 'range_pixel_spacing', + burst.range_pixel_spacing, + {'description':'Pixel spacing between slant range samples in the input burst SLC', + 'units':'meters'}) + add_dataset_and_attrs(s1ab_group, 'range_bandwidth', + burst.range_bandwidth, + {'description':'Slant range bandwidth of the signal', + 'units':'Hz'}) + add_dataset_and_attrs(s1ab_group, 'polarization', + burst.polarization, + {'description': 'Polarization of the burst'}) + add_dataset_and_attrs(s1ab_group, 'platform_id', + burst.platform_id, + {'description': 'Sensor platform identification string (e.g., S1A or S1B)'}) + # window parameters + add_dataset_and_attrs(s1ab_group, 'range_window_type', + burst.range_window_type, + {'description': 'name of the weighting window type used during processing'}) + add_dataset_and_attrs(s1ab_group, 'range_window_coefficient', + burst.range_window_coefficient, + {'description': 'value of the weighting window coefficient used during processing'}) + + def poly1d_to_h5(group, poly1d_name, poly1d): + '''Write isce3.core.Poly1d properties to hdf5 + + Parameters + ---------- + group: h5py.Group + h5py Group to store poly1d parameters in + poly1d_name: str + Name of Poly1d whose parameters are to be stored + poly1d: isce3.core.Poly1d + Poly1d ojbect whose parameters are to be stored + ''' + poly1d_group = group.require_group(poly1d_name) + add_dataset_and_attrs(poly1d_group, 'order', poly1d.order, + {'description': 'order of the polynomial'}) + add_dataset_and_attrs(poly1d_group, 'mean', poly1d.mean, + {'description': 'mean of the polynomial'}) + add_dataset_and_attrs(poly1d_group, 'std', poly1d.std, + {'description': 'standard deviation of the polynomial'}) + add_dataset_and_attrs(poly1d_group, 'coeffs', poly1d.coeffs, + {'description': 'coefficients of the polynomial'}) + poly1d_to_h5(s1ab_group, 'azimuth_fm_rate', burst.azimuth_fm_rate) + poly1d_to_h5(s1ab_group, 'doppler', burst.doppler.poly1d) + + # EAP metadata only written if it exists + if burst.burst_eap is not None: + eap_group = s1ab_group.require_group('elevation_antenna_pattern_correction') + eap = burst.burst_eap + add_dataset_and_attrs(eap_group, 'sampling_frequency', + eap.freq_sampling, + {'description': 'range sampling frequency', + 'units': 'Hz'}) + add_dataset_and_attrs(eap_group, 'eta_start', + eap.eta_start.strftime(time_str_fmt), + {'description': '?', + 'format': 'YYYY-MM-DD HH:MM:SS.6f'}) + add_dataset_and_attrs(eap_group, 'tau_0', eap.tau_0, + {'description': 'slant range time', + 'units': 'seconds'}) + add_dataset_and_attrs(eap_group, 'tau_sub', eap.tau_sub, + {'description': 'slant range time', + 'units': 'seconds'}) + add_dataset_and_attrs(eap_group, 'theta_sub', eap.theta_sub, + {'description': 'elevation angle', + 'units': 'radians'}) + add_dataset_and_attrs(eap_group, 'azimuth_time', + eap.azimuth_time.strftime(time_str_fmt), + {'description': '?', + 'format': 'YYYY-MM-DD HH:MM:SS.6f'}) + add_dataset_and_attrs(eap_group, 'ascending_node_time', + eap.ascending_node_time.strftime(time_str_fmt), + {'description': '?', + 'format': 'YYYY-MM-DD HH:MM:SS.6f'}) + + # save orbit + orbit_group = metadata_group.require_group('orbit') + ref_epoch = burst.orbit.reference_epoch.isoformat().replace('T', ' ') + add_dataset_and_attrs(orbit_group, 'ref_epoch', ref_epoch, + {'description': 'Reference epoch of the state vectors', + 'format': 'YYYY-MM-DD HH:MM:SS.6f'}) + orbit_time_obj = burst.orbit.time + add_dataset_and_attrs(orbit_group, 'time', + np.linspace(orbit_time_obj.first, + orbit_time_obj.last, + orbit_time_obj.size), + {'description': 'Time of the orbit state vectors relative to the reference epoch', + 'units': 'seconds'}) + for i_ax, axis in enumerate('xyz'): + add_dataset_and_attrs(orbit_group, f'position_{axis}', + burst.orbit.position[:, i_ax], + {'description': f'Platform position along {axis}-direction', + 'units': 'meters'}) + add_dataset_and_attrs(orbit_group, f'velocity_{axis}', + burst.orbit.velocity[:, i_ax], + {'description': f'Platform velocity along {axis}-direction', + 'units': 'meters/second'}) + add_dataset_and_attrs(orbit_group, 'orbit_direction', + burst.orbit_direction, + {'description':'Direction of sensor orbit ephermerides (e.g., ascending, descending)'}) + + # input params + input_group = metadata_group.require_group('input') + add_dataset_and_attrs(input_group, 'burst_file_path', + burst.tiff_path, + {'description': 'Path to TIFF file within the SAFE file containing the burst'}) + add_dataset_and_attrs(input_group, 'input_data_ipf_version', + str(burst.ipf_version), + {'description': 'Version of Instrument Processing Facility used to generate SAFE file'}) + add_dataset_and_attrs(input_group, 'dem_file', cfg.dem, + {'description': 'Path to DEM file'}) + + # processing params + processing_group = metadata_group.require_group('processing') + add_dataset_and_attrs(processing_group, 'burst_id', + str(burst.burst_id), + {'description': 'Unique burst identification string (ESA convention)'}) + + dataset_path_template = f'HDF5:%FILE_PATH%://complex_backscatter/{burst.polarization}' + geo_boundary = get_boundary_polygon(processing_group.file.filename, np.nan, + dataset_path_template) + center = geo_boundary.centroid + center_lon_lat = np.array([val[0] for val in center.coords.xy]) + add_dataset_and_attrs(processing_group, 'center', + center_lon_lat, + {'description': 'Burst geographical center coordinates in the projection system used for processing', + 'units': 'meters'}) + # list of coordinate tuples (in degrees) representing burst border + border_x, border_y = geo_boundary.exterior.coords.xy + border_lon_lat = np.array([[lon, lat] for lon, lat in zip(border_x, + border_y)]) + add_dataset_and_attrs(processing_group, 'border', border_lon_lat, + {'description': 'X- and Y- coordinates of the polygon including valid L2_CSLC_S1 data', + 'units': 'meters'}) + add_dataset_and_attrs(processing_group, 'i_burst', burst.i_burst, + {'description': 'Index of the burst of interest relative to other bursts in the S1-A/B SAFE file'}) + add_dataset_and_attrs(processing_group, 'start_x', geogrid.start_x, + {'description': 'X-coordinate of the L2_CSLC_S1 starting point in the coordinate system selected for processing', + 'units': 'meters'}) + add_dataset_and_attrs(processing_group, 'start_y', geogrid.start_y, + {'description': 'Y-coordinate of the L2_CSLC_S1 starting point in the coordinate system selected for processing', + 'units': 'meters'}) + add_dataset_and_attrs(processing_group, 'x_posting', + geogrid.spacing_x, + {'description': 'Spacing between product pixels along the X-direction ', + 'units': 'meters'}) + add_dataset_and_attrs(processing_group, 'y_posting', + geogrid.spacing_y, + {'description': 'Spacing between product pixels along the Y-direction ', + 'units': 'meters'}) + add_dataset_and_attrs(processing_group, 'width', geogrid.width, + {'description': 'Number of samples in the L2_CSLC_S1 product'}) + add_dataset_and_attrs(processing_group, 'length', geogrid.length, + {'description': 'Number of lines in the L2_CSLC_S1 product'}) + add_dataset_and_attrs(processing_group, 'epsg', + geogrid.epsg, + {'description': 'EPSG code identifying the coordinate system used for processing'}) + add_dataset_and_attrs(processing_group, 'no_data_value', 'NaN', + {'description': 'Value used when no data present'}) diff --git a/src/compass/utils/range_split_spectrum.py b/src/compass/utils/range_split_spectrum.py index 6380aeb4..02db8ee0 100644 --- a/src/compass/utils/range_split_spectrum.py +++ b/src/compass/utils/range_split_spectrum.py @@ -51,7 +51,7 @@ def range_split_spectrum(burst, ''' length, width = burst.shape lines_per_block = cfg_split_spectrum.lines_per_block - burst_id_pol = f'{burst.burst_id}_{burst.polarization}' + burst_id_pol = f'{str(burst.burst_id)}_{burst.polarization}' # In ISCE3, we can use raised cosine to implement S1-A/B Hamming window_type = burst.range_window_type diff --git a/src/compass/utils/raster_polygon.py b/src/compass/utils/raster_polygon.py index 819c451a..b9680fa0 100644 --- a/src/compass/utils/raster_polygon.py +++ b/src/compass/utils/raster_polygon.py @@ -8,7 +8,8 @@ from shapely.geometry import MultiPoint -def get_boundary_polygon(filename, invalid_value): +def get_boundary_polygon(filename, invalid_value=np.nan, + dataset_path_template=None): ''' Get boundary polygon for raster in 'filename'. Polygon includes only valid pixels @@ -19,6 +20,9 @@ def get_boundary_polygon(filename, invalid_value): Path to GDAL friendly file containing geocoded raster invalid_value: np.nan or float Invalid data value for raster in 'filename' + dataset_path_template: str + Template string with place holder string, %FILE_PATH%, to be replace by + actual path containing dataset. Will only be used if not None. Returns -------- @@ -30,7 +34,14 @@ def get_boundary_polygon(filename, invalid_value): f'{filename} not found') # Optimize this with block-processing? - ds = gdal.Open(filename, gdal.GA_ReadOnly) + if dataset_path_template is not None: + dataset_path = dataset_path_template.replace('%FILE_PATH%', filename) + else: + dataset_path = filename + try: + ds = gdal.Open(dataset_path, gdal.GA_ReadOnly) + except: + raise ValueError(f'GDAL unable to open: {dataset_path}') burst = ds.GetRasterBand(1).ReadAsArray() if np.isnan(invalid_value): diff --git a/src/compass/utils/runconfig.py b/src/compass/utils/runconfig.py index ed1c7e5f..af036e39 100644 --- a/src/compass/utils/runconfig.py +++ b/src/compass/utils/runconfig.py @@ -8,13 +8,13 @@ import journal import yamale from ruamel.yaml import YAML +from s1reader.s1_burst_slc import Sentinel1BurstSlc +from s1reader.s1_orbit import get_orbit_file_from_dir +from s1reader.s1_reader import load_bursts from compass.utils import helpers from compass.utils.radar_grid import file_to_rdr_grid from compass.utils.wrap_namespace import wrap_namespace, unwrap_to_dict -from s1reader.s1_burst_slc import Sentinel1BurstSlc -from s1reader.s1_orbit import get_orbit_file_from_dir -from s1reader.s1_reader import load_bursts def load_validate_yaml(yaml_path: str, workflow_name: str) -> dict: @@ -140,7 +140,7 @@ def validate_group_dict(group_cfg: dict, workflow_name) -> None: helpers.check_write_dir(product_path_group['sas_output_file']) -def runconfig_to_bursts(cfg: SimpleNamespace, auto_download: bool = False) -> list[Sentinel1BurstSlc]: +def runconfig_to_bursts(cfg: SimpleNamespace) -> list[Sentinel1BurstSlc]: '''Return bursts based on parameters in given runconfig Parameters @@ -201,7 +201,7 @@ def runconfig_to_bursts(cfg: SimpleNamespace, auto_download: bool = False) -> li # loop over burst objs extracted from SAFE zip for burst in load_bursts(safe_file, orbit_path, i_subswath, pol): # get burst ID - burst_id = burst.burst_id + burst_id = str(burst.burst_id) # include ALL bursts if no burst IDs given # is burst_id wanted? skip if not given in config @@ -242,15 +242,13 @@ def runconfig_to_bursts(cfg: SimpleNamespace, auto_download: bool = False) -> li return bursts -def get_ref_radar_grid_info(ref_path, burst_id): +def get_ref_radar_grid_info(ref_path): ''' Find all reference radar grids info Parameters ---------- ref_path: str Path where reference radar grids processing is stored - burst_id: str - Burst IDs for reference radar grids Returns ------- @@ -275,6 +273,35 @@ class ReferenceRadarInfo: grid: isce3.product.RadarGridParameters +def create_output_paths(sns, bursts): + # Generate scratch and output paths + output_paths = {} + product_paths = sns.product_path_group + for burst in bursts: + # Get burst ID and check if it already + burst_id = str(burst.burst_id) + date_str = burst.sensing_start.strftime("%Y%m%d") + + # Key for current burst ID + date combo + path_key = (burst_id, date_str) + + # Save output dir, output hdf5 and scratch dir to dict as + # SimpleNamespace + out_dir = f'{product_paths.product_path}/{burst_id}/{date_str}' + os.makedirs(out_dir, exist_ok=True) + + fname_stem = f"{burst_id}_{date_str}_{burst.polarization}" + h5_path = f"{out_dir}/{fname_stem}.h5" + + scratch_path = f'{product_paths.scratch_path}/{burst_id}/{date_str}' + os.makedirs(scratch_path, exist_ok=True) + + output_paths[path_key] = SimpleNamespace(output_directory=out_dir, + file_name_stem=fname_stem, + hdf5_path=h5_path, + scratch_directory=scratch_path) + return output_paths + @dataclass(frozen=True) class RunConfig: '''dataclass containing CSLC runconfig''' @@ -287,6 +314,12 @@ class RunConfig: # dict of reference radar paths and grids values keyed on burst ID # (empty/unused if rdr2geo) reference_radar_info: ReferenceRadarInfo + # entirety of yaml as string + yaml_string: str + # output paths: + # key = tuple[burst ID, date str] + # val = SimpleNamespace output directory path, HDF5 name, scratch directory path + output_paths: dict[tuple[str, str], SimpleNamespace] @classmethod def load_from_yaml(cls, yaml_path: str, workflow_name: str) -> RunConfig: @@ -313,7 +346,15 @@ def load_from_yaml(cls, yaml_path: str, workflow_name: str) -> RunConfig: sns.input_file_group.reference_burst.file_path, sns.input_file_group.burst_id) - return cls(cfg['runconfig']['name'], sns, bursts, ref_rdr_grid_info) + # For saving entire file as string to metadata. Stop gap for writing + # dict to individual elements to HDF5 metadata + with open(yaml_path, 'r') as f_yaml: + entire_yaml = f_yaml.read() + + output_paths = create_output_paths(sns, bursts) + + return cls(cfg['runconfig']['name'], sns, bursts, ref_rdr_grid_info, + entire_yaml, output_paths) @property def burst_id(self) -> list[str]: @@ -384,28 +425,33 @@ def gpu_id(self): return self.groups.worker.gpu_id def as_dict(self): - ''' Convert self to dict for write to YAML/JSON + '''Convert self to dict for write to YAML/JSON Unable to dataclasses.asdict() because isce3 objects can not be pickled ''' + # Convenience functions + def date_str(burst): + '''Burst datetime sensing_start to str conversion + ''' + return burst.sensing_start.date().strftime('%Y%m%d') + + def burst_as_key(burst): + '''Create an unique key of burst ID, date string, and polarization + ''' + return '_'.join([str(burst.burst_id), date_str(burst), burst.polarization]) + self_as_dict = {} for key, val in self.__dict__.items(): if key == 'groups': val = unwrap_to_dict(val) elif key == 'bursts': - # just date in datetime obj as string - date_str = lambda b : b.sensing_start.date().strftime('%Y%m%d') - - # create an unique burst key - burst_as_key = lambda b : '_'.join([b.burst_id, - date_str(b), - b.polarization]) - val = {burst_as_key(burst): burst.as_dict() for burst in val} self_as_dict[key] = val return self_as_dict def to_yaml(self): + '''Dump runconfig as string to sys.stdout + ''' self_as_dict = self.as_dict() yaml = YAML(typ='safe') - yaml.dump(self_as_dict, sys.stdout, indent=4) + yaml.dump(self_as_dict, sys.stdout) diff --git a/src/compass/utils/validate_cslc.py b/src/compass/utils/validate_cslc.py index 66be6459..ae803feb 100644 --- a/src/compass/utils/validate_cslc.py +++ b/src/compass/utils/validate_cslc.py @@ -1,7 +1,7 @@ import argparse -import json import os +import h5py import numpy as np from osgeo import gdal @@ -10,7 +10,6 @@ def cmd_line_parser(): """ Command line parser """ - parser = argparse.ArgumentParser(description="""Validate reference and generated (secondary) S1 CSLC products""", formatter_class=argparse.ArgumentDefaultsHelpFormatter) @@ -25,6 +24,44 @@ def cmd_line_parser(): return parser.parse_args() +def _gdal_nfo_retrieve(path_h5): + """ + Extract polarization, geotransform array, projection, and SLC from + given HDF5 + + Parameters + ---------- + path_h5: str + File path to CSLC HDF5 product + + Returns + ------- + pol: str + Polarization of CSLC product + geotransform: np.array + Array holding x/y start, x/y spacing, and geogrid dimensions + proj: str + EPSG projection of geogrid + slc: np.array + Array holding geocoded complex backscatter + """ + raster_key = 'complex_backscatter' + + # Extract polarization with h5py + with h5py.File(path_h5) as h: + # convert keyview to list to get polarization key + pol = list(h[raster_key].keys())[0] + + # Extract some info from reference/secondary CSLC products with GDAL + h5_gdal_path = f'HDF5:{path_h5}://{raster_key}/{pol}' + dataset = gdal.Open(h5_gdal_path, gdal.GA_ReadOnly) + geotransform = dataset.GetGeoTransform() + proj = dataset.GetProjection() + slc = dataset.GetRasterBand(1).ReadAsArray() + + return pol, geotransform, proj, slc + + def compare_cslc_products(file_ref, file_sec): ''' Compare a reference and a newly generated @@ -48,38 +85,28 @@ def compare_cslc_products(file_ref, file_sec): return # Extract some info from reference/secondary CSLC products - dataset_ref = gdal.Open(file_ref, gdal.GA_ReadOnly) - geotransform_ref = dataset_ref.GetGeoTransform() - ref_proj = dataset_ref.GetProjection() - nbands_ref = dataset_ref.RasterCount - - dataset_sec = gdal.Open(file_sec, gdal.GA_ReadOnly) - geotransform_sec = dataset_sec.GetGeoTransform() - sec_proj = dataset_sec.GetProjection() - nbands_sec = dataset_sec.RasterCount + pol_ref, geotransform_ref, proj_ref, slc_ref = _gdal_nfo_retrieve(file_ref) + pol_sec, geotransform_sec, proj_sec, slc_sec = _gdal_nfo_retrieve(file_sec) # Compare number of bands print('Comparing CSLC number of bands ...') - if not nbands_ref == nbands_sec: - print(f'ERROR Number of bands in reference CSLC {nbands_ref} differs' - f'from number of bands {nbands_sec} in secondary CSLC') + if not pol_ref == pol_sec: + print(f'ERROR Polarization in reference CSLC {pol_ref} differs ' + f'from polarization {pol_sec} in secondary CSLC') return print('Comparing CSLC projection ...') - if not ref_proj == sec_proj: - print(f'ERROR projection in reference CSLC {ref_proj} differs' - f'from projection in secondary CSLC {sec_proj}') + if not proj_ref == proj_sec: + print(f'ERROR projection in reference CSLC {proj_ref} differs ' + f'from projection in secondary CSLC {proj_sec}') print('Comparing geo transform arrays ...') if not np.array_equal(geotransform_ref, geotransform_sec): - print(f'ERROR Reference geo transform array {dataset_ref} differs' - f'from secondary CSLC geo transform array {dataset_sec}') + print(f'ERROR Reference geo transform array {geotransform_ref} differs' + f'from secondary CSLC geo transform array {geotransform_sec}') return # Compare amplitude of reference and generated CSLC products - slc_ref = dataset_ref.GetRasterBand(1).ReadAsArray() - slc_sec = dataset_sec.GetRasterBand(1).ReadAsArray() - print('Check mean real part difference between CSLC products is < 1.0e-5') assert np.allclose(slc_ref.real, slc_sec.real, atol=0.0, rtol=1.0e-5, equal_nan=True) @@ -87,7 +114,35 @@ def compare_cslc_products(file_ref, file_sec): assert np.allclose(slc_ref.imag, slc_sec.imag, atol=0.0, rtol=1.0e-5, equal_nan=True) - return + +def _get_metadata_keys(path_h5): + """ + Extract metadata group/dataset names of a given HDF5 + + Parameters + ---------- + path_h5: str + File path to CSLC HDF5 product + + Returns + ------- + metadata_dict: str + Dict holding metadata group and dataset names where: + keys: metadata key names representing datasets or groups + values: set of key names belonging to each metadata key names + """ + metadata_dict = {} + with h5py.File(path_h5, 'r') as h: + metadata = h['metadata'] + # get metadata keys and iterate + metadata_keys = set(metadata.keys()) + for metadata_key in metadata_keys: + if not isinstance(metadata[metadata_key], h5py.Group): + continue + # save keys to current metadata key + metadata_dict[metadata_key] = set(metadata[metadata_key].keys()) + + return metadata_dict def compare_cslc_metadata(file_ref, file_sec): @@ -111,12 +166,9 @@ def compare_cslc_metadata(file_ref, file_sec): print(f'ERROR CSLC metadata not found: {file_sec}') return - # Load metadata - with open(file_ref, 'r') as f: - metadata_ref = json.load(f) - - with open(file_sec, 'r') as f: - metadata_sec = json.load(f) + # Get metadata keys + metadata_ref = _get_metadata_keys(file_ref) + metadata_sec = _get_metadata_keys(file_sec) # Intersect metadata keys set_ref_minus_sec = metadata_ref.keys() - metadata_sec.keys() @@ -127,15 +179,24 @@ def compare_cslc_metadata(file_ref, file_sec): err_str += f'\nReference CSLC metadata extra entries: {set_ref_minus_sec}' if set_sec_minus_ref: err_str += f'\nSecondary CSLC metadata extra entries: {set_sec_minus_ref}' + # Check if metadata key differ assert not set_ref_minus_sec or not set_sec_minus_ref, err_str - # Check remaining metadatakeys - for k_ref, v_ref in metadata_ref.items(): - if metadata_sec[k_ref] != v_ref: - print(f'ERROR the content of metadata key {k_ref} from' - f'reference CSLC metadata has a value {v_ref} whereas the same' - f'key in the secondary CSLC metadata has value {metadata_sec[k_ref]}') + # Check sub metadatakeys (after establishing top level metadata matches) + for key in metadata_ref.keys(): + # Intersect metadata keys + set_ref_minus_sec = metadata_ref[key] - metadata_sec[key] + set_sec_minus_ref = metadata_sec[key] - metadata_ref[key] + + err_str = "Metadata keys do not match.\n" + if set_ref_minus_sec: + err_str += f'\nReference CSLC {key} metadata extra entries: {set_ref_minus_sec}' + if set_sec_minus_ref: + err_str += f'\nSecondary CSLC {key} metadata extra entries: {set_sec_minus_ref}' + + # Check if metadata key differ + assert not set_ref_minus_sec or not set_sec_minus_ref, err_str if __name__ == '__main__': @@ -146,5 +207,5 @@ def compare_cslc_metadata(file_ref, file_sec): print('All CSLC product checks have passed') # Check CSLC metadata - compare_cslc_metadata(cmd.ref_metadata, cmd.sec_metadata) + compare_cslc_metadata(cmd.ref_product, cmd.sec_product) print('All CSLC metadata checks have passed')