diff --git a/cxx/isce3/geometry/loadDem.h b/cxx/isce3/geometry/loadDem.h index d4c7ae57d..eb35139c8 100644 --- a/cxx/isce3/geometry/loadDem.h +++ b/cxx/isce3/geometry/loadDem.h @@ -61,7 +61,7 @@ isce3::geometry::DEMInterpolator DEMRasterToInterpolator( * (when the DEM is in geographic coordinates and `proj` is in polar stereo) * @param[in] minY Minimum Y/northing position * @param[in] maxY Maximum Y/northing position -* @param[out] dem_interp_block DEM interpolation object +* @param[out] dem_interp DEM interpolation object * @param[in] proj Projection object (nullptr to use same * DEM projection) * @param[in] dem_margin_x_in_pixels DEM X/easting margin in pixels diff --git a/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp b/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp index a82cab5d4..a9cf51f98 100644 --- a/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp +++ b/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp @@ -218,3 +218,84 @@ void addbinding_DEM_raster2interpolator(py::module& m) )") ; } + + +void addbinding_load_dem_from_proj(py::module& m) +{ + m.def("load_dem_from_proj", + [](isce3::io::Raster &dem_raster, + const double x0, + const double xf, + const double minY, + const double maxY, + const isce3::core::dataInterpMethod dem_interp_method, + isce3::core::ProjectionBase* proj, + const int dem_margin_x_in_pixels, + const int dem_margin_y_in_pixels, + const int dem_raster_band) { + + DEMInterp dem_interp(0, dem_interp_method); + + isce3::geometry::loadDemFromProj(dem_raster, + x0, + xf, + minY, + maxY, + &dem_interp, + proj, + dem_margin_x_in_pixels, + dem_margin_y_in_pixels, + dem_raster_band); + + return dem_interp; + }, + py::arg("dem_raster"), + py::arg("x0"), + py::arg("xf"), + py::arg("min_y"), + py::arg("max_y"), + py::arg("dem_interp_method") = isce3::core::BIQUINTIC_METHOD, + py::arg("proj") = nullptr, + py::arg("dem_margin_x_in_pixels") = 100, + py::arg("dem_margin_y_in_pixels") = 200, + py::arg("dem_raster_band") = 1, + R"( + Load DEM raster into a DEMInterpolator object around a given bounding box + in the same or different coordinate system as the DEM raster + Parameters + ---------- + dem_raster: isce3.io.Raster + Raster of the DEM + x0: double + If the DEM is in geographic coordinates and the `x0` coordinate is not + from the polar stereo system EPSG 3031 or EPSG 3413, this point represents + the minimum X coordinate value. In this case, the maximum + longitude span that this function can handle is 180 degrees + (when the DEM is in geographic coordinates and `proj` is in polar stereo + xf: double + Easting/longitude of eastern edge of bounding box + If the DEM is in geographic coordinates and the `xf` coordinate is not + from the polar stereo system EPSG 3031 or EPSG 3413, this point represents + the maximum X coordinate value. In this case, the maximum + longitude span that this function can handle is 180 degrees + (when the DEM is in geographic coordinates and `proj` is in polar stereo) + min_y: double + Minimum Y/northing position + max_y: double + Maximum Y/northing position + dem_interp_method: isce3.core.DataInterpMethod + DEM interpolation method + proj: + Projection object (nullptr to use same DEM projection) + dem_margin_x_in_pixels, int + DEM X/easting margin in pixels + dem_margin_y_in_pixels, int + DEM Y/northing margin in pixels + dem_raster_band: int + DEM raster band (starting from 1) + Returns + ------- + dem_interp: isce3.geometry.DEMInterpolator + DEM interpolator for given DEM raster and geo grid. + )"); +} diff --git a/python/extensions/pybind_isce3/geometry/DEMInterpolator.h b/python/extensions/pybind_isce3/geometry/DEMInterpolator.h index f5ccd765a..1fb9efca9 100644 --- a/python/extensions/pybind_isce3/geometry/DEMInterpolator.h +++ b/python/extensions/pybind_isce3/geometry/DEMInterpolator.h @@ -5,3 +5,4 @@ void addbinding(pybind11::class_&); void addbinding_DEM_raster2interpolator(pybind11::module&); +void addbinding_load_dem_from_proj(pybind11::module&); diff --git a/python/extensions/pybind_isce3/geometry/geometry.cpp b/python/extensions/pybind_isce3/geometry/geometry.cpp index 8193aa3bd..822985ed4 100644 --- a/python/extensions/pybind_isce3/geometry/geometry.cpp +++ b/python/extensions/pybind_isce3/geometry/geometry.cpp @@ -72,4 +72,5 @@ void addsubmodule_geometry(py::module & m) addbinding_pnt_intersect(geometry); addbinding_look_inc_from_sr(geometry); addbinding_DEM_raster2interpolator(geometry); + addbinding_load_dem_from_proj(geometry); } diff --git a/python/packages/isce3/geometry/bounding_radar_grid.py b/python/packages/isce3/geometry/bounding_radar_grid.py index d583da79f..2c06b6247 100644 --- a/python/packages/isce3/geometry/bounding_radar_grid.py +++ b/python/packages/isce3/geometry/bounding_radar_grid.py @@ -303,6 +303,11 @@ def get_bounding_radar_grid( Azimuth time spacing of the output grid, in seconds. Must be > 0. rg_spacing : float Slant range spacing of the output grid, in meters. Must be > 0. + sensing_start : float + The sensing start time of the input `geo_grid`, in seconds since the epoch of + `orbit`. + starting_range : float + The starting slant range of the input `geo_grid`, in meters. orbit : isce3.core.Orbit The trajectory of the radar antenna phase center over a time interval that includes the observation times of each point in `geo_grid` at each height diff --git a/python/packages/nisar/static/ephemeris.py b/python/packages/nisar/static/ephemeris.py index 5fc249d83..ad6eb445e 100644 --- a/python/packages/nisar/static/ephemeris.py +++ b/python/packages/nisar/static/ephemeris.py @@ -4,7 +4,8 @@ from datetime import datetime from nisar.products.readers.attitude import load_attitude_from_xml -from nisar.products.readers.orbit import load_orbit_from_xml +from nisar.products.readers.orbit import load_orbit_from_xml, load_orbit +from nisar.products.readers import SLC import isce3 @@ -13,10 +14,11 @@ def get_cropped_orbit_and_attitude( - orbit_xml_file: str | os.PathLike, - pointing_xml_file: str | os.PathLike, - start_time: str | datetime | None, - end_time: str | datetime | None, + input_file_path: str | os.PathLike | None = None, + orbit_xml_file: str | os.PathLike | None = None, + pointing_xml_file: str | os.PathLike | None = None, + start_time: str | datetime | None = None, + end_time: str | datetime | None = None, *, padding: float = 0.0, ) -> tuple[isce3.core.Orbit, isce3.core.Attitude]: @@ -31,10 +33,12 @@ def get_cropped_orbit_and_attitude( Parameters ---------- - orbit_xml_file : path-like + input_file_path : str or os.PathLike or None + Path to the input NISAR L1 RSLC formatted HDF5 file. + orbit_xml_file : path-like or None Path to the input orbit ephemeris XML file. Must be an existing XML file conforming to the NISAR Orbit Ephemeris Product Specification\ [1]_. - pointing_xml_file : path-like + pointing_xml_file : path-like or None Path to the input radar pointing XML file. Must be an existing XML file conforming to the NISAR Radar Pointing Product Specification\ [2]_. start_time : str or datetime.datetime or None @@ -70,16 +74,48 @@ def get_cropped_orbit_and_attitude( logger = get_logger() - # Load ephemeris data from input XML files. - logger.info(f"Load orbit data from file {orbit_xml_file}") - orbit_full = load_orbit_from_xml(orbit_xml_file) + if orbit_xml_file is not None: + # Load ephemeris data from input XML files. + logger.info(f"Load orbit data from file {orbit_xml_file}") + + if input_file_path is not None: + # Ensure the orbit is referenced to the RSLC radar grid + # reference epoch. + rslc_product = SLC(hdf5file=str(input_file_path)) + rslc_radar_grid = rslc_product.getRadarGrid() + orbit_full = load_orbit(rslc_product, orbit_xml_file, + rslc_radar_grid.ref_epoch) + else: + orbit_full = load_orbit_from_xml(orbit_xml_file) + + elif input_file_path is not None: + # Load ephemeris data from input RSLC HDF5 file. + logger.info(f"Load orbit data from RSLC file {input_file_path}") + rslc_product = SLC(hdf5file=str(input_file_path)) + orbit_full = rslc_product.getOrbit() + else: + raise ValueError( + "Either the RSLC HDF5 or the orbit XML file must be provided" + ) + logger.info( "Original orbit data spans time interval" f" [{orbit_full.start_datetime, orbit_full.end_datetime}]" ) - logger.info(f"Load attitude data from file {pointing_xml_file}") - attitude_full = load_attitude_from_xml(pointing_xml_file) + if pointing_xml_file is not None: + logger.info(f"Load attitude data from file {pointing_xml_file}") + attitude_full = load_attitude_from_xml(pointing_xml_file) + elif input_file_path is not None: + # Load attitude data from input RSLC HDF5 file. + logger.info(f"Load attitude data from RSLC file {input_file_path}") + rslc_product = SLC(hdf5file=str(input_file_path)) + attitude_full = rslc_product.getAttitude() + else: + raise ValueError( + "Either the RSLC HDF5 or the pointing XML file must be provided" + ) + logger.info( "Original attitude data spans time interval" f" [{attitude_full.start_datetime, attitude_full.end_datetime}]" diff --git a/python/packages/nisar/static/product.py b/python/packages/nisar/static/product.py index 4a8176a7c..7352a22a2 100644 --- a/python/packages/nisar/static/product.py +++ b/python/packages/nisar/static/product.py @@ -322,7 +322,7 @@ def populate_grids_group( local_incidence_angle: isce3.io.Raster, line_of_sight_x: isce3.io.Raster, line_of_sight_y: isce3.io.Raster, - water_mask: isce3.io.Raster, + water_mask: isce3.io.Raster | None, rtc_gamma_to_sigma_factor: isce3.io.Raster, rtc_gamma_to_beta_factor: isce3.io.Raster, geo_grid: isce3.product.GeoGridParameters, @@ -439,8 +439,11 @@ def create_raster_layer_dataset(name: str, raster: isce3.io.Raster) -> h5py.Data dem_dataset = create_raster_layer_dataset("digitalElevationModel", reprojected_dem) dem_dataset.attrs["disclaimer"] = to_bytes(dem_disclaimer) - water_mask_dataset = create_raster_layer_dataset("waterMask", water_mask) - water_mask_dataset.attrs["disclaimer"] = to_bytes(water_mask_disclaimer) + if water_mask is not None: + water_mask_dataset = create_raster_layer_dataset("waterMask", + water_mask) + water_mask_dataset.attrs["disclaimer"] = to_bytes( + water_mask_disclaimer) create_raster_layer_dataset("layoverShadowMask", layover_shadow_mask) create_raster_layer_dataset("localIncidenceAngle", local_incidence_angle) diff --git a/python/packages/nisar/workflows/static.py b/python/packages/nisar/workflows/static.py index e7165519b..c340ba553 100644 --- a/python/packages/nisar/workflows/static.py +++ b/python/packages/nisar/workflows/static.py @@ -7,6 +7,7 @@ from collections.abc import Sequence from datetime import datetime, timezone from pathlib import Path +from xmlrpc.client import DateTime import h5py import nisar @@ -15,7 +16,8 @@ from nisar.static.geo_grid import get_output_geo_grid from nisar.static.geometry_layers import compute_geometry_layers from nisar.static.granule_id import form_granule_id -from nisar.static.layover_shadow_mask import compute_geocoded_layover_shadow_mask +from nisar.static.layover_shadow_mask import \ + compute_geocoded_layover_shadow_mask from nisar.static.logging import get_logger, log_elapsed_time from nisar.static.product import ( build_hdf5_dataset_creation_kwds_dict, @@ -25,16 +27,21 @@ ) from nisar.static.rtc_anf_layers import compute_rtc_anf_layers from nisar.static.runconfig import get_runconfig_params -from nisar.static.util import get_raster_dataset_metadata_item, scratch_directory +from nisar.static.util import get_raster_dataset_metadata_item, \ + scratch_directory from nisar.static.water_mask import binarize_and_reproject_water_mask import isce3 -from isce3.geometry import make_geo_grid_bounding_polygon +from isce3.geometry import make_geo_grid_bounding_polygon, load_dem_from_proj +from isce3.core import normalize_look_side, normalize_data_interp_method +from nisar.products.readers import SLC +import numpy as np def run_static_layers_workflow(config_file: os.PathLike | str) -> None: """ - Run the NISAR Static Layers workflow with the specified run configuration file. + Run the NISAR Static Layers workflow with the specified run configuration + file. Will generate a single STATIC HDF5 granule, as specified in the runconfig. @@ -58,6 +65,7 @@ def run_static_layers_workflow(config_file: os.PathLike | str) -> None: output_params = groups["output"] # Open the input DEM and water mask raster datasets. + input_file_path = dynamic_ancillary_files["input_file_path"] dem_raster_file = dynamic_ancillary_files["dem_raster_file"] water_mask_raster_file = dynamic_ancillary_files["water_mask_raster_file"] logger.info(f"Open DEM raster file {dem_raster_file}") @@ -65,62 +73,221 @@ def run_static_layers_workflow(config_file: os.PathLike | str) -> None: # Construct a DEM interpolator. dem_interp_method = processing_params["dem"]["interp_method"] - dem = isce3.geometry.DEMInterpolator(dem_raster) - dem.interp_method = dem_interp_method # Construct the output geocoded coordinate grid. geo_grid_params = processing_params["geo_grid"] geo_grid = get_output_geo_grid(dem_raster=dem_raster, **geo_grid_params) logger.info(f"Output geo grid: {geo_grid}") - # Parse the orbit and attitude data from the input XML files. Crop the data to the - # time interval of interest to avoid possible geo2rdr convergence errors due to - # ambiguity between orbit periods. + flag_save_water_mask = output_params["layers"]["save_water_mask"] + + proj = isce3.core.make_projection(geo_grid.epsg) + + # dem = isce3.geometry.DEMInterpolator(dem_raster) + # dem.interp_method = dem_interp_method + dem_interp = load_dem_from_proj( + dem_raster, + geo_grid.start_x, + geo_grid.end_x, + geo_grid.end_y, + geo_grid.start_y, + normalize_data_interp_method(dem_interp_method), + proj) + + # Load the orbit and attitude data from the input RSLC or XML files. + # Crop the data to the time interval of interest to avoid possible + # geo2rdr convergence errors due to ambiguity between orbit periods. orbit, attitude = get_cropped_orbit_and_attitude( + input_file_path=input_file_path, orbit_xml_file=dynamic_ancillary_files["orbit_xml_file"], pointing_xml_file=dynamic_ancillary_files["pointing_xml_file"], **processing_params["ephemeris"], ) - # Get the Doppler centroid associated with the radar grid. NISAR image grids are - # always zero-Doppler. + # Get the Doppler centroid associated with the radar grid. NISAR image + # grids are always zero-Doppler. img_grid_doppler = isce3.core.LUT2d() - # Estimate the required radar grid spacing necessary to avoid undersampling the - # output geocoded grid. + # Estimate the required radar grid spacing necessary to avoid + # undersampling the output geocoded grid. # XXX: We deliberately don't pass geo2rdr parameters to either - # `infer_radar_grid_spacing_from_geo_grid()` or `get_bounding_radar_grid()` because - # these functions use `geo2rdr_bracket`, which takes different parameters than the - # legacy `geo2rdr` routine that's used by most of the workflow. Exposing both sets - # of parameters would introduce a lot of additional bookkeeping for seemingly little - # benefit. - logger.info("Estimate maximum required radar grid spacing") + # `infer_radar_grid_spacing_from_geo_grid()` or `get_bounding_radar_grid()` + # because these functions use `geo2rdr_bracket`, which takes different + # parameters than the legacy `geo2rdr` routine that's used by most of the + # workflow. Exposing both sets of parameters would introduce a lot of + # additional bookkeeping for seemingly little benefit. radar_grid_params = processing_params["radar_grid"] look_side = radar_grid_params["look_side"] wavelength = radar_grid_params["wavelength"] - az_spacing, rg_spacing = isce3.geometry.infer_radar_grid_spacing_from_geo_grid( - geo_grid=geo_grid, - dem=dem, - orbit=orbit, - doppler=img_grid_doppler, - look_side=look_side, - wavelength=wavelength, - **radar_grid_params["spacing"], - ) - # Compute a radar grid whose footprint on the ground encloses the geocoded grid on - # which each output layer is defined. - logger.info("Compute a radar grid spanning the region of interest") - radar_grid = isce3.geometry.get_bounding_radar_grid( - geo_grid=geo_grid, - az_spacing=az_spacing, - rg_spacing=rg_spacing, - orbit=orbit, - look_side=look_side, - wavelength=wavelength, - doppler=img_grid_doppler, - **radar_grid_params["bounding_box"], - ) + radar_grid_spacing_params = radar_grid_params["spacing"] + az_spacing = radar_grid_spacing_params["az_spacing"] + rg_spacing = radar_grid_spacing_params["rg_spacing"] + + pts_per_side = radar_grid_spacing_params["pts_per_side"] + + bounding_box_params = radar_grid_params["bounding_box"] + start_datetime_str = bounding_box_params["start_time"] + end_datetime_str = bounding_box_params["end_time"] + start_range = bounding_box_params["start_range"] + end_range = bounding_box_params["end_range"] + min_height = bounding_box_params["min_height"] + max_height = bounding_box_params["max_height"] + pts_per_edge = bounding_box_params["pts_per_edge"] + az_margin = bounding_box_params["az_margin"] + rg_margin = bounding_box_params["rg_margin"] + + start_time = None + end_time = None + + # Print user radar grid bounding box parameters, if provided. + if (start_datetime_str is not None or start_range is not None or + end_datetime_str is not None or end_range is not None): + + logger.info("Using user-specified radar grid bounding box parameters") + if start_datetime_str is not None: + logger.info(f' start time: {start_datetime_str}') + start_time = isce3.core.DateTime(start_datetime_str) + + if end_datetime_str is not None: + logger.info(f' end time: {end_datetime_str}') + end_time = isce3.core.DateTime(end_datetime_str) + + if start_range is not None: + logger.info(f' start range: {start_range}') + + if end_range is not None: + logger.info(f' end range: {end_range}') + + if rg_spacing is not None: + logger.info(f' range spacing: {rg_spacing}') + if az_spacing is not None: + logger.info(f' azimuth time interval: {az_spacing}') + + # Load radar grid parameters from RSLC (if provided) + if (input_file_path is not None and + (start_time is None or end_time is None or + start_range is None or end_range is None or + rg_spacing is None or az_spacing is None)): + logger.info("Load radar grid parameters from input RSLC file:") + rslc_product = SLC(hdf5file=str(input_file_path)) + rslc_radar_grid = rslc_product.getRadarGrid() + rslc_orbit = rslc_product.getOrbit() + + if start_time is None: + start_time = (rslc_orbit.reference_epoch + + isce3.core.TimeDelta( + rslc_radar_grid.sensing_start)) + logger.info(f" start time: {start_time.isoformat()}") + + if end_time is None: + end_time = (rslc_orbit.reference_epoch + + isce3.core.TimeDelta( + rslc_radar_grid.sensing_stop)) + logger.info(f" end time: {end_time.isoformat()}") + + if start_range is None: + start_range = rslc_radar_grid.starting_range + logger.info(f" start range: {start_range}") + if end_range is None: + end_range = rslc_radar_grid.end_range + logger.info(f" end range: {end_range}") + + if rg_spacing is None: + rg_spacing = rslc_radar_grid.range_pixel_spacing + logger.info(f" range spacing: {rg_spacing}") + + if az_spacing is None: + az_spacing = 1.0 / rslc_radar_grid.prf + logger.info(f" azimuth time interval: {az_spacing}") + + if start_time is not None and az_margin != 0.0: + start_time -= isce3.core.TimeDelta(az_margin) + logger.info(f' start time (adjusted for az. margin) {az_margin}:' + f' {start_time}') + if end_time is not None and az_margin != 0.0: + end_time += isce3.core.TimeDelta(az_margin) + logger.info(f' end time (adjusted for az. margin) {az_margin}:' + f' {end_time}') + if start_range is not None and rg_margin != 0.0: + start_range -= rg_margin + logger.info(f' start range (adjusted for rg margin) {rg_margin}:' + f' {start_range}') + if end_range is not None and rg_margin != 0.0: + end_range += rg_margin + logger.info(f' end range (adjusted for rg margin) {rg_margin}:' + f' {end_range}') + + if rg_spacing is None or az_spacing is None: + az_spacing_inferred, rg_spacing_inferred = \ + isce3.geometry.infer_radar_grid_spacing_from_geo_grid( + geo_grid=geo_grid, + dem=dem_interp, + orbit=orbit, + doppler=img_grid_doppler, + look_side=look_side, + wavelength=wavelength, + pts_per_side=pts_per_side + ) + if rg_spacing is None: + rg_spacing = rg_spacing_inferred + logger.info(f' inferred range spacing: {rg_spacing}') + if az_spacing is None: + az_spacing = az_spacing_inferred + logger.info(f' inferred azimuth time interval: {az_spacing}') + + if (start_time is not None and start_range is not None and + end_time is not None and end_range is not None): + + epoch = orbit.reference_epoch + + t0 = (isce3.core.DateTime(start_time) - epoch).total_seconds() + tf = (isce3.core.DateTime(end_time) - epoch).total_seconds() + + num_az = round((tf - t0) / az_spacing) + num_rg = round((end_range - start_range) / rg_spacing) + + logger.info(f' number of lines: {num_az}') + logger.info(f' number of range samples: {num_rg}') + + radar_grid = isce3.product.RadarGridParameters( + sensing_start=t0, + wavelength=wavelength, + prf=1.0 / az_spacing, + starting_range=start_range, + range_pixel_spacing=rg_spacing, + lookside=normalize_look_side(look_side), + length=num_az, + width=num_rg, + ref_epoch=epoch, + ) + + elif (start_time is not None or start_range is not None or + end_time is not None or end_range is not None): + raise ValueError( + "If specifying radar grid bounding box parameters, must provide" + " all of 'start_time', 'start_range', 'end_time', and" + " 'end_range'" + ) + else: + + # Compute a radar grid whose footprint on the ground encloses the geocoded + # grid on which each output layer is defined. + logger.info("Compute a radar grid spanning the region of interest") + radar_grid = isce3.geometry.get_bounding_radar_grid( + geo_grid=geo_grid, + az_spacing=az_spacing, + rg_spacing=rg_spacing, + orbit=orbit, + look_side=look_side, + wavelength=wavelength, + doppler=img_grid_doppler, + min_height=min_height, + max_height=max_height, + pts_per_edge=pts_per_edge, + az_margin=az_margin, + rg_margin=rg_margin + ) logger.info(f"Using radar grid: {radar_grid}") # Get the native Doppler LUT. @@ -129,18 +296,19 @@ def run_static_layers_workflow(config_file: os.PathLike | str) -> None: radar_grid=radar_grid, orbit=orbit, attitude=attitude, - dem=dem, + dem=dem_interp, **processing_params["doppler"], ) - # Create a (possibly temporary) scratch directory to store intermediate files. + # Create a (possibly temporary) scratch directory to store intermediate + # files. logger.info("Create scratch directory") with scratch_directory( product_paths["scratch_dir"], delete=product_paths["delete_scratch_dir"] ) as scratch_dir: - # Compute static geometry layers (height above ellipsoid, line-of-sight X and Y, - # local incidence angle). Results are stored as GeoTIFF files in the scratch - # directory. + # Compute static geometry layers (height above ellipsoid, + # line-of-sight X and Y, local incidence angle). Results are stored as + # GeoTIFF files in the scratch directory. logger.info("Compute static geometry layers") geo2rdr_params = processing_params["geo2rdr"] with log_elapsed_time(logger.info, "Computing static geometry layers"): @@ -155,13 +323,16 @@ def run_static_layers_workflow(config_file: os.PathLike | str) -> None: dem_interp_method=dem_interp_method, geo2rdr_params=geo2rdr_params, ) - reprojected_dem, los_east, los_north, local_inc_angle = geometry_layers + reprojected_dem, los_east, los_north, local_inc_angle = \ + geometry_layers - # Compute static mask layers (geocoded layover/shadow mask and water mask). + # Compute static mask layers (geocoded layover/shadow mask and water + # mask). # Results are stored as GeoTIFF files in the scratch directory. logger.info("Compute geocoded layover/shadow mask layer") geocode_params = processing_params["geocode"] - with log_elapsed_time(logger.info, "Computing geocoded layover/shadow mask"): + with log_elapsed_time(logger.info, + "Computing geocoded layover/shadow mask"): layover_shadow_mask = compute_geocoded_layover_shadow_mask( radar_grid=radar_grid, orbit=orbit, @@ -178,52 +349,61 @@ def run_static_layers_workflow(config_file: os.PathLike | str) -> None: max_block_size=geocode_params["max_block_size"], ) - logger.info("Compute re-projected binary water mask layer") - with log_elapsed_time(logger.info, "Computing re-projected binary water mask"): - binary_water_mask = binarize_and_reproject_water_mask( - water_distance_raster_file=water_mask_raster_file, - geo_grid=geo_grid, - scratch_dir=scratch_dir, - **processing_params["water_mask"], - ) + if flag_save_water_mask: + logger.info("Compute re-projected binary water mask layer") + with log_elapsed_time(logger.info, + "Computing re-projected binary water mask"): + binary_water_mask = binarize_and_reproject_water_mask( + water_distance_raster_file=water_mask_raster_file, + geo_grid=geo_grid, + scratch_dir=scratch_dir, + **processing_params["water_mask"], + ) + else: + binary_water_mask = None - # Compute radiometric terrain correction (RTC) area normalization factor (ANF) - # layers. Results are stored as GeoTIFF files in the scratch directory. + # Compute radiometric terrain correction (RTC) area normalization + # factor (ANF) layers. Results are stored as GeoTIFF files in the + # scratch directory. logger.info("Compute RTC area normalization factor layers") rtc_params = processing_params["rtc"] - with log_elapsed_time(logger.info, "Computing RTC area normalization layers"): - gamma0_to_beta0_factor, gamma0_to_sigma0_factor = compute_rtc_anf_layers( - radar_grid=radar_grid, - orbit=orbit, - native_doppler=native_doppler, - img_grid_doppler=img_grid_doppler, - geo_grid=geo_grid, - dem_raster=dem_raster, - scratch_dir=scratch_dir, - dem_interp_method=dem_interp_method, - geo2rdr_params=geo2rdr_params, - **geocode_params, - **rtc_params, - ) + with log_elapsed_time(logger.info, + "Computing RTC area normalization layers"): + gamma0_to_beta0_factor, gamma0_to_sigma0_factor = \ + compute_rtc_anf_layers( + radar_grid=radar_grid, + orbit=orbit, + native_doppler=native_doppler, + img_grid_doppler=img_grid_doppler, + geo_grid=geo_grid, + dem_raster=dem_raster, + scratch_dir=scratch_dir, + dem_interp_method=dem_interp_method, + geo2rdr_params=geo2rdr_params, + **geocode_params, + **rtc_params, + ) # Infer the orbit pass direction from the orbit velocity vectors. orbit_pass_direction = isce3.core.get_orbit_pass_direction(orbit) - # Pop 'product_counter' from the dict. This parameter is used to form the - # granule ID but doesn't correspond to any dataset in the 'identification' group - # of the product. The other dict contents will be passed as keyword arguments to - # `populate_identification_group()` below. + # Pop 'product_counter' from the dict. This parameter is used to form + # the granule ID but doesn't correspond to any dataset in the + # 'identification' group of the product. The other dict contents will + # be passed as keyword arguments to `populate_identification_group()` + # below. product_counter = primary_executable_params.pop("product_counter") # Get `validity_start_datetime` from the input parameters as a - # `datetime.datetime` object. If it was passed as a non-quoted string in ISO - # 8601 format, `ruamel.yaml` will have already converted it. Otherwise, manually - # convert it here. + # `datetime.datetime` object. If it was passed as a non-quoted string + # in ISO 8601 format, `ruamel.yaml` will have already converted it. + # Otherwise, manually convert it here. validity_start_datetime = primary_executable_params.pop( "validity_start_datetime" ) if not isinstance(validity_start_datetime, datetime): - validity_start_datetime = datetime.fromisoformat(validity_start_datetime) + validity_start_datetime = \ + datetime.fromisoformat(validity_start_datetime) # Get the unique ID of the granule based on the input parameters. radar_band = primary_executable_params["radar_band"] @@ -237,7 +417,8 @@ def run_static_layers_workflow(config_file: os.PathLike | str) -> None: x_posting=abs(geo_grid.spacing_x), y_posting=abs(geo_grid.spacing_y), validity_start_datetime=validity_start_datetime, - composite_release_id=primary_executable_params["composite_release_id"], + composite_release_id=primary_executable_params[ + "composite_release_id"], processing_center=primary_executable_params["processing_center"], product_counter=product_counter, **geometry_params, @@ -262,28 +443,34 @@ def run_static_layers_workflow(config_file: os.PathLike | str) -> None: # Populate global attributes in the root group of the file. logger.info("Populate global HDF5 attributes") product_spec = nisar.products.get_product_spec("STATIC") - nisar.products.populate_global_attrs_from_spec(hdf5_file, product_spec) + nisar.products.populate_global_attrs_from_spec(hdf5_file, + product_spec) # Get the current processing datetime (truncated to integer seconds # precision). - processing_datetime = datetime.now(timezone.utc).replace(microsecond=0) + processing_datetime = \ + datetime.now(timezone.utc).replace(microsecond=0) - # XXX: It's not really obvious what should go in the `zeroDopplerStartTime` - # and `zeroDopplerEndTime` datasets in the 'identification' group. For now, - # we'll use the start & stop time of the radar grid, which is roughly + # XXX: It's not really obvious what should go in the + # `zeroDopplerStartTime` and `zeroDopplerEndTime` datasets in the + # 'identification' group. For now, we'll use the start & stop time + # of the radar grid, which is roughly # analogous what they represent in other NISAR L2 products. - img_grid_start_datetime = radar_grid.ref_epoch + isce3.core.TimeDelta( - radar_grid.sensing_start - ) - img_grid_end_datetime = radar_grid.ref_epoch + isce3.core.TimeDelta( - radar_grid.sensing_stop - ) + img_grid_start_datetime = (radar_grid.ref_epoch + + isce3.core.TimeDelta( + radar_grid.sensing_start)) + img_grid_end_datetime = (radar_grid.ref_epoch + + isce3.core.TimeDelta( + radar_grid.sensing_stop)) # Populate the 'identification' group. logger.info("Populate identification metadata in output HDF5 file") - instrument_group = hdf5_file.create_group(f"/science/{radar_band}SAR") - identification_group = instrument_group.create_group("identification") - bounding_polygon = make_geo_grid_bounding_polygon(geo_grid, dem=dem) + instrument_group = \ + hdf5_file.create_group(f"/science/{radar_band}SAR") + identification_group = \ + instrument_group.create_group("identification") + bounding_polygon = make_geo_grid_bounding_polygon(geo_grid, + dem=dem_interp) populate_identification_group( identification_group=identification_group, product_spec=product_spec, @@ -304,20 +491,25 @@ def run_static_layers_workflow(config_file: os.PathLike | str) -> None: dem_description = get_raster_dataset_metadata_item( dem_raster_file, "dem_description", default="(NOT SPECIFIED)" ) - water_mask_description = get_raster_dataset_metadata_item( - water_mask_raster_file, - "water_mask_description", - default="(NOT SPECIFIED)", - ) + if flag_save_water_mask: + water_mask_description = get_raster_dataset_metadata_item( + water_mask_raster_file, + "water_mask_description", + default="(NOT SPECIFIED)", + ) + else: + water_mask_description = "(NOT APPLICABLE)" # Populate the 'grids' group. - logger.info("Populate raster layers and grid coordinates in output HDF5") + logger.info("Populate raster layers and grid coordinates in output" + " HDF5") grids_group = instrument_group.create_group("STATIC/grids") dataset_creation_kwds = build_hdf5_dataset_creation_kwds_dict( dataset_shape=(geo_grid.length, geo_grid.width), **output_params["dataset"] ) - with log_elapsed_time(logger.info, "Writing raster layers to output HDF5"): + with log_elapsed_time(logger.info, "Writing raster layers to" + " output HDF5"): populate_grids_group( grids_group=grids_group, product_spec=product_spec, @@ -360,18 +552,20 @@ def main(args: Sequence[str] | None = None) -> None: Parameters ---------- args : sequence of str or None, optional - The list of arguments. If None, the argument list is taken from `sys.argv`. - Defaults to None. + The list of arguments. If None, the argument list is taken from + `sys.argv`. Defaults to None. """ # Setup the argument parser. - parser = argparse.ArgumentParser(description="Run the NISAR Static Layers workflow") + parser = argparse.ArgumentParser( + description="Run the NISAR Static Layers workflow") parser.add_argument( "config_file", type=Path, help="Run configuration YAML file for the STATIC workflow", ) - # Parse the arguments and convert the result to a dict of keyword arguments. + # Parse the arguments and convert the result to a dict of keyword + # arguments. kwargs = vars(parser.parse_args(args)) # Run the workflow with the unpacked keyword arguments. diff --git a/share/nisar/defaults/static.yaml b/share/nisar/defaults/static.yaml index 83077376f..25f1768bf 100644 --- a/share/nisar/defaults/static.yaml +++ b/share/nisar/defaults/static.yaml @@ -2,6 +2,9 @@ runconfig: groups: dynamic_ancillary_file_group: + + # [OPTIONAL] One NISAR L1 RSLC formatted HDF5 file + input_file_path: # [REQUIRED] File path or URL of the input Digital Elevation Model (DEM) raster # file. # Must be an existing file in a GDAL-compatible raster format that spans the @@ -14,12 +17,14 @@ runconfig: # region of interest and conforms to the NISAR Water Mask Product Specification # (JPL D-107710). water_mask_raster_file: - # [REQUIRED] Path to the input orbit ephemeris XML file. + # [REQUIRED] if `input_file_path` is not provided, otherwise [OPTIONAL]. + # Path to the input orbit ephemeris XML file. # Must be an existing XML file conforming to the NISAR Orbit Ephemeris Product # Specification (JPL D-102253) and spanning the desired radar observation time # interval. orbit_xml_file: - # [REQUIRED] Path to the input radar pointing XML file. + # [REQUIRED] if `input_file_path` is not provided, otherwise [OPTIONAL]. + # Path to the input radar pointing XML file. # Must be an existing XML file conforming to the NISAR Radar Pointing Product # Specification (JPL D-102264) and spanning the desired radar observation time # interval. @@ -225,6 +230,16 @@ runconfig: # Defaults to 0.24. wavelength: 0.24 spacing: + # [OPTIONAL] Azimuth interval, in seconds. Equivalent to the Pulse + # Repetition Interval (PRI). + # If not provided, it will be inferred from the specified geographic + # grid. If provided, must be a positive value. + az_spacing: + # [OPTIONAL] Slant-range spacing, in meters, of the radar grid. + # Must be a positive value. + # If not provided, it will be inferred from the specified geographic + # grid. If provided, must be a positive value. + rg_spacing: # [OPTIONAL] Side length of the NxN grid of samples used to estimate the # required radar grid pixel spacing necessary to avoid undersampling the # output geocoded grid. @@ -232,6 +247,31 @@ runconfig: # Defaults to 5. pts_per_side: 5 bounding_box: + # [OPTIONAL] Azimuth start UTC date and time of the start of the radar + # observation, as a string in ISO 8601 format with up to nanosecond precision. + # If not provided, it will be inferred from the specified geographic grid. + # If provided, must be a valid azimuth time within the observation time + # interval. + start_time: + + # [OPTIONAL] Azimuth ending UTC date and time of the end of the radar + # observation, as a string in ISO 8601 + # format with up to nanosecond precision. Must be >= `start_time`. + # If not provided, it will be inferred from the specified geographic grid. + # If provided, must be a valid azimuth time within the observation time + # interval. + end_time: + + # [OPTIONAL] Starting range, in meters, of the first column of the radar grid. + # If not provided, it will be inferred from the specified geographic grid. + # If provided, must be a positive value. + start_range: + + # [OPTIONAL] Ending range, in meters, of the last column of the radar grid. + # If not provided, it will be inferred from the specified geographic grid. + # If provided, must be a positive value. + end_range: + # [OPTIONAL] Lower bound on the height of targets within the region of interest, # in meters above the reference ellipsoid of the output product. # Used to estimate the bounds of a radar grid that spans the region of interest. @@ -399,6 +439,10 @@ runconfig: extraiter: 10 output: + layers: + # [OPTIONAL] Whether to save the water mask layer in the output product. + # Defaults to true. + save_water_mask: true dataset: # [OPTIONAL] Chunk dimensions of 2-D raster datasets in the output product. # Setting `chunk_size` to [-1, -1] will disable chunked storage. diff --git a/share/nisar/schemas/static.yaml b/share/nisar/schemas/static.yaml index 7d20a5670..5f29c5674 100644 --- a/share/nisar/schemas/static.yaml +++ b/share/nisar/schemas/static.yaml @@ -1,10 +1,11 @@ runconfig: groups: dynamic_ancillary_file_group: + input_file_path: str(required=False) dem_raster_file: str() - water_mask_raster_file: str() - orbit_xml_file: str() - pointing_xml_file: str() + water_mask_raster_file: str(required=False) + orbit_xml_file: str(required=False) + pointing_xml_file: str(required=False) product_path_group: include('product_path_group_options', required=False) primary_executable: include('primary_executable_options', required=False) geometry: include('geometry_options', required=False) @@ -94,9 +95,23 @@ radar_grid_options: bounding_box: include('radar_grid_bounding_box_options', required=False) radar_grid_spacing_options: + az_spacing: num(min=0.0, required=False) + rg_spacing: num(min=0.0, required=False) pts_per_side: int(min=2, required=False) radar_grid_bounding_box_options: + start_time: any( + timestamp(), + regex(r'^\d{4}-\d{2}-\d{2}[T ]\d{2}:\d{2}:\d{2}(\.\d{1,9})?$'), + required=False + ) + end_time: any( + timestamp(), + regex(r'^\d{4}-\d{2}-\d{2}[T ]\d{2}:\d{2}:\d{2}(\.\d{1,9})?$'), + required=False + ) + start_range: num(required=False) + end_range: num(required=False) min_height: num(required=False) max_height: num(required=False) pts_per_edge: int(min=2, required=False) @@ -138,9 +153,13 @@ rdr2geo_options: extraiter: int(min=1, required=False) output_options: + layers: include('output_layers_options', required=False) dataset: include('dataset_options', required=False) file: include('file_options', required=False) +output_layers_options: + save_water_mask: bool(required=False) + dataset_options: chunk_size: list(int(min=-1), min=2, max=2, required=False) compression_enabled: bool(required=False)