diff --git a/cxx/isce3/geometry/RTC.cpp b/cxx/isce3/geometry/RTC.cpp index 80db426f1..5a65e7e9e 100644 --- a/cxx/isce3/geometry/RTC.cpp +++ b/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); diff --git a/cxx/isce3/geometry/loadDem.cpp b/cxx/isce3/geometry/loadDem.cpp index 4ad801508..44eb57001 100644 --- a/cxx/isce3/geometry/loadDem.cpp +++ b/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 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 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(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::max(); + max_x = std::numeric_limits::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; diff --git a/cxx/isce3/geometry/loadDem.h b/cxx/isce3/geometry/loadDem.h index eb35139c8..ce95d1809 100644 --- a/cxx/isce3/geometry/loadDem.h +++ b/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 diff --git a/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp b/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp index 47167b30f..b58c8988e 100644 --- a/python/extensions/pybind_isce3/geometry/DEMInterpolator.cpp +++ b/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 -------