diff --git a/src/compass/s1_geocode_metadata.py b/src/compass/s1_geocode_metadata.py index f9dffb72..ef2a0054 100755 --- a/src/compass/s1_geocode_metadata.py +++ b/src/compass/s1_geocode_metadata.py @@ -45,7 +45,6 @@ def run(cfg, burst, fetch_from_scratch=False): threshold = cfg.geo2rdr_params.threshold iters = cfg.geo2rdr_params.numiter lines_per_block = cfg.geo2rdr_params.lines_per_block - output_format = cfg.geocoding_params.output_format # process one burst only date_str = burst.sensing_start.strftime("%Y%m%d") @@ -126,8 +125,8 @@ def run(cfg, burst, fetch_from_scratch=False): f"{module_name} burst successfully ran in {dt} (hr:min:sec)") -def geocode_calibration_luts(geo_burst_h5, burst, cfg, - dec_factor=40): +def geocode_luts(geo_burst_h5, burst, cfg, dst_group_path, item_dict, + dec_factor=40): ''' Geocode the radiometric calibratio paremeters, and write them into output HDF5. @@ -140,6 +139,10 @@ def geocode_calibration_luts(geo_burst_h5, burst, cfg, Sentinel-1 burst SLC cfg: GeoRunConfig GeoRunConfig object with user runconfig options + dst_group_path: str + Path in HDF5 where geocode rasters will be placed + item_dict: dict + Dict containing item names and values to be geocoded dec_factor: int Decimation factor to downsample the slant range pixels for LUT ''' @@ -160,17 +163,9 @@ def geocode_calibration_luts(geo_burst_h5, burst, cfg, iters = cfg.geo2rdr_params.numiter scratch_path = out_paths.scratch_directory - # Designate radiometric calibration parameter to geocode - calibration_dict = { - 'gamma':burst.burst_calibration.gamma, - 'sigma_naught':burst.burst_calibration.sigma_naught, - 'dn':burst.burst_calibration.dn - } - - # define the geogrid for calbration LUT - calibration_radargrid = radar_grid.multilook(dec_factor, - dec_factor) - calibration_geogrid = isce3.product.GeoGridParameters( + # generate decimated radar and geo grids for LUT(s) + decimated_radargrid = radar_grid.multilook(dec_factor, dec_factor) + decimated_geogrid = isce3.product.GeoGridParameters( geo_grid.start_x, geo_grid.start_y, geo_grid.spacing_x * dec_factor, @@ -186,70 +181,59 @@ def geocode_calibration_luts(geo_burst_h5, burst, cfg, geocode_obj.doppler = isce3.core.LUT2d() geocode_obj.threshold_geo2rdr = threshold geocode_obj.numiter_geo2rdr = iters - geocode_obj.geogrid(calibration_geogrid.start_x, - calibration_geogrid.start_y, - calibration_geogrid.spacing_x, - calibration_geogrid.spacing_y, - calibration_geogrid.width, - calibration_geogrid.length, - calibration_geogrid.epsg) - calibration_group_path =\ - f'{ROOT_PATH}/metadata/calibration_information' - calibration_group =\ - geo_burst_h5.require_group(calibration_group_path) + geocode_obj.geogrid(decimated_geogrid.start_x, + decimated_geogrid.start_y, + decimated_geogrid.spacing_x, + decimated_geogrid.spacing_y, + decimated_geogrid.width, + decimated_geogrid.length, + decimated_geogrid.epsg) + dst_group =\ + geo_burst_h5.require_group(dst_group_path) gdal_envi_driver = gdal.GetDriverByName('ENVI') - for calibration_key, _ in calibration_dict.items(): + for item_name, _ in item_dict.items(): # prepare input dataset in output HDF5 - init_geocoded_dataset(calibration_group, - calibration_key, - calibration_geogrid, + init_geocoded_dataset(dst_group, + item_name, + decimated_geogrid, 'float32', - f'geocoded {calibration_key}') + f'geocoded {item_name}') - calibration_dataset =\ - geo_burst_h5[f'{calibration_group_path}/{calibration_key}'] + dst_dataset = geo_burst_h5[f'{dst_group_path}/{item_name}'] # prepare output raster geocoded_cal_lut_raster =\ isce3.io.Raster( - f"IH5:::ID={calibration_dataset.id.id}".encode("utf-8"), - update=True) + f"IH5:::ID={dst_dataset.id.id}".encode("utf-8"), update=True) # populate and prepare radargrid LUT input raster # NOTE: `lut_arr` below is a placeholder, which will be # eventually replaced by LUTs for geocoded calibration parameters. - lut_arr = np.zeros((calibration_radargrid.length, - calibration_radargrid.width)) + lut_arr = np.zeros((decimated_radargrid.length, + decimated_radargrid.width)) + lut_path = f'{scratch_path}/{item_name}_radargrid.rdr' lut_gdal_raster = gdal_envi_driver.Create( - f'{scratch_path}/{calibration_key}_radargrid.rdr', - calibration_radargrid.width, - calibration_radargrid.length, - 1, - gdal.GDT_Float32) + lut_path, decimated_radargrid.width, + decimated_radargrid.length, 1, gdal.GDT_Float32) lut_band = lut_gdal_raster.GetRasterBand(1) lut_band.WriteArray(lut_arr) lut_band.FlushCache() lut_gdal_raster = None - input_raster =\ - isce3.io.Raster(f'{scratch_path}/' - f'{calibration_key}_radargrid.rdr') + input_raster = isce3.io.Raster(lut_path) # geocode then set transfrom and EPSG in output raster - geocode_obj.geocode(radar_grid=calibration_radargrid, + geocode_obj.geocode(radar_grid=decimated_radargrid, input_raster=input_raster, output_raster=geocoded_cal_lut_raster, dem_raster=dem_raster, output_mode=isce3.geocode.GeocodeOutputMode.INTERP) - geotransform=[calibration_geogrid.start_x, - calibration_geogrid.spacing_x, - 0, - calibration_geogrid.start_y, - 0, - calibration_geogrid.spacing_y] + geotransform = [decimated_geogrid.start_x, decimated_geogrid.spacing_x, + 0, decimated_geogrid.start_y, 0, + decimated_geogrid.spacing_y] geocoded_cal_lut_raster.set_geotransform(geotransform) geocoded_cal_lut_raster.set_epsg(epsg) @@ -258,6 +242,53 @@ def geocode_calibration_luts(geo_burst_h5, burst, cfg, del geocoded_cal_lut_raster +def geocode_calibration_luts(geo_burst_h5, burst, cfg, + dec_factor=40): + ''' + Geocode the radiometric calibratio paremeters, + and write them into output HDF5. + + Parameters + ---------- + geo_burst_h5: h5py.files.File + HDF5 object as the output product + burst: s1reader.Sentinel1BurstSlc + Sentinel-1 burst SLC + cfg: GeoRunConfig + GeoRunConfig object with user runconfig options + dec_factor: int + Decimation factor to downsample the slant range pixels for LUT + ''' + dst_group_path = f'{ROOT_PATH}/metadata/calibration_information' + item_dict = {'gamma':burst.burst_calibration.gamma, + 'sigma_naught':burst.burst_calibration.sigma_naught, + 'dn':burst.burst_calibration.dn} + geocode_luts(geo_burst_h5, burst, cfg, dst_group_path, item_dict, + dec_factor) + + +def geocode_noise_luts(geo_burst_h5, burst, cfg, + dec_factor=40): + ''' + Geocode the noise LUT, and write that into output HDF5. + + Parameters + ---------- + geo_burst_h5: h5py.files.File + HDF5 object as the output product + burst: s1reader.Sentinel1BurstSlc + Sentinel-1 burst SLC + cfg: GeoRunConfig + GeoRunConfig object with user runconfig options + dec_factor: int + Decimation factor to downsample the slant range pixels for LUT + ''' + dst_group_path = f'{ROOT_PATH}/metadata/noise_information' + item_dict = {'thermal_noise_lut': None} + geocode_luts(geo_burst_h5, burst, cfg, dst_group_path, item_dict, + dec_factor) + + if __name__ == "__main__": ''' run geocode metadata layers from command line''' parser = YamlArgparse() diff --git a/src/compass/s1_geocode_slc.py b/src/compass/s1_geocode_slc.py index 96ab941d..6474f2bb 100755 --- a/src/compass/s1_geocode_slc.py +++ b/src/compass/s1_geocode_slc.py @@ -224,6 +224,12 @@ def run(cfg: GeoRunConfig): burst, cfg) + if burst.burst_noise is not None: + # Geocode the calibration parameters and write them into HDF5 + s1_geocode_metadata.geocode_noise_luts(geo_burst_h5, + burst, + cfg) + 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/utils/h5_helpers.py b/src/compass/utils/h5_helpers.py index 5dc03875..ab15cf30 100644 --- a/src/compass/utils/h5_helpers.py +++ b/src/compass/utils/h5_helpers.py @@ -422,6 +422,19 @@ def metadata_to_h5group(parent_group, burst, cfg): for meta_item in cal_items: add_dataset_and_attrs(cal_group, meta_item) + # write out noise metadata, if present + if burst.burst_noise is not None: + noise = burst.burst_noise + noise_items = [ + Meta('basename', noise.basename_nads, ''), + Meta('range_azimuth_time', + noise.range_azimuth_time.strftime(TIME_STR_FMT), + 'Start time', {'format': 'YYYY-MM-DD HH:MM:SS.6f'}) + ] + noise_group = meta_group.require_group('noise_information') + for meta_item in noise_items: + add_dataset_and_attrs(noise_group, meta_item) + # runconfig yaml text processing_group['runconfig'] = cfg.yaml_string @@ -699,33 +712,6 @@ def corrections_to_h5group(parent_group, burst, cfg, rg_lut, az_lut, for meta_item in eap_items: add_dataset_and_attrs(eap_group, meta_item) - # write out noise metadata, if present - if burst.burst_noise is not None: - noise = burst.burst_noise - noise_items = [ - Meta('basename', noise.basename_nads, ''), - Meta('range_azimuth_time', - noise.range_azimuth_time.strftime(TIME_STR_FMT), - 'Start time', {'format': 'YYYY-MM-DD HH:MM:SS.6f'}), - Meta('range_line', noise.range_line, 'Range line'), - Meta('range_pixel', noise.range_pixel, 'Range array in pixel for LUT'), - Meta('range_lut', noise.range_lut, 'Range noise lookup table data'), - Meta('azimuth_first_azimuth_line', noise.azimuth_first_azimuth_line, - 'First line of the burst in subswath. NaN if not available in annotation.'), - Meta('azimuth_first_range_sample', noise.azimuth_first_range_sample, - 'First range sample of the burst. NaN if not available in annotation.'), - Meta('azimuth_last_azimuth_line', noise.azimuth_last_azimuth_line, - 'Last line of the burst in subswatn. NaN if not available in annotation.'), - Meta('azimuth_last_range_sample', noise.azimuth_last_range_sample, - 'Last range of the burst. NaN if not available in annotation.'), - Meta('azimuth_line', noise.azimuth_line, 'azimuth line index for noise LUT'), - Meta('azimuth_lut', noise.azimuth_lut, 'azimuth noise lookup table data') - ] - noise_group = correction_group.require_group('noise') - for meta_item in noise_items: - add_dataset_and_attrs(noise_group, meta_item) - - def get_cslc_geotransform(filename, pol: str = "VV"): gdal_str = f'NETCDF:{filename}:/{GRID_PATH}/{pol}' return gdal.Info(gdal_str, format='json')['geoTransform']