ENH: Add Transformer.transform_bounds

docs/advanced_examples.rst

@@ -161,6 +161,22 @@ In PROJ 6+ you neeed to explictly change your CRS to 3D if you have
     >>> transformer_3d.transform(8.37909, 47.01987, 1000)
     (2671499.8913080636, 1208075.1135782297, 951.4265527743846)
 
+Projected CRS Bounds
+----------------------
+
+.. versionadded:: 3.1
+
+The boundary of the CRS is given in geographic coordinates.
+This is the recommended method for calculating the projected bounds.
+
+.. code-block:: python
+
+    >>> from pyproj import CRS, Transformer
+    >>> crs = CRS("EPSG:3857")
+    >>> transformer = Transformer.from_crs(crs.geodetic_crs, crs, always_xy=True)
+    >>> transformer.transform_bounds(*crs.area_of_use.bounds)
+    (-20037508.342789244, -20048966.104014594, 20037508.342789244, 20048966.104014594)
+
 
 Multithreading
 --------------

docs/history.rst

@@ -11,6 +11,7 @@ Change Log
 * ENH: Added authority, accuracy, and allow_ballpark kwargs to :meth:`pyproj.transformer.Transformer.from_crs` (issue #754)
 * ENH: Added support for "AUTH:CODE" input to :meth:`pyproj.transformer.Transformer.from_pipeline` (issue #755)
 * ENH: Added :meth:`pyproj.crs.CRS.to_3d` (pull #808)
+* ENH: Added :meth:`pyproj.transformer.Transformer.transform_bounds` (issue #809)
 
 3.0.1
 -----

pyproj/_transformer.pyx

@@ -2,6 +2,7 @@ include "base.pxi"
 
 cimport cython
 from cpython cimport array
+from cpython.mem cimport PyMem_Free, PyMem_Malloc
 
 import copy
 import re
@@ -302,6 +303,48 @@ cdef PJ* proj_create_crs_to_crs(
     return transform
 
 
+cdef class PySimpleArray:
+    cdef double* data
+    cdef public Py_ssize_t len
+
+    def __cinit__(self):
+        self.data = NULL
+
+    def __init__(self, Py_ssize_t arr_len):
+        self.len = arr_len
+        self.data = <double*> PyMem_Malloc(arr_len * sizeof(double))
+        if self.data == NULL:
+            raise MemoryError("error creating array for pyproj")
+
+    def __dealloc__(self):
+        PyMem_Free(self.data)
+        self.data = NULL
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef double simple_min(double* data, Py_ssize_t arr_len) nogil:
+    cdef int iii = 0
+    cdef double min_value = data[0]
+    for iii in range(1, arr_len):
+        if data[iii] < min_value:
+            min_value = data[iii]
+    return min_value
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef double simple_max(double* data, Py_ssize_t arr_len) nogil:
+    cdef int iii = 0
+    cdef double max_value = data[0]
+    for iii in range(1, arr_len):
+        if (data[iii] > max_value or
+            (max_value == HUGE_VAL and data[iii] != HUGE_VAL)
+        ):
+            max_value = data[iii]
+    return max_value
+
+
 cdef class _Transformer(Base):
     def __cinit__(self):
         self._area_of_use = None
@@ -743,6 +786,113 @@ cdef class _Transformer(Base):
 
         ProjError.clear()
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def _transform_bounds(
+        self,
+        double left,
+        double bottom,
+        double right,
+        double top,
+        int densify_pts,
+        object direction,
+        bint radians,
+        bint errcheck,
+    ):
+        if (
+            self.projections_exact_same
+            or (self.projections_equivalent and self.skip_equivalent)
+        ):
+            return
+
+        if densify_pts < 0:
+            raise ProjError("densify_pts must be positive")
+        cdef int side_pts = densify_pts + 1  # add one because we are densifying
+        tmp_pj_direction = _PJ_DIRECTION_MAP[TransformDirection.create(direction)]
+        cdef PJ_DIRECTION pj_direction = <PJ_DIRECTION>tmp_pj_direction
+        cdef Py_ssize_t boundary_len = side_pts * 4
+        cdef PySimpleArray x_boundary_array = PySimpleArray(boundary_len)
+        cdef PySimpleArray y_boundary_array = PySimpleArray(boundary_len)
+        cdef double delta_x = 0
+        cdef double delta_y = 0
+        cdef int iii = 0
+
+        with nogil:
+            # degrees to radians
+            if not radians and proj_angular_input(self.projobj, pj_direction):
+                left *= _DG2RAD
+                bottom *= _DG2RAD
+                right *= _DG2RAD
+                top *= _DG2RAD
+            # radians to degrees
+            elif radians and proj_degree_input(self.projobj, pj_direction):
+                left *= _RAD2DG
+                bottom *= _RAD2DG
+                right *= _RAD2DG
+                top *= _RAD2DG
+
+            if proj_degree_input(self.projobj, pj_direction) and right < left:
+                # handle antimeridian
+                delta_x = (right - left + 360.0) / side_pts
+            else:
+                delta_x = (right - left) / side_pts
+            delta_y = (top - bottom) / side_pts
+
+            # build densified bounding box
+            for iii in range(side_pts):
+                # left boundary
+                y_boundary_array.data[iii] = top - iii * delta_y
+                x_boundary_array.data[iii] = left
+                # bottom boundary
+                y_boundary_array.data[iii + side_pts] = bottom
+                x_boundary_array.data[iii + side_pts] = left + iii * delta_x
+                # right boundary
+                y_boundary_array.data[iii + side_pts * 2] = bottom + iii * delta_y
+                x_boundary_array.data[iii + side_pts * 2] = right
+                # top boundary
+                y_boundary_array.data[iii + side_pts * 3] = top
+                x_boundary_array.data[iii + side_pts * 3] = right - iii * delta_x
+
+            proj_errno_reset(self.projobj)
+            proj_trans_generic(
+                self.projobj,
+                pj_direction,
+                x_boundary_array.data, _DOUBLESIZE, x_boundary_array.len,
+                y_boundary_array.data, _DOUBLESIZE, y_boundary_array.len,
+                NULL, 0, 0,
+                NULL, 0, 0,
+            )
+            errno = proj_errno(self.projobj)
+            if errcheck and errno:
+                with gil:
+                    raise ProjError(
+                        f"itransform error: {pyproj_errno_string(self.context, errno)}"
+                    )
+            elif errcheck:
+                with gil:
+                    if ProjError.internal_proj_error is not None:
+                        raise ProjError("itransform error")
+
+            left = simple_min(x_boundary_array.data, x_boundary_array.len)
+            right = simple_max(x_boundary_array.data, x_boundary_array.len)
+            bottom = simple_min(y_boundary_array.data, y_boundary_array.len)
+            top = simple_max(y_boundary_array.data, y_boundary_array.len)
+            # radians to degrees
+            if not radians and proj_angular_output(self.projobj, pj_direction):
+                left *= _RAD2DG
+                bottom *= _RAD2DG
+                right *= _RAD2DG
+                top *= _RAD2DG
+            # degrees to radians
+            elif radians and proj_degree_output(self.projobj, pj_direction):
+                left *= _DG2RAD
+                bottom *= _DG2RAD
+                right *= _DG2RAD
+                top *= _DG2RAD
+
+        ProjError.clear()
+        return left, bottom, right, top
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _get_factors(self, longitude, latitude, bint radians, bint errcheck):

pyproj/transformer.py

@@ -806,6 +806,67 @@ def itransform(
             for pt in zip(*([iter(buff)] * stride)):
                 yield pt
 
+    def transform_bounds(
+        self,
+        left: float,
+        bottom: float,
+        right: float,
+        top: float,
+        densify_pts: int = 21,
+        radians: bool = False,
+        errcheck: bool = False,
+        direction: Union[TransformDirection, str] = TransformDirection.FORWARD,
+    ) -> Tuple[float, float, float, float]:
+        """
+        .. versionadded:: 3.1.0
+
+        Transform boundary densifying the edges to account for nonlinear
+        transformations along these edges and extracting the outermost bounds.
+
+        .. note:: Does not account for the antimeridian on destination bounds.
+
+        Parameters
+        ----------
+        left: float
+            Left bounding coordinate in source CRS.
+        bottom: float
+            Bottom bounding coordinate in source CRS.
+        right: float
+            Right bounding coordinate in source CRS.
+        top: float
+            Top bounding coordinate in source CRS.
+        densify_points: uint, optional
+            Number of points to add to each edge to account for nonlinear edges
+            produced by the transform process. Large numbers will produce worse
+            performance. Default: 21 (gdal default).
+        radians: boolean, optional
+            If True, will expect input data to be in radians and will return radians
+            if the projection is geographic. Default is False (degrees). Ignored for
+            pipeline transformations.
+        errcheck: boolean, optional
+            If True an exception is raised if the transformation is invalid.
+            By default errcheck=False and an invalid transformation
+            returns ``inf`` and no exception is raised.
+        direction: pyproj.enums.TransformDirection, optional
+            The direction of the transform.
+            Default is :attr:`pyproj.enums.TransformDirection.FORWARD`.
+
+        Returns
+        -------
+        left, bottom, right, top: float
+            Outermost coordinates in target coordinate reference system.
+        """
+        return self._transformer._transform_bounds(
+            left=left,
+            bottom=bottom,
+            right=right,
+            top=top,
+            direction=direction,
+            radians=radians,
+            errcheck=errcheck,
+            densify_pts=densify_pts,
+        )
+
     def to_proj4(
         self,
         version: Union[ProjVersion, str] = ProjVersion.PROJ_5,

test/test_transformer.py

@@ -11,7 +11,7 @@
 from numpy.testing import assert_almost_equal, assert_array_equal
 
 import pyproj
-from pyproj import Proj, Transformer, itransform, transform
+from pyproj import CRS, Proj, Transformer, itransform, transform
 from pyproj.datadir import append_data_dir
 from pyproj.enums import TransformDirection
 from pyproj.exceptions import ProjError
@@ -1220,3 +1220,66 @@ def test_transformer_from_pipeline__wkt_json(method_name):
         ).description
         == "RGF93 to WGS 84 (1)"
     )
+
+
+@pytest.mark.parametrize(
+    "density,expected",
+    [
+        (0, (-1684649.41338, -350356.81377, 1684649.41338, 2234551.18559)),
+        (100, (-1684649.41338, -555777.79210, 1684649.41338, 2234551.18559)),
+    ],
+)
+def test_transform_bounds_densify(density, expected):
+    transformer = Transformer.from_crs(
+        "EPSG:4326",
+        "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 "
+        "+a=6370997 +b=6370997 +units=m +no_defs",
+        always_xy=True,
+    )
+    assert np.allclose(
+        transformer.transform_bounds(-120, 40, -80, 64, densify_pts=density),
+        expected,
+    )
+
+
+def test_transform_bounds_densify_out_of_bounds():
+    transformer = Transformer.from_crs(
+        "EPSG:4326",
+        "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 "
+        "+a=6370997 +b=6370997 +units=m +no_defs",
+        always_xy=True,
+    )
+    with pytest.raises(ProjError):
+        transformer.transform_bounds(-120, 40, -80, 64, densify_pts=-1)
+
+
+def test_transform_bounds_radians():
+    transformer = Transformer.from_crs(
+        "EPSG:4326",
+        "+proj=geocent +ellps=WGS84 +datum=WGS84",
+        always_xy=True,
+    )
+    assert_almost_equal(
+        transformer.transform_bounds(
+            -2704026.010,
+            -4253051.810,
+            -2704025.010,
+            -4253050.810,
+            radians=True,
+            direction="INVERSE",
+        ),
+        (-2.1371136, 0.0, -2.1371133, 0.0),
+    )
+
+
+def test_transform_bounds__antimeridian():
+    crs = CRS("EPSG:3851")
+    transformer = Transformer.from_crs(
+        crs.geodetic_crs,
+        crs,
+        always_xy=True,
+    )
+    assert_almost_equal(
+        transformer.transform_bounds(*crs.area_of_use.bounds),
+        (1722483.900174921, 5228058.6143420935, 4624385.494808555, 8692574.544944234),
+    )