Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Placeholder for geocoded noise LUT #119

Merged
merged 7 commits into from
Mar 22, 2023
Merged

Placeholder for geocoded noise LUT #119

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
0729b20
noise info moved to `metadata/noise_information`
Mar 21, 2023
2958e6a
`geocode_noise_luts` implemented
Mar 21, 2023
2c492be
removed radargrid-based information in `noise_information`
Mar 21, 2023
50f0826
codacy issue
Mar 21, 2023
0ced2cf
reduce duplicate code
LiangJYu Mar 22, 2023
2137a4c
fix on `lut_path` definition
Mar 22, 2023
ab80ff0
Apply suggestions from code review
seongsujeong Mar 22, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
133 changes: 82 additions & 51 deletions src/compass/s1_geocode_metadata.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down Expand Up @@ -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.
Expand All @@ -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
'''
Expand All @@ -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,
Expand All @@ -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)
Expand All @@ -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,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function has a lot of code repetition with the calibration LUT function. I am wondering whether is possible to consolidate the two in a single function

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree with your concern. The code redundancy has reduced thanks to @LiangJYu.

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()
Expand Down
6 changes: 6 additions & 0 deletions src/compass/s1_geocode_slc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)")

Expand Down
40 changes: 13 additions & 27 deletions src/compass/utils/h5_helpers.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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']
Expand Down