Improve DEM bounding box selection when product coordinates are in a different projection than the DEM

cxx/isce3/geometry/RTC.cpp

@@ -1174,7 +1174,7 @@ void _RunBlock(const int jmax, const int block_size,
     the bounding-box corners.
     */
     int dem_margin_x_in_pixels = 100;
-    int dem_margin_y_in_pixels = 200;
+    int dem_margin_y_in_pixels = 100;
     auto error_code = loadDemFromProj(
             dem_raster, minX, maxX, minY, maxY, &dem_interp_block, proj,
             dem_margin_x_in_pixels, dem_margin_y_in_pixels);

cxx/isce3/geometry/loadDem.cpp

@@ -132,143 +132,111 @@ DEMInterpolator DEMRasterToInterpolator(
     return demInterp;
 }
 
-
 isce3::error::ErrorCode loadDemFromProj(
     isce3::io::Raster& dem_raster, const double x0, const double xf,
-    const double minY, const double maxY,
+    const double y0, const double yf,
     DEMInterpolator* dem_interp,
     isce3::core::ProjectionBase* proj, const int dem_margin_x_in_pixels,
-    const int dem_margin_y_in_pixels, const int dem_raster_band) {
-
+    const int dem_margin_y_in_pixels, const int dem_raster_band,
+    const int n_edge_samples){
     double min_x, max_x, min_y, max_y;
 
-    if (proj == nullptr || proj->code() == dem_raster.getEPSG()) {
+    min_y = std::min(y0, yf);
+    max_y = std::max(y0, yf);
+
+    min_x = std::min(x0, xf);
+    max_x = std::max(x0, xf);
+
+    // Test for antimeridian crossing (EPSG 4326) and unwrap X (longitude)
+    // coordinates if `max_x - min_x` is greater than 180 degrees
+    if ((((proj == nullptr) && (dem_raster.getEPSG() == 4326)) ||
+         ((proj != nullptr) && (proj->code() == 4326))) &&
+        (max_x - min_x > 180)) {
+            const double x0_unwrapped = x0 < 0 ? x0 + 360.0 : x0;
+            const double xf_unwrapped = xf < 0 ? xf + 360.0 : xf;
+            min_x = std::min(x0_unwrapped, xf_unwrapped);
+            max_x = std::max(x0_unwrapped, xf_unwrapped);
+    }
 
-        min_x = x0;
-        max_x = xf;
-        min_y = std::min(minY, maxY);
-        max_y = std::max(minY, maxY);
+    if (proj != nullptr && proj->code() != dem_raster.getEPSG()) {
 
-    } else {
         std::unique_ptr<isce3::core::ProjectionBase> dem_proj(
                 isce3::core::createProj(dem_raster.getEPSG()));
-        auto p_west_1_llh = proj->inverse({x0, minY, 0});
-        auto p_west_2_llh = proj->inverse({x0, maxY, 0});
-        auto p_east_1_llh = proj->inverse({xf, minY, 0});
-        auto p_east_2_llh = proj->inverse({xf, maxY, 0});
-
-        auto p_west_1_xy = dem_proj->forward(p_west_1_llh);
-        auto p_west_2_xy = dem_proj->forward(p_west_2_llh);
-        auto p_east_1_xy = dem_proj->forward(p_east_1_llh);
-        auto p_east_2_xy = dem_proj->forward(p_east_2_llh);
-
-        min_y = std::min(std::min(p_west_1_xy[1], p_west_2_xy[1]),
-                         std::min(p_east_1_xy[1], p_east_2_xy[1]));
-        max_y = std::max(std::max(p_west_1_xy[1], p_west_2_xy[1]),
-                         std::max(p_east_1_xy[1], p_east_2_xy[1]));
-
-        /* We address two cases in this if statement below:
-           1. If the DEM projection is NOT geographic:
-              No antimeridian crossing, compute `min_x` and `max_x`
-              directly
-           2. The user projection is in polar stereographic AND 
-              the DEM projection is geographic:
-              In this case we need to check for antimeridian crossing.
-        */
-        if (dem_raster.getEPSG() != 4326 or proj->code() == 3031 or
-                proj->code() == 3413) {
-
-            // Compute X min/max directly
-            min_x = std::min(std::min(p_west_1_xy[0], p_west_2_xy[0]),
-                             std::min(p_east_1_xy[0], p_east_2_xy[0]));
-            max_x = std::max(std::max(p_west_1_xy[0], p_west_2_xy[0]),
-                             std::max(p_east_1_xy[0], p_east_2_xy[0]));
-
-            if (dem_raster.getEPSG() == 4326 and
-                    max_x - min_x > 180 and
-                    (proj->code() == 3031 or proj->code() == 3413)) {
-
-                /*
-                If (DEM is in geographic (EPSG: 4326) and
-                the difference between max and min longitudes is greater
-                than 180 and the map grid is in polar stereo (i.e., proj
-                epsg == 3031 or 3413), we cannot assume that `x0` is at
-                the western side of `xf`.
-                In that case, we also compute the (min/max using longitudes in
-                the [0, 360] range */
-
-                /* The conversion of longitude values from the [-180, 180]
-                domain to the [0, 360] domain is done by adding 360 to
-                negative longitude values. */
-                const double p1_0_360 = \
-                    p_west_1_xy[0] < 0 ? p_west_1_xy[0] + 360 : p_west_1_xy[0];
-                const double p2_0_360 = \
-                    p_west_2_xy[0] < 0 ? p_west_2_xy[0] + 360 : p_west_2_xy[0];
-                const double p3_0_360 = \
-                    p_east_1_xy[0] < 0 ? p_east_1_xy[0] + 360 : p_east_1_xy[0];
-                const double p4_0_360 = \
-                    p_east_2_xy[0] < 0 ? p_east_2_xy[0] + 360 : p_east_2_xy[0];
-
-                // Compute min/max longitudes in the [0, 360] domain
-                min_x = std::min(std::min(p1_0_360, p2_0_360),
-                                 std::min(p3_0_360, p4_0_360));
-                max_x = std::max(std::max(p1_0_360, p2_0_360),
-                                 std::max(p3_0_360, p4_0_360));
-                }
-
-       } else {
-            /*
-            X-coordinates may be wrapped due to the antimeridian
-            crossing. In this case, we compute western and eastern boundaries
-            separately.
-            We just need to make sure that there's no antimeridian crossing
-            in between the western and eastern edges
-            */
-
-            // Western edge
-            if (std::abs(p_west_1_xy[0] - p_west_2_xy[0]) < 180) {
-
-                // Normal case
-                min_x = std::min(p_west_1_xy[0], p_west_2_xy[0]);
-            }
-            else {
-
-                // Antimeridian crossing
-                const double p1_0_360 = \
-                    p_west_1_xy[0] < 0 ? p_west_1_xy[0] + 360 : p_west_1_xy[0];
-                const double p2_0_360 = \
-                    p_west_2_xy[0] < 0 ? p_west_2_xy[0] + 360 : p_west_2_xy[0];
-                min_x = std::min(p1_0_360, p2_0_360);
-            }
 
-            // Eastern edge
-            if (std::abs(p_east_1_xy[0] - p_east_2_xy[0]) < 180) {
-
-                // Normal case
-                max_x = std::max(p_east_1_xy[0], p_east_2_xy[0]);
-            }
-
-            else {
+        const int N = std::max(n_edge_samples, 2);
+
+        // Densely sample all four edges of the input bounding box
+        // to capture curvature introduced by reprojection
+        // (e.g. UTM -> geographic).
+        std::vector<double> all_x, all_y;
+        all_x.reserve(4 * N);
+        all_y.reserve(4 * N);
+
+        for (int i = 0; i < N; ++i) {
+            double t = static_cast<double>(i) / (N - 1);
+            double x_mid = min_x + t * (max_x - min_x);
+            double y_mid = min_y + t * (max_y - min_y);
+
+            // Left edge (x = min_x, y varies)
+            auto left_llh = proj->inverse({min_x, y_mid, 0});
+            auto left_xy  = dem_proj->forward(left_llh);
+            all_x.push_back(left_xy[0]);
+            all_y.push_back(left_xy[1]);
+
+            // Right edge (x = max_x, y varies)
+            auto right_llh = proj->inverse({max_x, y_mid, 0});
+            auto right_xy  = dem_proj->forward(right_llh);
+            all_x.push_back(right_xy[0]);
+            all_y.push_back(right_xy[1]);
+
+            // Bottom edge (y = min_y, x varies)
+            auto bottom_llh = proj->inverse({x_mid, min_y, 0});
+            auto bottom_xy  = dem_proj->forward(bottom_llh);
+            all_x.push_back(bottom_xy[0]);
+            all_y.push_back(bottom_xy[1]);
+
+            // Top edge (y = max_y, x varies)
+            auto top_llh = proj->inverse({x_mid, max_y, 0});
+            auto top_xy  = dem_proj->forward(top_llh);
+            all_x.push_back(top_xy[0]);
+            all_y.push_back(top_xy[1]);
+        }
 
-                // Antimeridian crossing
-                const double p3_0_360 = \
-                    p_east_1_xy[0] < 0 ? p_east_1_xy[0] + 360 : p_east_1_xy[0];
-                const double p4_0_360 = \
-                    p_east_2_xy[0] < 0 ? p_east_2_xy[0] + 360 : p_east_2_xy[0];
-                max_x = std::max(p3_0_360, p4_0_360);
+        min_y = *std::min_element(all_y.begin(), all_y.end());
+        max_y = *std::max_element(all_y.begin(), all_y.end());
+
+        min_x = *std::min_element(all_x.begin(), all_x.end());
+        max_x = *std::max_element(all_x.begin(), all_x.end());
+
+        // If the DEM is in geographic coordinates and the X range
+        // exceeds 180 degrees, an antimeridian crossing is likely.
+        // Retry in [0, 360] domain.
+        if (dem_raster.getEPSG() == 4326 && (max_x - min_x) > 180.0) {
+            min_x = std::numeric_limits<double>::max();
+            max_x = std::numeric_limits<double>::lowest();
+            for (double lon : all_x) {
+                double lon_360 = lon < 0 ? lon + 360.0 : lon;
+                min_x = std::min(min_x, lon_360);
+                max_x = std::max(max_x, lon_360);
             }
         }
     }
 
-    float margin_y = dem_margin_y_in_pixels * std::abs(dem_raster.dy());
+    double margin_y = dem_margin_y_in_pixels * std::abs(dem_raster.dy());
     min_y -= margin_y;
     max_y += margin_y;
 
-
-    float margin_x = dem_margin_x_in_pixels * dem_raster.dx();
+    double margin_x = dem_margin_x_in_pixels * std::abs(dem_raster.dx());
     min_x -= margin_x;
     max_x += margin_x;
 
+    // If DEM coordinates are in geographic, ensure latitude values
+    // fall between [-90.0, 90.0] after applying `margin_y`
+    if (dem_raster.getEPSG() == 4326) {
+        min_y = std::clamp(min_y, -90.0, 90.0);
+        max_y = std::clamp(max_y, -90.0, 90.0);
+    }
+
     isce3::error::ErrorCode error_code;
     _Pragma("omp critical")
     {

cxx/isce3/geometry/loadDem.h

@@ -47,35 +47,33 @@ isce3::geometry::DEMInterpolator DEMRasterToInterpolator(
 * in the same or different coordinate system as the DEM raster
 *
 * @param[in]  dem_raster              DEM raster
-* @param[in]  x0                      Easting/longitude of western edge of bounding box,
-* 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)
-* @param[in]  xf                      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)
-* @param[in]  minY                    Minimum Y/northing position
-* @param[in]  maxY                    Maximum Y/northing position
+* @param[in]  x0                      Starting X/easting position of the input
+* bounding box in the coordinate system of `proj`
+* @param[in]  xf                      Ending X/easting position of the input
+* bounding box in the coordinate system of `proj`
+* @param[in]  y0                      Starting Y/northing position of the input
+* bounding box in the coordinate system of `proj`
+* @param[in]  yf                      Ending Y/northing position of the input
+* bounding box in the coordinate system of `proj`
 * @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
 * @param[in]  dem_margin_y_in_pixels  DEM Y/northing margin in pixels
 * @param[in]  dem_raster_band         DEM raster band (starting from 1)
+* @param[in]  n_edge_samples          Number of points sampled along each edge
+* of the bounding box when reprojecting. Must be >= 2 to include corners; 
+* values below 2 are clamped to 2. Default: 11.
 */
 isce3::error::ErrorCode loadDemFromProj(
     isce3::io::Raster& dem_raster,
-    const double minX, const double maxX, const double minY,
-    const double maxY, isce3::geometry::DEMInterpolator* dem_interp,
+    const double x0, const double xf, const double y0,
+    const double yf, isce3::geometry::DEMInterpolator* dem_interp,
     isce3::core::ProjectionBase* proj = nullptr,
     const int dem_margin_x_in_pixels = 100,
     const int dem_margin_y_in_pixels = 100,
-    const int dem_raster_band = 1);
+    const int dem_raster_band = 1,
+    const int n_edge_samples = 11);
 
 /*
  Interpolate DEM at position (x, y) considering that input_proj and

python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp

@@ -226,39 +226,42 @@ void addbinding_load_dem_from_proj(py::module& m)
         [](isce3::io::Raster &dem_raster,
             const double x0,
             const double xf,
-            const double minY,
-            const double maxY,
+            const double y0,
+            const double yf,
             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) {
+            const int dem_raster_band,
+            const int n_edge_samples) {
 
             DEMInterp dem_interp(0, dem_interp_method);
 
             isce3::geometry::loadDemFromProj(dem_raster,
                                              x0,
                                              xf,
-                                             minY,
-                                             maxY,
+                                             y0,
+                                             yf,
                                              &dem_interp,
                                              proj,
                                              dem_margin_x_in_pixels,
                                              dem_margin_y_in_pixels,
-                                             dem_raster_band);
+                                             dem_raster_band,
+                                             n_edge_samples);
 
             return dem_interp;
         },
         py::arg("dem_raster"),
         py::arg("x0"),
         py::arg("xf"),
-        py::arg("min_y"),
-        py::arg("max_y"),
+        py::arg("y0"),
+        py::arg("yf"),
         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,
+        py::arg("n_edge_samples") = 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
@@ -268,22 +271,17 @@ void addbinding_load_dem_from_proj(py::module& m)
     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        
+        Minimum X/easting position of the input bounding box in the coordinate
+        system of `proj`
     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)
+        Maximum X/easting position of the input bounding box in the coordinate
+        system of `proj`
     min_y: double
-        Minimum Y/northing position
+        Minimum Y/northing position of the input bounding box in the coordinate
+        system of `proj`
     max_y: double
-        Maximum Y/northing position
+        Maximum Y/northing position of the input bounding box in the coordinate
+        system of `proj`
     dem_interp_method: isce3.core.DataInterpMethod
         DEM interpolation method
     proj: 
@@ -294,6 +292,10 @@ void addbinding_load_dem_from_proj(py::module& m)
         DEM Y/northing margin in pixels
     dem_raster_band: int
         DEM raster band (starting from 1)
+    n_edge_samples: int
+        Number of points sampled along each edge of the bounding box when
+        reprojecting. Must be >= 2 to include corners.  Values below 2
+        are clamped to 2. Default: 11.
 
     Returns
     -------