diff --git a/docs/advanced_examples.rst b/docs/advanced_examples.rst index b4f3912ab..e4ac750a9 100644 --- a/docs/advanced_examples.rst +++ b/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 -------------- diff --git a/docs/history.rst b/docs/history.rst index 3152d2609..9218d0ded 100644 --- a/docs/history.rst +++ b/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 ----- diff --git a/pyproj/_transformer.pyx b/pyproj/_transformer.pyx index 8155d7f3f..d4157026c 100644 --- a/pyproj/_transformer.pyx +++ b/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 = 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 = 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): diff --git a/pyproj/transformer.py b/pyproj/transformer.py index 57527172e..d7fb5fed1 100644 --- a/pyproj/transformer.py +++ b/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, diff --git a/test/test_transformer.py b/test/test_transformer.py index 4da2540f2..26f4c9716 100644 --- a/test/test_transformer.py +++ b/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), + )