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,14 +132,13 @@ 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,
     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()) {
@@ -152,110 +151,62 @@ isce3::error::ErrorCode loadDemFromProj(
     } 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 xMid = x0 + t * (xf - x0);
+            double yMid = minY + t * (maxY - minY);
+
+            // West edge (x = x0, y varies)
+            auto west_llh = proj->inverse({x0, yMid, 0});
+            auto west_xy  = dem_proj->forward(west_llh);
+            all_x.push_back(west_xy[0]);
+            all_y.push_back(west_xy[1]);
+
+            // East edge (x = xf, y varies)
+            auto east_llh = proj->inverse({xf, yMid, 0});
+            auto east_xy  = dem_proj->forward(east_llh);
+            all_x.push_back(east_xy[0]);
+            all_y.push_back(east_xy[1]);
+
+            // South edge (y = minY, x varies)
+            auto south_llh = proj->inverse({xMid, minY, 0});
+            auto south_xy  = dem_proj->forward(south_llh);
+            all_x.push_back(south_xy[0]);
+            all_y.push_back(south_xy[1]);
+
+            // North edge (y = maxY, x varies)
+            auto north_llh = proj->inverse({xMid, maxY, 0});
+            auto north_xy  = dem_proj->forward(north_llh);
+            all_x.push_back(north_xy[0]);
+            all_y.push_back(north_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);
             }
         }
     }
@@ -264,7 +215,6 @@ isce3::error::ErrorCode loadDemFromProj(
     min_y -= margin_y;
     max_y += margin_y;
 
-
     float margin_x = dem_margin_x_in_pixels * dem_raster.dx();
     min_x -= margin_x;
     max_x += margin_x;

cxx/isce3/geometry/loadDem.h

@@ -47,26 +47,23 @@ 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                      Minimum X/easting position of the input
+# bounding box in the coordinate system of `proj`
+* @param[in]  xf                      Maximum X/easting position of the input
+# bounding box in the coordinate system of `proj`
+* @param[in]  minY                    Minimum Y/northing position of the input
+# bounding box in the coordinate system of `proj`
+* @param[in]  maxY                    Maximum 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,
@@ -75,7 +72,8 @@ isce3::error::ErrorCode loadDemFromProj(
     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

@@ -232,7 +232,8 @@ void addbinding_load_dem_from_proj(py::module& m)
             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);
 
@@ -245,7 +246,8 @@ void addbinding_load_dem_from_proj(py::module& m)
                                              proj,
                                              dem_margin_x_in_pixels,
                                              dem_margin_y_in_pixels,
-                                             dem_raster_band);
+                                             dem_raster_band,
+                                             n_edge_samples);
 
             return dem_interp;
         },
@@ -259,6 +261,7 @@ void addbinding_load_dem_from_proj(py::module& m)
         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
     -------