Add proj_trans_bounds to compute the image of a input bounding box through a transformation (#2882)

Fixes #2779

7 files changed, 1089 insertions, 3 deletions

diff --git a/Doxyfile b/Doxyfile
index 37776f9f..accc13f0 100644
--- a/Doxyfile
+++ b/Doxyfile
@@ -766,7 +766,7 @@ WARN_LOGFILE           =
 # spaces. See also FILE_PATTERNS and EXTENSION_MAPPING
 # Note: If this tag is empty the current directory is searched.
 
-INPUT                  = src/iso19111 src/iso19111/operation include/proj src/proj.h src/proj_experimental.h src/general_doc.dox src/filemanager.cpp src/networkfilemanager.cpp
+INPUT                  = src/iso19111 src/iso19111/operation include/proj src/proj.h src/proj_experimental.h src/general_doc.dox src/filemanager.cpp src/networkfilemanager.cpp src/4D_api.cpp
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/docs/source/development/reference/functions.rst b/docs/source/development/reference/functions.rst
index bf5bb94b..564135d2 100644
--- a/docs/source/development/reference/functions.rst
+++ b/docs/source/development/reference/functions.rst
@@ -391,6 +391,11 @@ Coordinate transformation
               reasons.
 
 
+
+.. doxygenfunction:: proj_trans_bounds
+   :project: doxygen_api
+
+
 Error reporting
 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
diff --git a/scripts/doxygen.sh b/scripts/doxygen.sh
index 517447d2..7de60d52 100755
--- a/scripts/doxygen.sh
+++ b/scripts/doxygen.sh
@@ -35,7 +35,7 @@ fi
 
 rm -rf docs/build/xml/
 
-(cat Doxyfile; printf "GENERATE_HTML=NO\nGENERATE_XML=YES\nINPUT= src/iso19111 src/iso19111/operation include/proj src/proj.h src/filemanager.cpp src/networkfilemanager.cpp src/general_doc.dox") | doxygen -  > docs/build/docs_log.txt 2>&1
+(cat Doxyfile; printf "GENERATE_HTML=NO\nGENERATE_XML=YES\nINPUT= src/iso19111 src/iso19111/operation include/proj src/proj.h src/filemanager.cpp src/networkfilemanager.cpp src/4D_api.cpp src/general_doc.dox") | doxygen -  > docs/build/docs_log.txt 2>&1
 if grep -i warning docs/build/docs_log.txt; then
     echo "Doxygen warnings found" && cat docs/build/docs_log.txt && /bin/false;
 else
diff --git a/scripts/reference_exported_symbols.txt b/scripts/reference_exported_symbols.txt
index 682f5d39..250ac793 100644
--- a/scripts/reference_exported_symbols.txt
+++ b/scripts/reference_exported_symbols.txt
@@ -1048,6 +1048,7 @@ proj_todeg
 proj_torad
 proj_trans
 proj_trans_array
+proj_trans_bounds
 proj_trans_generic
 proj_unit_list_destroy
 proj_uom_get_info_from_database
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index 132c21d3..6f8c3027 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -966,6 +966,534 @@ std::string pj_add_type_crs_if_needed(const std::string& str)
     return ret;
 }
 
