Enable the use of an optional RSLC product to generate static layers (WIP)

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

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.
+        )");
+}

python/extensions/pybind_isce3/geometry/DEMInterpolator.h

@@ -5,3 +5,4 @@
 
 void addbinding(pybind11::class_<isce3::geometry::DEMInterpolator>&);
 void addbinding_DEM_raster2interpolator(pybind11::module&);
+void addbinding_load_dem_from_proj(pybind11::module&);

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);
 }

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

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}]"

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)