Skip to content

Commit

Permalink
HDF5 output (opera-adt#62)
Browse files Browse the repository at this point in the history 
* 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 <[email protected]>

* equal spacing in y

Co-authored-by: Scott Staniewicz <[email protected]>

* fix paths
correctly place CF attr

* fix indentation

* add type at h5 root
extension back to h5

Co-authored-by: Scott Staniewicz <[email protected]>
Co-authored-by: Scott Staniewicz <[email protected]>
  • Loading branch information
@LiangJYu @scottstanie
3 people authored Dec 5, 2022
1 parent b7de1d7 commit b5e9341
Show file tree
Hide file tree
Showing 14 changed files with 766 additions and 512 deletions.
13 changes: 3 additions & 10 deletions src/compass/s1_geo2rdr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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
Expand All @@ -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)")
Expand Down
99 changes: 48 additions & 51 deletions src/compass/s1_geocode_metadata.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,26 +6,27 @@
import os
import time

import h5py
import isce3
import journal
import numpy as np

from osgeo import gdal
from compass.utils.runconfig import RunConfig
from compass.utils.helpers import get_module_name
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
----------
cfg: dict
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")
Expand All @@ -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()
Expand All @@ -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(
Expand Down
108 changes: 59 additions & 49 deletions src/compass/s1_geocode_slc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
'''
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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',
Expand All @@ -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

Expand All @@ -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)")
Expand Down
Loading

0 comments on commit b5e9341

Please sign in to comment.