diff --git a/docs/history.rst b/docs/history.rst index 5d507ea2a..95eff55ca 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -5,6 +5,7 @@ Latest ------- - BUG: Prepend "Derived" to CRS type name if CRS is derived (issue #932) - BUG: Improved handling of inf values in :attr:`pyproj.transformer.Transformer.transform_bounds` (pull #961) +- ENH: Add support for transforming bounds at the poles in :attr:`pyproj.transformer.Transformer.transform_bounds` (pull #962) 3.2.1 ------ diff --git a/pyproj/_transformer.pyx b/pyproj/_transformer.pyx index 5f4fafc41..91c68ef25 100644 --- a/pyproj/_transformer.pyx +++ b/pyproj/_transformer.pyx @@ -512,6 +512,77 @@ cdef double antimeridian_max(double* data, Py_ssize_t arr_len) nogil: return max_value +cdef bint contains_north_pole( + PJ* projobj, + PJ_DIRECTION pj_direction, + double left, + double bottom, + double right, + double top, +) nogil: + """ + Check if the original projected bounds contains + the north pole. + + This assumes that the destination CRS is geographic. + """ + # North Pole + cdef double pole_y = 90 + cdef double pole_x = 0 + cdef PJ_DIRECTION direction = PJ_IDENT + if pj_direction == PJ_INV: + direction = PJ_FWD + elif pj_direction == PJ_FWD: + direction = PJ_INV + proj_trans_generic( + projobj, + direction, + &pole_x, _DOUBLESIZE, 1, + &pole_y, _DOUBLESIZE, 1, + NULL, _DOUBLESIZE, 0, + NULL, _DOUBLESIZE, 0, + ) + if left < pole_x < right and top > pole_y > bottom: + return True + return False + + +cdef bint contains_south_pole( + PJ* projobj, + PJ_DIRECTION pj_direction, + double left, + double bottom, + double right, + double top, +) nogil: + """ + Check if the original projected bounds contains + the south pole. + + This assumes that the destination CRS is geographic. + """ + # South Pole + cdef double pole_y = -90 + cdef double pole_x = 0 + cdef PJ_DIRECTION direction = PJ_IDENT + if pj_direction == PJ_INV: + direction = PJ_FWD + elif pj_direction == PJ_FWD: + direction = PJ_INV + + proj_trans_generic( + projobj, + direction, + &pole_x, _DOUBLESIZE, 1, + &pole_y, _DOUBLESIZE, 1, + NULL, _DOUBLESIZE, 0, + NULL, _DOUBLESIZE, 0, + ) + if left < pole_x < right and top > pole_y > bottom: + return True + return False + + cdef class _Transformer(Base): def __cinit__(self): self._area_of_use = None @@ -954,8 +1025,27 @@ cdef class _Transformer(Base): cdef double delta_x = 0 cdef double delta_y = 0 cdef int iii = 0 - + cdef bint north_pole_in_bounds = False + cdef bint south_pole_in_bounds = False with nogil: + if degree_output: + north_pole_in_bounds = contains_north_pole( + self.projobj, + pj_direction, + left, + bottom, + right, + top, + ) + south_pole_in_bounds = contains_south_pole( + self.projobj, + pj_direction, + left, + bottom, + right, + top, + ) + # degrees to radians if not radians and proj_angular_input(self.projobj, pj_direction): left *= _DG2RAD @@ -1018,7 +1108,22 @@ cdef class _Transformer(Base): if ProjError.internal_proj_error is not None: raise ProjError("itransform error") - if degree_output: + if degree_output and (north_pole_in_bounds or south_pole_in_bounds): + # only works with lon/lat axis order + # need a way to test axis order to support both + if north_pole_in_bounds: + # counter intuitive, but the max actually represents the min + # https://github.com/OSGeo/PROJ/issues/2779#issuecomment-926253439 + bottom = simple_max(y_boundary_array.data, y_boundary_array.len) + top = 90 + elif south_pole_in_bounds: + bottom = -90 + # counter intuitive, but the min actually represents the max + # https://github.com/OSGeo/PROJ/issues/2779#issuecomment-926253439 + top = simple_min(y_boundary_array.data, y_boundary_array.len) + left = -180 + right = 180 + elif degree_output: left = antimeridian_min(x_boundary_array.data, x_boundary_array.len) right = antimeridian_max(x_boundary_array.data, x_boundary_array.len) # depending on the axis order, longitude has the potential diff --git a/pyproj/transformer.py b/pyproj/transformer.py index 3a029492b..57e2c8d96 100644 --- a/pyproj/transformer.py +++ b/pyproj/transformer.py @@ -895,6 +895,10 @@ def bounding_polygon(left, bottom, right, top): ) return shapely.geometry.box(left, bottom, right, top) + When projecting from polar projections to geographic, + lon, lat output order is required. The 'always_xy' flag + in the 'from_crs' constructor can help ensure that is the case. + Parameters ---------- diff --git a/test/test_transformer.py b/test/test_transformer.py index b9f27914c..2c355612c 100644 --- a/test/test_transformer.py +++ b/test/test_transformer.py @@ -1374,6 +1374,39 @@ def test_transform_bounds__noop_geographic(): ) +def test_transform_bounds__north_pole(): + crs = CRS("EPSG:32661") + transformer = Transformer.from_crs(crs, "EPSG:4326", always_xy=True) + bounds = transformer.transform_bounds(*crs.area_of_use.bounds, direction="INVERSE") + assert_almost_equal( + bounds, + (-1371213.76, -1405880.72, 5371213.76, 5405880.72), + decimal=0, + ) + print(transformer.transform_bounds(*bounds)) + assert_almost_equal( + transformer.transform_bounds(*bounds), + (-180.0, 60.29, 180.0, 90.0), + decimal=0, + ) + + +def test_transform_bounds__south_pole(): + crs = CRS("EPSG:32761") + transformer = Transformer.from_crs(crs, "EPSG:4326", always_xy=True) + bounds = transformer.transform_bounds(*crs.area_of_use.bounds, direction="INVERSE") + assert_almost_equal( + bounds, + (-1371213.76, -1405880.72, 5371213.76, 5405880.72), + decimal=0, + ) + assert_almost_equal( + transformer.transform_bounds(*bounds), + (-180.0, -90.0, 180.0, -60.29), + decimal=0, + ) + + @pytest.mark.parametrize("inplace", [True, False]) def test_transform__fortran_order(inplace): lons, lats = np.arange(-180, 180, 20), np.arange(-90, 90, 10)