BUG: Handle axis order with poles in transform bounds (#966)

pyproj/_transformer.pyx

@@ -11,6 +11,7 @@ from collections import namedtuple
 from pyproj._compat cimport cstrencode
 from pyproj._crs cimport (
     _CRS,
+    Axis,
     Base,
     CoordinateOperation,
     _get_concatenated_operations,
@@ -519,6 +520,7 @@ cdef bint contains_north_pole(
     double bottom,
     double right,
     double top,
+    bint lon_lat_order,
 ) nogil:
     """
     Check if the original projected bounds contains
@@ -529,11 +531,14 @@ cdef bint contains_north_pole(
     # North Pole
     cdef double pole_y = 90
     cdef double pole_x = 0
-    cdef PJ_DIRECTION direction = PJ_IDENT
+    if not lon_lat_order:
+        pole_x = 90
+        pole_y = 0
+
+    cdef PJ_DIRECTION direction = PJ_INV
     if pj_direction == PJ_INV:
         direction = PJ_FWD
-    elif pj_direction == PJ_FWD:
-        direction = PJ_INV
+
     proj_trans_generic(
         projobj,
         direction,
@@ -554,6 +559,7 @@ cdef bint contains_south_pole(
     double bottom,
     double right,
     double top,
+    bint lon_lat_order,
 ) nogil:
     """
     Check if the original projected bounds contains
@@ -564,11 +570,13 @@ cdef bint contains_south_pole(
     # South Pole
     cdef double pole_y = -90
     cdef double pole_x = 0
-    cdef PJ_DIRECTION direction = PJ_IDENT
+    if not lon_lat_order:
+        pole_y = 0
+        pole_x = -90
+
+    cdef PJ_DIRECTION direction = PJ_INV
     if pj_direction == PJ_INV:
         direction = PJ_FWD
-    elif pj_direction == PJ_FWD:
-        direction = PJ_INV
 
     proj_trans_generic(
         projobj,
@@ -583,6 +591,36 @@ cdef bint contains_south_pole(
     return False
 
 
+cdef bint target_crs_lon_lat_order(
+    PJ_CONTEXT* transformer_ctx,
+    PJ* transformer_pj,
+    PJ_DIRECTION pj_direction,
+) except *:
+    cdef PJ* target_crs = NULL
+    cdef PJ* coord_system_pj = NULL
+    cdef Axis first_axis = None
+    if pj_direction == PJ_FWD:
+        target_crs = proj_get_target_crs(transformer_ctx, transformer_pj)
+    elif pj_direction == PJ_INV:
+        target_crs = proj_get_source_crs(transformer_ctx, transformer_pj)
+    if target_crs == NULL:
+        raise ProjError("Unable to retrieve target CRS")
+    try:
+        coord_system_pj = proj_crs_get_coordinate_system(
+            transformer_ctx,
+            target_crs,
+        )
+        if coord_system_pj == NULL:
+            raise ProjError("Unable to get target CRS coordinate system.")
+        first_axis = Axis.create(transformer_ctx, coord_system_pj, 0)
+    finally:
+        proj_destroy(target_crs)
+        if coord_system_pj != NULL:
+            proj_destroy(coord_system_pj)
+    ProjError.clear()
+    return first_axis.abbrev.lower() == "lon"
+
+
 cdef class _Transformer(Base):
     def __cinit__(self):
         self._area_of_use = None
@@ -1004,20 +1042,25 @@ cdef class _Transformer(Base):
         bint errcheck,
         object direction,
     ):
-        if self.id == "noop":
+        tmp_pj_direction = _PJ_DIRECTION_MAP[TransformDirection.create(direction)]
+        cdef PJ_DIRECTION pj_direction = <PJ_DIRECTION>tmp_pj_direction
+
+        if self.id == "noop" or pj_direction == PJ_IDENT:
             return (left, bottom, right, top)
 
         if densify_pts < 0:
             raise ProjError("densify_pts must be positive")
 
-        tmp_pj_direction = _PJ_DIRECTION_MAP[TransformDirection.create(direction)]
-        cdef PJ_DIRECTION pj_direction = <PJ_DIRECTION>tmp_pj_direction
-
         cdef bint degree_output = proj_degree_output(self.projobj, pj_direction)
         cdef bint degree_input = proj_degree_input(self.projobj, pj_direction)
         if degree_output and densify_pts < 2:
             raise ProjError("densify_pts must be 2+ for degree output")
 
+        cdef bint output_lon_lat_order = False
+        if degree_output:
+            output_lon_lat_order = target_crs_lon_lat_order(
+                self.context, self.projobj, pj_direction,
+            )
         cdef int side_pts = densify_pts + 1  # add one because we are densifying
         cdef Py_ssize_t boundary_len = side_pts * 4
         cdef PySimpleArray x_boundary_array = PySimpleArray(boundary_len)
@@ -1036,6 +1079,7 @@ cdef class _Transformer(Base):
                     bottom,
                     right,
                     top,
+                    output_lon_lat_order,
                 )
                 south_pole_in_bounds = contains_south_pole(
                     self.projobj,
@@ -1044,6 +1088,7 @@ cdef class _Transformer(Base):
                     bottom,
                     right,
                     top,
+                    output_lon_lat_order,
                 )
 
             # degrees to radians
@@ -1111,19 +1156,34 @@ cdef class _Transformer(Base):
             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:
+                if north_pole_in_bounds and output_lon_lat_order:
+                    left = -180
                     bottom = simple_min(y_boundary_array.data, y_boundary_array.len)
+                    right = 180
                     top = 90
-                elif south_pole_in_bounds:
+                elif north_pole_in_bounds:
+                    left = simple_min(x_boundary_array.data, x_boundary_array.len)
+                    bottom = -180
+                    right = 90
+                    top = 180
+                elif south_pole_in_bounds and output_lon_lat_order:
+                    left = -180
                     bottom = -90
+                    right = 180
                     top = simple_max(y_boundary_array.data, y_boundary_array.len)
-                left = -180
-                right = 180
-            elif degree_output:
+                else:
+                    left = -90
+                    right = simple_max(x_boundary_array.data, x_boundary_array.len)
+                    bottom = -180
+                    top = 180
+            elif degree_output and output_lon_lat_order:
                 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
-                # to be on the y axis. It shouldn't cause troubles if it is latitude.
+                bottom = simple_min(y_boundary_array.data, y_boundary_array.len)
+                top = simple_max(y_boundary_array.data, y_boundary_array.len)
+            elif degree_output:
+                left = simple_min(x_boundary_array.data, x_boundary_array.len)
+                right = simple_max(x_boundary_array.data, x_boundary_array.len)
                 bottom = antimeridian_min(y_boundary_array.data, y_boundary_array.len)
                 top = antimeridian_max(y_boundary_array.data, y_boundary_array.len)
             else:

pyproj/transformer.py

@@ -895,21 +895,21 @@ 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
         ----------
         left: float
-            Left bounding coordinate in source CRS.
+            Minimum bounding coordinate of the first axis in source CRS
+            (or the target CRS if using the reverse direction).
         bottom: float
-            Bottom bounding coordinate in source CRS.
+            Minimum bounding coordinate of the second axis in source CRS.
+            (or the target CRS if using the reverse direction).
         right: float
-            Right bounding coordinate in source CRS.
+            Maximum bounding coordinate of the first axis in source CRS.
+            (or the target CRS if using the reverse direction).
         top: float
-            Top bounding coordinate in source CRS.
+            Maximum bounding coordinate of the second axis in source CRS.
+            (or the target CRS if using the reverse direction).
         densify_points: uint, default=21
             Number of points to add to each edge to account for nonlinear edges
             produced by the transform process. Large numbers will produce worse

test/test_transformer.py

@@ -1375,6 +1375,28 @@ def test_transform_bounds__noop_geographic():
 
 
 def test_transform_bounds__north_pole():
+    crs = CRS("EPSG:32661")
+    transformer = Transformer.from_crs(crs, "EPSG:4326")
+    minx, miny, maxx, maxy = crs.area_of_use.bounds
+    bounds = transformer.transform_bounds(miny, minx, maxy, maxx, direction="INVERSE")
+    assert_almost_equal(
+        bounds,
+        (
+            -1405880.72,
+            -1371213.76,
+            5405880.72,
+            5371213.76,
+        ),
+        decimal=0,
+    )
+    assert_almost_equal(
+        transformer.transform_bounds(*bounds),
+        (48.656, -180.0, 90.0, 180.0),
+        decimal=0,
+    )
+
+
+def test_transform_bounds__north_pole__xy():
     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")
@@ -1391,6 +1413,28 @@ def test_transform_bounds__north_pole():
 
 
 def test_transform_bounds__south_pole():
+    crs = CRS("EPSG:32761")
+    transformer = Transformer.from_crs(crs, "EPSG:4326")
+    minx, miny, maxx, maxy = crs.area_of_use.bounds
+    bounds = transformer.transform_bounds(miny, minx, maxy, maxx, direction="INVERSE")
+    assert_almost_equal(
+        bounds,
+        (
+            -1405880.72,
+            -1371213.76,
+            5405880.72,
+            5371213.76,
+        ),
+        decimal=0,
+    )
+    assert_almost_equal(
+        transformer.transform_bounds(*bounds),
+        (-90, -180.0, -48.656, 180.0),
+        decimal=0,
+    )
+
+
+def test_transform_bounds__south_pole__xy():
     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")