BUG: Improved handling of inf values in pyproj.transformer.Transformer.transform_bounds

docs/history.rst

@@ -4,6 +4,7 @@ Change Log
 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)
 
 3.2.1
 ------

pyproj/_transformer.pyx

@@ -340,6 +340,20 @@ cdef double simple_max(double* data, Py_ssize_t arr_len) nogil:
     return max_value
 
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef int _find_previous_index(int iii, double* data, int arr_len) nogil:
+    # find index of nearest valid previous value if exists
+    cdef int prev_iii = iii - 1
+    if prev_iii == -1:  # handle wraparound
+        prev_iii = arr_len - 1
+    while data[prev_iii] == HUGE_VAL and prev_iii != iii:
+        prev_iii -= 1
+        if prev_iii == -1:  # handle wraparound
+            prev_iii = arr_len - 1
+    return prev_iii
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef double antimeridian_min(double* data, Py_ssize_t arr_len) nogil:
@@ -410,19 +424,19 @@ cdef double antimeridian_min(double* data, Py_ssize_t arr_len) nogil:
     cdef bint positive_meridian = False
 
     for iii in range(0, arr_len):
-        prev_iii = iii - 1
-        if prev_iii == -1:
-            prev_iii = arr_len - 1
+        if data[iii] == HUGE_VAL:
+            continue
+        prev_iii = _find_previous_index(iii, data, arr_len)
         # check if crossed meridian
         delta = data[prev_iii] - data[iii]
         # 180 -> -180
-        if delta >= 200:
+        if delta >= 200 and delta != HUGE_VAL:
             if crossed_meridian_count == 0:
                 positive_min = min_value
             crossed_meridian_count += 1
             positive_meridian = False
         # -180 -> 180
-        elif delta <= -200:
+        elif delta <= -200 and delta != HUGE_VAL:
             if crossed_meridian_count == 0:
                 positive_min = data[iii]
             crossed_meridian_count += 1
@@ -433,6 +447,7 @@ cdef double antimeridian_min(double* data, Py_ssize_t arr_len) nogil:
         # track genral min value
         if data[iii] < min_value:
             min_value = data[iii]
+
     if crossed_meridian_count == 2:
         return positive_min
     elif crossed_meridian_count == 4:
@@ -463,19 +478,19 @@ cdef double antimeridian_max(double* data, Py_ssize_t arr_len) nogil:
     cdef bint negative_meridian = False
     cdef int crossed_meridian_count = 0
     for iii in range(0, arr_len):
-        prev_iii = iii - 1
-        if prev_iii == -1:
-            prev_iii = arr_len - 1
+        if data[iii] == HUGE_VAL:
+            continue
+        prev_iii = _find_previous_index(iii, data, arr_len)
         # check if crossed meridian
         delta = data[prev_iii] - data[iii]
         # 180 -> -180
-        if delta >= 200:
+        if delta >= 200 and delta != HUGE_VAL:
             if crossed_meridian_count == 0:
                 negative_max = data[iii]
             crossed_meridian_count += 1
             negative_meridian = True
         # -180 -> 180
-        elif delta <= -200:
+        elif delta <= -200 and delta != HUGE_VAL:
             if crossed_meridian_count == 0:
                 negative_max = max_value
             negative_meridian = False

test/test_transformer.py

@@ -1345,6 +1345,26 @@ def test_transform_bounds__ignore_inf(input_crs, expected_bounds):
     )
 
 
+def test_transform_bounds__ignore_inf_geographic():
+    crs_wkt = (
+        'PROJCS["Interrupted_Goode_Homolosine",'
+        'GEOGCS["GCS_unnamed ellipse",DATUM["D_unknown",'
+        'SPHEROID["Unknown",6378137,298.257223563]],'
+        'PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],'
+        'PROJECTION["Interrupted_Goode_Homolosine"],'
+        'UNIT["metre",1,AUTHORITY["EPSG","9001"]],'
+        'AXIS["Easting",EAST],AXIS["Northing",NORTH]]'
+    )
+    transformer = Transformer.from_crs(crs_wkt, "EPSG:4326", always_xy=True)
+    assert_almost_equal(
+        transformer.transform_bounds(
+            left=-15028000.0, bottom=7515000.0, right=-14975000.0, top=7556000.0
+        ),
+        (-179.2133, 70.9345, -177.9054, 71.4364),
+        decimal=0,
+    )
+
+
 def test_transform_bounds__noop_geographic():
     crs = CRS("Pulkovo 1942")
     transformer = Transformer.from_crs(crs.geodetic_crs, crs, always_xy=True)