ENH: Add support for transforming bounds at the poles

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
 ------

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,20 @@ 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
+                    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
+                    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

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
         ----------

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)