+
+// ---------------------------------------------------------------------------
+static double simple_min(const double* data, const int arr_len) {
+    double min_value = data[0];
+    for( int iii = 1; iii < arr_len; iii++ ) {
+        if (data[iii] < min_value)
+            min_value = data[iii];
+    }
+    return min_value;
+}
+
+
+// ---------------------------------------------------------------------------
+static double simple_max(const double* data, const int arr_len) {
+    double max_value = data[0];
+    for( int iii = 1; iii < arr_len; iii++ ) {
+        if ((data[iii] > max_value || max_value == HUGE_VAL) && data[iii] != HUGE_VAL)
+            max_value = data[iii];
+    }
+    return max_value;
+ }
+
+
+// ---------------------------------------------------------------------------
+static int _find_previous_index(const int iii, const double* data, const int arr_len) {
+    // find index of nearest valid previous value if exists
+    int prev_iii = iii - 1;
+    if (prev_iii == -1)  // handle wraparound
+        prev_iii = arr_len - 1;
+    while (data[prev_iii] == HUGE_VAL && prev_iii != iii) {
+        prev_iii --;
+        if (prev_iii == -1)  // handle wraparound
+            prev_iii = arr_len - 1;
+    }
+    return prev_iii;
+}
+
+
+// ---------------------------------------------------------------------------
+/******************************************************************************
+Handles the case when longitude values cross the antimeridian
+when calculating the minimum.
+Note: The data array must be in a linear ring.
+Note: This requires a densified ring with at least 2 additional
+        points per edge to correctly handle global extents.
+If only 1 additional point:
+    |        |
+    |RL--x0--|RL--
+    |        |
+-180    180|-180
+If they are evenly spaced and it crosses the antimeridian:
+x0 - L = 180
+R - x0 = -180
+For example:
+Let R = -179.9, x0 = 0.1, L = -179.89
+x0 - L = 0.1 - -179.9 = 180
+R - x0 = -179.89 - 0.1 ~= -180
+This is the same in the case when it didn't cross the antimeridian.
+If you have 2 additional points:
+    |            |
+    |RL--x0--x1--|RL--
+    |            |
+-180        180|-180
+If they are evenly spaced and it crosses the antimeridian:
+x0 - L = 120
+x1 - x0 = 120
+R - x1 = -240
+For example:
+Let R = -179.9, x0 = -59.9, x1 = 60.1 L = -179.89
+x0 - L = 59.9 - -179.9 = 120
+x1 - x0 = 60.1 - 59.9 = 120
+R - x1 = -179.89 - 60.1 ~= -240
+However, if they are evenly spaced and it didn't cross the antimeridian:
+x0 - L = 120
+x1 - x0 = 120
+R - x1 = 120
+From this, we have a delta that is guaranteed to be significantly
+large enough to tell the difference reguarless of the direction
+the antimeridian was crossed.
+However, even though the spacing was even in the source projection, it isn't
+guaranteed in the target geographic projection. So, instead of 240, 200 is used
+as it significantly larger than 120 to be sure that the antimeridian was crossed
+but smalller than 240 to account for possible irregularities in distances
+when re-projecting. Also, 200 ensures latitudes are ignored for axis order handling.
+******************************************************************************/
+static double antimeridian_min(const double* data, const int arr_len) {
+    double positive_min = HUGE_VAL;
+    double min_value = HUGE_VAL;
+    int crossed_meridian_count = 0;
+    bool positive_meridian = false;
+
+    for( int iii = 0; iii < arr_len; iii++ ) {
+        if (data[iii] == HUGE_VAL)
+            continue;
+        int prev_iii = _find_previous_index(iii, data, arr_len);
+        // check if crossed meridian
+        double delta = data[prev_iii] - data[iii];
+        // 180 -> -180
+        if (delta >= 200 && delta != HUGE_VAL) {
+            if (crossed_meridian_count == 0)
+                positive_min = min_value;
+            crossed_meridian_count ++;
+            positive_meridian = false;
+        // -180 -> 180
+        } else if (delta <= -200 && delta != HUGE_VAL) {
+            if (crossed_meridian_count == 0)
+                positive_min = data[iii];
+            crossed_meridian_count ++;
+            positive_meridian = true;
+        }
+        // positive meridian side min
+        if (positive_meridian && data[iii] < positive_min)
+            positive_min = data[iii];
+        // track general min value
+        if (data[iii] < min_value)
+            min_value = data[iii];
+    }
+
+    if (crossed_meridian_count == 2)
+        return positive_min;
+    else if (crossed_meridian_count == 4)
+        // bounds extends beyond -180/180
+        return -180;
+    return min_value;
+}
+
+
+// ---------------------------------------------------------------------------
+// Handles the case when longitude values cross the antimeridian
+// when calculating the minimum.
+// Note: The data array must be in a linear ring.
+// Note: This requires a densified ring with at least 2 additional
+//       points per edge to correctly handle global extents.
+// See antimeridian_min docstring for reasoning.
+static double antimeridian_max(const double* data, const int arr_len) {
+    double negative_max = -HUGE_VAL;
+    double max_value = -HUGE_VAL;
+    bool negative_meridian = false;
+    int crossed_meridian_count = 0;
+
+    for( int iii = 0; iii < arr_len; iii++ ) {
+        if (data[iii] == HUGE_VAL)
+            continue;
+        int prev_iii = _find_previous_index(iii, data, arr_len);
+        // check if crossed meridian
+        double delta = data[prev_iii] - data[iii];
+        // 180 -> -180
+        if (delta >= 200 && delta != HUGE_VAL) {
+            if (crossed_meridian_count == 0)
+                negative_max = data[iii];
+            crossed_meridian_count ++;
+            negative_meridian = true;
+        // -180 -> 180
+        } else if (delta <= -200 && delta != HUGE_VAL){
+            if (crossed_meridian_count == 0)
+                negative_max = max_value;
+            negative_meridian = false;
+            crossed_meridian_count++;
+        }
+        // negative meridian side max
+        if (negative_meridian
+            && (data[iii] > negative_max || negative_max == HUGE_VAL)
+            && data[iii] != HUGE_VAL
+        )
+            negative_max = data[iii];
+        // track general max value
+        if ((data[iii] > max_value || max_value == HUGE_VAL) && data[iii] != HUGE_VAL)
+            max_value = data[iii];
+    }
+    if (crossed_meridian_count == 2)
+        return negative_max;
+    else if (crossed_meridian_count == 4)
+        // bounds extends beyond -180/180
+        return 180;
+    return max_value;
+}
+
+
+// ---------------------------------------------------------------------------
+// Check if the original projected bounds contains
+// the north pole.
+// This assumes that the destination CRS is geographic.
+static bool contains_north_pole(
+    PJ* projobj,
+    PJ_DIRECTION pj_direction,
+    const double xmin,
+    const double ymin,
+    const double xmax,
+    const double ymax,
+    bool lon_lat_order
+) {
+    double pole_y = 90;
+    double pole_x = 0;
+    if (!lon_lat_order) {
+        pole_y = 0;
+        pole_x = 90;
+    }
+    proj_trans_generic(
+        projobj,
+        opposite_direction(pj_direction),
+        &pole_x, sizeof(double), 1,
+        &pole_y, sizeof(double), 1,
+        nullptr, sizeof(double), 0,
+        nullptr, sizeof(double), 0
+    );
+    if (xmin < pole_x && pole_x < xmax && ymax > pole_y && pole_y > ymin)
+        return true;
+    return false;
+}
+
+
+// ---------------------------------------------------------------------------
+// Check if the original projected bounds contains
+// the south pole.
+// This assumes that the destination CRS is geographic.
+static bool contains_south_pole(
+    PJ* projobj,
+    PJ_DIRECTION pj_direction,
+    const double xmin,
+    const double ymin,
+    const double xmax,
+    const double ymax,
+    bool lon_lat_order
+) {
+    double pole_y = -90;
+    double pole_x = 0;
+    if (!lon_lat_order) {
+        pole_y = 0;
+        pole_x = -90;
+    }
+    proj_trans_generic(
+        projobj,
+        opposite_direction(pj_direction),
+        &pole_x, sizeof(double), 1,
+        &pole_y, sizeof(double), 1,
+        nullptr, sizeof(double), 0,
+        nullptr, sizeof(double), 0
+    );
+    if (xmin < pole_x && pole_x < xmax && ymax > pole_y && pole_y > ymin)
+        return true;
+    return false;
+}
+
+
+// ---------------------------------------------------------------------------
+// Check if the target CRS of the transformation
+// has the longitude latitude axis order.
+// This assumes that the destination CRS is geographic.
+static int target_crs_lon_lat_order(
+    PJ_CONTEXT* transformer_ctx,
+    PJ* transformer_pj,
+    PJ_DIRECTION pj_direction
+) {
+    PJ* target_crs = nullptr;
+    if (pj_direction == PJ_FWD)
+        target_crs = proj_get_target_crs(transformer_ctx, transformer_pj);
+    else if (pj_direction == PJ_INV)
+        target_crs = proj_get_source_crs(transformer_ctx, transformer_pj);
+    if (target_crs == nullptr) {
+        proj_context_log_debug(transformer_ctx, "Unable to retrieve target CRS");
+        return -1;
+    }
+    PJ* coord_system_pj = proj_crs_get_coordinate_system(
+        transformer_ctx,
+        target_crs
+    );
+    proj_destroy(target_crs);
+    if (coord_system_pj == nullptr) {
+        proj_context_log_debug(transformer_ctx, "Unable to get target CRS coordinate system.");
+        return -1;
+    }
+    const char* abbrev = nullptr;
+    int success = proj_cs_get_axis_info(
+        transformer_ctx,
+        coord_system_pj,
+        0,
+        nullptr,
+        &abbrev,
+        nullptr,
+        nullptr,
+        nullptr,
+        nullptr,
+        nullptr
+    );
+    proj_destroy(coord_system_pj);
+    if (success != 1)
+        return -1;
+    return strcmp(abbrev, "lon") == 0 || strcmp(abbrev, "Lon") == 0;
+}
+
+// ---------------------------------------------------------------------------
+
+/** \brief Transform boundary,
+ *
+ * Transform boundary densifying the edges to account for nonlinear
+ * transformations along these edges and extracting the outermost bounds.
+ *
+ * If the destination CRS is geographic, the first axis is longitude,
+ * and xmax < xmin then the bounds crossed the antimeridian.
+ * In this scenario there are two polygons, one on each side of the antimeridian.
+ * The first polygon should be constructed with (xmin, ymin, 180, ymax)
+ * and the second with (-180, ymin, xmax, ymax).
+ *
+ * If the destination CRS is geographic, the first axis is latitude,
+ * and ymax < ymin then the bounds crossed the antimeridian.
+ * In this scenario there are two polygons, one on each side of the antimeridian.
+ * The first polygon should be constructed with (ymin, xmin, ymax, 180)
+ * and the second with (ymin, -180, ymax, xmax).
+ *
+ * @param context The PJ_CONTEXT object.
+ * @param P The PJ object representing the transformation.
+ * @param direction The direction of the transformation.
+ * @param xmin Minimum bounding coordinate of the first axis in source CRS
+ *             (target CRS if direction is inverse).
+ * @param ymin Minimum bounding coordinate of the second axis in source CRS.
+ *             (target CRS if direction is inverse).
+ * @param xmax Maximum bounding coordinate of the first axis in source CRS.
+ *             (target CRS if direction is inverse).
+ * @param ymax Maximum bounding coordinate of the second axis in source CRS.
+ *             (target CRS if direction is inverse).
+ * @param out_xmin Minimum bounding coordinate of the first axis in target CRS
+ *             (source CRS if direction is inverse).
+ * @param out_ymin Minimum bounding coordinate of the second axis in target CRS.
+ *             (source CRS if direction is inverse).
+ * @param out_xmax Maximum bounding coordinate of the first axis in target CRS.
+ *             (source CRS if direction is inverse).
+ * @param out_ymax Maximum bounding coordinate of the second axis in target CRS.
+ *             (source CRS if direction is inverse).
+ * @param densify_pts Recommended to use 21. This is the number of points
+ *     to use to densify the bounding polygon in the transformation.
+ * @return an integer. 1 if successful. 0 if failures encountered.
+ * @since 8.2
+ */
+int proj_trans_bounds(PJ_CONTEXT* context,
+                      PJ *P,
+                      PJ_DIRECTION direction,
+                      const double xmin,
+                      const double ymin,
+                      const double xmax,
+                      const double ymax,
+                      double* out_xmin,
+                      double* out_ymin,
+                      double* out_xmax,
+                      double* out_ymax,
+                      int densify_pts
+) {
+    *out_xmin = HUGE_VAL;
+    *out_ymin = HUGE_VAL;
+    *out_xmax = HUGE_VAL;
+    *out_ymax = HUGE_VAL;
+
+    if (P == nullptr) {
+        proj_log_error(P, _("NULL P object not allowed."));
+        proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+        return false;
+    }
+    if (densify_pts < 0 || densify_pts > 10000) {
+        proj_log_error(P, _("densify_pts must be between 0-10000."));
+        proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+        return false;
+    }
+
+    PJ_PROJ_INFO pj_info = proj_pj_info(P);
+    if (pj_info.id == nullptr) {
+        proj_log_error(P, _("NULL transformation not allowed,"));
+        proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+        return false;
+    }
+    if (strcmp(pj_info.id, "noop") == 0 || direction == PJ_IDENT) {
+        *out_xmin = xmin;
+        *out_xmax = xmax;
+        *out_ymin = ymin;
+        *out_ymax = ymax;
+        return true;
+    }
+
+    bool degree_output = proj_degree_output(P, direction) != 0;
+    bool degree_input = proj_degree_input(P, direction) != 0;
+    if (degree_output && densify_pts < 2) {
+        proj_log_error(P, _("densify_pts must be at least 2 if the output is geograpic."));
+        proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+        return false;
+    }
+
+    int side_pts = densify_pts + 1;  // add one because we are densifying
+    const int boundary_len = side_pts * 4;
+    std::vector<double> x_boundary_array;
+    std::vector<double> y_boundary_array;
+    try
+    {
+        x_boundary_array.resize(boundary_len);
+        y_boundary_array.resize(boundary_len);
+    }
+    catch( const std::exception & e ) // memory allocation failure
+    {
+        proj_log_error(P, e.what());
+        proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+        return false;
+    }
+    double delta_x = 0;
+    double delta_y = 0;
+    bool north_pole_in_bounds = false;
+    bool south_pole_in_bounds = false;
+    bool input_lon_lat_order = false;
+    bool output_lon_lat_order = false;
+    if (degree_input) {
+        int in_order_lon_lat = target_crs_lon_lat_order(
+            context, P, opposite_direction(direction)
+        );
+        if (in_order_lon_lat == -1)
+            return false;
+        input_lon_lat_order = in_order_lon_lat != 0;
+    }
+    if (degree_output) {
+        int out_order_lon_lat = target_crs_lon_lat_order(context, P, direction);
+        if (out_order_lon_lat == -1)
+            return false;
+        output_lon_lat_order = out_order_lon_lat != 0;
+        north_pole_in_bounds = contains_north_pole(
+            P,
+            direction,
+            xmin,
+            ymin,
+            xmax,
+            ymax,
+            output_lon_lat_order
+        );
+        south_pole_in_bounds = contains_south_pole(
+            P,
+            direction,
+            xmin,
+            ymin,
+            xmax,
+            ymax,
+            output_lon_lat_order
+        );
+    }
+
+    if (degree_input && xmax < xmin) {
+        if (!input_lon_lat_order) {
+            proj_log_error(P, _("latitude max < latitude min."));
+            proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+            return false;
+        }
+        // handle antimeridian
+        delta_x = (xmax - xmin + 360.0) / side_pts;
+    } else {
+        delta_x = (xmax - xmin) / side_pts;
+    }
+    if (degree_input && ymax < ymin) {
+        if (input_lon_lat_order) {
+            proj_log_error(P, _("latitude max < latitude min."));
+            proj_errno_set (P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
+            return false;
+        }
+        // handle antimeridian
+        delta_y = (ymax - ymin + 360.0) / side_pts;
+    } else {
+        delta_y = (ymax - ymin) / side_pts;
+    }
+
+
+    // build densified bounding box
+    // Note: must be a linear ring for antimeridian logic
+    for( int iii = 0; iii < side_pts; iii++ )
+    {
+        // xmin boundary
+        y_boundary_array[iii] = ymax - iii * delta_y;
+        x_boundary_array[iii] = xmin;
+        // ymin boundary
+        y_boundary_array[iii + side_pts] = ymin;
+        x_boundary_array[iii + side_pts] = xmin + iii * delta_x;
+        // xmax boundary
+        y_boundary_array[iii + side_pts * 2] = ymin + iii * delta_y;
+        x_boundary_array[iii + side_pts * 2] = xmax;
+        // ymax boundary
+        y_boundary_array[iii + side_pts * 3] = ymax;
+        x_boundary_array[iii + side_pts * 3] = xmax - iii * delta_x;
+    }
+    proj_trans_generic (
+        P,
+        direction,
+        &x_boundary_array[0], sizeof(double), boundary_len,
+        &y_boundary_array[0], sizeof(double), boundary_len,
+        nullptr, 0, 0,
+        nullptr, 0, 0
+    );
+
+    if (!degree_output) {
+        *out_xmin = simple_min(&x_boundary_array[0], boundary_len);
+        *out_xmax = simple_max(&x_boundary_array[0], boundary_len);
+        *out_ymin = simple_min(&y_boundary_array[0], boundary_len);
+        *out_ymax = simple_max(&y_boundary_array[0], boundary_len);
+    } else if (north_pole_in_bounds && output_lon_lat_order) {
+        *out_xmin = -180;
+        *out_ymin = simple_min(&y_boundary_array[0], boundary_len);
+        *out_xmax = 180;
+        *out_ymax = 90;
+    } else if (north_pole_in_bounds) {
+        *out_xmin = simple_min(&x_boundary_array[0], boundary_len);
+        *out_ymin = -180;
+        *out_xmax = 90;
+        *out_ymax = 180;
+    } else if (south_pole_in_bounds && output_lon_lat_order) {
+        *out_xmin = -180;
+        *out_ymin = -90;
+        *out_xmax = 180;
+        *out_ymax = simple_max(&y_boundary_array[0], boundary_len);
+    } else if (south_pole_in_bounds) {
+        *out_xmin = -90;
+        *out_ymin = -180;
+        *out_xmax = simple_max(&x_boundary_array[0], boundary_len);
+        *out_ymax = 180;
+    } else if (output_lon_lat_order) {
+        *out_xmin = antimeridian_min(&x_boundary_array[0], boundary_len);
+        *out_xmax = antimeridian_max(&x_boundary_array[0], boundary_len);
+        *out_ymin = simple_min(&y_boundary_array[0], boundary_len);
+        *out_ymax = simple_max(&y_boundary_array[0], boundary_len);
+    } else {
+        *out_xmin = simple_min(&x_boundary_array[0], boundary_len);
+        *out_xmax = simple_max(&x_boundary_array[0], boundary_len);
+        *out_ymin = antimeridian_min(&y_boundary_array[0], boundary_len);
+        *out_ymax = antimeridian_max(&y_boundary_array[0], boundary_len);
+    }
+    return true;
+}
+
+
 /*****************************************************************************/
 static void reproject_bbox(PJ* pjGeogToCrs,
                           double west_lon, double south_lat,
diff --git a/src/proj.h b/src/proj.h
index 1af1440f..b378cd00 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -599,7 +599,22 @@ size_t PROJ_DLL proj_trans_generic (
     double *z, size_t sz, size_t nz,
     double *t, size_t st, size_t nt
 );
-
+/*! @endcond */
+int PROJ_DLL proj_trans_bounds(
+    PJ_CONTEXT* context,
+    PJ *P,
+    PJ_DIRECTION direction,
+    const double xmin,
+    const double ymin,
+    const double xmax,
+    const double ymax,
+    double* out_xmin,
+    double* out_ymin,
+    double* out_xmax,
+    double* out_ymax,
+    int densify_pts
+);
+/*! @cond Doxygen_Suppress */
 
 /* Initializers */
 PJ_COORD PROJ_DLL proj_coord (double x, double y, double z, double t);
diff --git a/test/unit/test_c_api.cpp b/test/unit/test_c_api.cpp
index aceae656..a334ccff 100644
--- a/test/unit/test_c_api.cpp
+++ b/test/unit/test_c_api.cpp
@@ -5678,6 +5678,543 @@ TEST_F(CApi, proj_get_geoid_models_from_database) {
 
 // ---------------------------------------------------------------------------
 
+TEST_F(CApi, proj_trans_bounds_densify_0) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:4326",
+        "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 "
+        "+a=6370997 +b=6370997 +units=m +no_defs",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, P, PJ_FWD,
+        40, -120, 64, -80,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        0
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, -1684649.41338, 1);
+    EXPECT_NEAR(out_bottom, -350356.81377, 1);
+    EXPECT_NEAR(out_right, 1684649.41338, 1);
+    EXPECT_NEAR(out_top, 2234551.18559, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds_densify_100) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:4326",
+        "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 "
+        "+a=6370997 +b=6370997 +units=m +no_defs",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, P, PJ_FWD,
+        40, -120, 64, -80,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        100
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, -1684649.41338, 1);
+    EXPECT_NEAR(out_bottom, -555777.79210, 1);
+    EXPECT_NEAR(out_right, 1684649.41338, 1);
+    EXPECT_NEAR(out_top, 2234551.18559, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds_normalized) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:4326",
+        "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 "
+        "+a=6370997 +b=6370997 +units=m +no_defs",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto normalized_p = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper normal_keeper_P(normalized_p);
+    ASSERT_NE(normalized_p, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, normalized_p, PJ_FWD,
+        -120, 40, -80, 64,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        100
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, -1684649.41338, 1);
+    EXPECT_NEAR(out_bottom, -555777.79210, 1);
+    EXPECT_NEAR(out_right, 1684649.41338, 1);
+    EXPECT_NEAR(out_top, 2234551.18559, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds_antimeridian_xy) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:4167",
+        "EPSG:3851",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto normalized_p = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper normal_keeper_P(normalized_p);
+    ASSERT_NE(normalized_p, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, normalized_p, PJ_FWD,
+        160.6, -55.95, -171.2, -25.88,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, 1722483.900174921, 1);
+    EXPECT_NEAR(out_bottom, 5228058.6143420935, 1);
+    EXPECT_NEAR(out_right, 4624385.494808555, 1);
+    EXPECT_NEAR(out_top, 8692574.544944234, 1);
+    double out_left_inv;
+    double out_bottom_inv;
+    double out_right_inv;
+    double out_top_inv;
+    int success_inv = proj_trans_bounds(
+        m_ctxt, normalized_p, PJ_INV,
+        1722483.900174921, 5228058.6143420935, 4624385.494808555, 8692574.544944234,
+        &out_left_inv,
+        &out_bottom_inv,
+        &out_right_inv,
+        &out_top_inv,
+        21
+    );
+    EXPECT_TRUE(success_inv == 1);
+    EXPECT_NEAR(out_left_inv, 153.2799922, 1);
+    EXPECT_NEAR(out_bottom_inv, -56.7471249, 1);
+    EXPECT_NEAR(out_right_inv, -162.1813873, 1);
+    EXPECT_NEAR(out_top_inv, -24.6148194, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds_antimeridian) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:4167",
+        "EPSG:3851",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, P, PJ_FWD,
+        -55.95, 160.6, -25.88, -171.2,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, 5228058.6143420935, 1);
+    EXPECT_NEAR(out_bottom, 1722483.900174921, 1);
+    EXPECT_NEAR(out_right, 8692574.544944234, 1);
+    EXPECT_NEAR(out_top, 4624385.494808555, 1);
+    double out_left_inv;
+    double out_bottom_inv;
+    double out_right_inv;
+    double out_top_inv;
+    int success_inv = proj_trans_bounds(
+        m_ctxt, P, PJ_INV,
+        5228058.6143420935, 1722483.900174921, 8692574.544944234, 4624385.494808555,
+        &out_left_inv,
+        &out_bottom_inv,
+        &out_right_inv,
+        &out_top_inv,
+        21
+    );
+    EXPECT_TRUE(success_inv == 1);
+    EXPECT_NEAR(out_left_inv, -56.7471249, 1);
+    EXPECT_NEAR(out_bottom_inv, 153.2799922, 1);
+    EXPECT_NEAR(out_right_inv, -24.6148194, 1);
+    EXPECT_NEAR(out_top_inv, -162.1813873, 1);
+}
+
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds_beyond_global_bounds) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:6933",
+        "EPSG:4326",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto normalized_p = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper normal_keeper_P(normalized_p);
+    ASSERT_NE(normalized_p, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, normalized_p, PJ_FWD,
+        -17367531.3203125, -7314541.19921875, 17367531.3203125, 7314541.19921875,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, -180, 1);
+    EXPECT_NEAR(out_bottom, -85.0445994113099, 1);
+    EXPECT_NEAR(out_right, 180, 1);
+    EXPECT_NEAR(out_top, 85.0445994113099, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds_ignore_inf) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "OGC:CRS84",
+        "ESRI:102036",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, P, PJ_FWD,
+        -180.0, -90.0, 180.0, 0.0,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, 0, 1);
+    EXPECT_NEAR(out_bottom, -89178008, 1);
+    EXPECT_NEAR(out_right, 0, 1);
+    EXPECT_NEAR(out_top, 0, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds_ignore_inf_geographic) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "PROJCS[\"Interrupted_Goode_Homolosine\","
+        "GEOGCS[\"GCS_unnamed ellipse\",DATUM[\"D_unknown\","
+        "SPHEROID[\"Unknown\",6378137,298.257223563]],"
+        "PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]],"
+        "PROJECTION[\"Interrupted_Goode_Homolosine\"],"
+        "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],"
+        "AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]",
+        "OGC:CRS84",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, P, PJ_FWD,
+        -15028000.0, 7515000.0, -14975000.0, 7556000.0,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, -179.2133, 1);
+    EXPECT_NEAR(out_bottom, 70.9345, 1);
+    EXPECT_NEAR(out_right, -177.9054, 1);
+    EXPECT_NEAR(out_top, 71.4364, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds_noop_geographic) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:4284",
+        "EPSG:4284",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, P, PJ_FWD,
+        19.57, 35.14, -168.97, 81.91,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, 19.57, 1);
+    EXPECT_NEAR(out_bottom, 35.14, 1);
+    EXPECT_NEAR(out_right, -168.97, 1);
+    EXPECT_NEAR(out_top, 81.91, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds__north_pole_xy) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:32661",
+        "EPSG:4326",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto normalized_p = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper normal_keeper_P(normalized_p);
+    ASSERT_NE(normalized_p, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, normalized_p, PJ_FWD,
+        -1371213.7625429356, -1405880.71737131, 5371213.762542935, 5405880.71737131,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, -180.0, 1);
+    EXPECT_NEAR(out_bottom, 48.656, 1);
+    EXPECT_NEAR(out_right, 180.0, 1);
+    EXPECT_NEAR(out_top, 90.0, 1);
+    double out_left_inv;
+    double out_bottom_inv;
+    double out_right_inv;
+    double out_top_inv;
+    int success_inv = proj_trans_bounds(
+        m_ctxt, normalized_p, PJ_INV,
+        -180.0, 60.0, 180.0, 90.0,
+        &out_left_inv,
+        &out_bottom_inv,
+        &out_right_inv,
+        &out_top_inv,
+        21
+    );
+    EXPECT_TRUE(success_inv == 1);
+    EXPECT_NEAR(out_left_inv, -1371213.76, 1);
+    EXPECT_NEAR(out_bottom_inv, -1405880.72, 1);
+    EXPECT_NEAR(out_right_inv, 5371213.76, 1);
+    EXPECT_NEAR(out_top_inv, 5405880.72, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds__north_pole) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:32661",
+        "EPSG:4326",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, P, PJ_FWD,
+        -1405880.71737131, -1371213.7625429356, 5405880.71737131, 5371213.762542935,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, 48.656, 1);
+    EXPECT_NEAR(out_bottom, -180.0, 1);
+    EXPECT_NEAR(out_right, 90.0, 1);
+    EXPECT_NEAR(out_top, 180.0, 1);
+    double out_left_inv;
+    double out_bottom_inv;
+    double out_right_inv;
+    double out_top_inv;
+    int success_inv = proj_trans_bounds(
+        m_ctxt, P, PJ_INV,
+        60.0, -180.0, 90.0, 180.0,
+        &out_left_inv,
+        &out_bottom_inv,
+        &out_right_inv,
+        &out_top_inv,
+        21
+    );
+    EXPECT_TRUE(success_inv == 1);
+    EXPECT_NEAR(out_left_inv, -1405880.72, 1);
+    EXPECT_NEAR(out_bottom_inv, -1371213.76, 1);
+    EXPECT_NEAR(out_right_inv, 5405880.72, 1);
+    EXPECT_NEAR(out_top_inv, 5371213.76, 1);
+}
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds__south_pole_xy) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:32761",
+        "EPSG:4326",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    auto normalized_p = proj_normalize_for_visualization(m_ctxt, P);
+    ObjectKeeper normal_keeper_P(normalized_p);
+    ASSERT_NE(normalized_p, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, normalized_p, PJ_FWD,
+        -1371213.7625429356, -1405880.71737131, 5371213.762542935, 5405880.71737131,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, -180.0, 1);
+    EXPECT_NEAR(out_bottom, -90, 1);
+    EXPECT_NEAR(out_right, 180.0, 1);
+    EXPECT_NEAR(out_top, -48.656, 1);
+    double out_left_inv;
+    double out_bottom_inv;
+    double out_right_inv;
+    double out_top_inv;
+    int success_inv = proj_trans_bounds(
+        m_ctxt, normalized_p, PJ_INV,
+        -180.0, -90.0, 180.0, -60.0,
+        &out_left_inv,
+        &out_bottom_inv,
+        &out_right_inv,
+        &out_top_inv,
+        21
+    );
+    EXPECT_TRUE(success_inv == 1);
+    EXPECT_NEAR(out_left_inv, -1371213.76, 1);
+    EXPECT_NEAR(out_bottom_inv, -1405880.72, 1);
+    EXPECT_NEAR(out_right_inv, 5371213.76, 1);
+    EXPECT_NEAR(out_top_inv, 5405880.72, 1);
+}
+
+
+// ---------------------------------------------------------------------------
+
+TEST_F(CApi, proj_trans_bounds__south_pole) {
+    auto P = proj_create_crs_to_crs(
+        m_ctxt,
+        "EPSG:32761",
+        "EPSG:4326",
+        nullptr
+    );
+    ObjectKeeper keeper_P(P);
+    ASSERT_NE(P, nullptr);
+    double out_left;
+    double out_bottom;
+    double out_right;
+    double out_top;
+    int success = proj_trans_bounds(
+        m_ctxt, P, PJ_FWD,
+        -1405880.71737131, -1371213.7625429356, 5405880.71737131, 5371213.762542935,
+        &out_left,
+        &out_bottom,
+        &out_right,
+        &out_top,
+        21
+    );
+    EXPECT_TRUE(success == 1);
+    EXPECT_NEAR(out_left, -90.0, 1);
+    EXPECT_NEAR(out_bottom, -180.0, 1);
+    EXPECT_NEAR(out_right, -48.656, 1);
+    EXPECT_NEAR(out_top, 180.0, 1);
+    double out_left_inv;
+    double out_bottom_inv;
+    double out_right_inv;
+    double out_top_inv;
+    int success_inv = proj_trans_bounds(
+        m_ctxt, P, PJ_INV,
+        -90.0, -180.0, -60.0, 180.0,
+        &out_left_inv,
+        &out_bottom_inv,
+        &out_right_inv,
+        &out_top_inv,
+        21
+    );
+    EXPECT_TRUE(success_inv == 1);
+    EXPECT_NEAR(out_left_inv, -1405880.72, 1);
+    EXPECT_NEAR(out_bottom_inv, -1371213.76, 1);
+    EXPECT_NEAR(out_right_inv, 5405880.72, 1);
+    EXPECT_NEAR(out_top_inv, 5371213.76, 1);
+}
+
+// ---------------------------------------------------------------------------
+
 #if !defined(_WIN32)
 TEST_F(CApi, open_plenty_of_contexts) {
     // Test that we only consume 1 file handle for the connection to the