BUG: fix Geod.npts NaN handling (#1288)

docs/history.rst

@@ -7,6 +7,7 @@ Latest
 - ENH: Added allow_superseded kwargs to :class:`pyproj.transformer.TransformerGroup` (pull #1269)
 - ENH: Added :meth:`CRS.to_2d` to demote 3D CRS to 2D (issue #1266)
 - ENH: Added parameter `output_axis_rule` to :meth:`CRS.to_wkt` (pull #1287)
+- BUG: fix Geod.npts NaN handling (issue #1282)
 
 3.5.0
 ------

pyproj/_geod.pyi

@@ -63,7 +63,7 @@ class Geod:
         lon1: float,
         lat1: float,
         lon2_or_azi1: float,
-        lat2_or_nan: float,
+        lat2: float,
         npts: int,
         del_s: float,
         radians: bool,
@@ -74,6 +74,7 @@ class Geod:
         out_lats: Any,
         out_azis: Any,
         return_back_azimuth: bool,
+        is_fwd: bool,
     ) -> GeodIntermediateReturn: ...
     def _line_length(self, lons: Any, lats: Any, radians: bool = False) -> float: ...
     def _polygon_area_perimeter(

pyproj/_geod.pyx

@@ -365,7 +365,7 @@ cdef class Geod:
         double lon1,
         double lat1,
         double lon2_or_azi1,
-        double lat2_or_nan,
+        double lat2,
         int npts,
         double del_s,
         bint radians,
@@ -376,6 +376,7 @@ cdef class Geod:
         object out_lats,
         object out_azis,
         bint return_back_azimuth,
+        bint is_fwd,
     ) -> GeodIntermediateReturn:
         """
         .. versionadded:: 3.1.0
@@ -397,7 +398,6 @@ cdef class Geod:
             PyBuffWriteManager lons_buff
             PyBuffWriteManager lats_buff
             PyBuffWriteManager azis_buff
-            bint is_fwd = isnan(lat2_or_nan)
 
         if not is_fwd and (del_s == 0) == (npts == 0):
             raise GeodError("inv_intermediate: "
@@ -409,7 +409,7 @@ cdef class Geod:
                 lat1 *= _RAD2DG
                 lon2_or_azi1 *= _RAD2DG
                 if not is_fwd:
-                    lat2_or_nan *= _RAD2DG
+                    lat2 *= _RAD2DG
 
             if is_fwd:
                 # do fwd computation to set azimuths, distance.
@@ -418,7 +418,7 @@ cdef class Geod:
             else:
                 # do inverse computation to set azimuths, distance.
                 geod_inverseline(&line, &self._geod_geodesic, lat1, lon1,
-                                 lat2_or_nan, lon2_or_azi1, 0u)
+                                 lat2, lon2_or_azi1, 0u)
 
                 if npts == 0:
                     # calc the number of required points by the distance increment

pyproj/geod.py

@@ -524,7 +524,7 @@ def npts(
             lon1=lon1,
             lat1=lat1,
             lon2_or_azi1=lon2,
-            lat2_or_nan=lat2,
+            lat2=lat2,
             npts=npts,
             del_s=0,
             radians=radians,
@@ -535,6 +535,7 @@ def npts(
             out_lats=None,
             out_azis=None,
             return_back_azimuth=False,
+            is_fwd=False,
         )
         return list(zip(res.lons, res.lats))
 
@@ -688,7 +689,7 @@ def inv_intermediate(
             lon1=lon1,
             lat1=lat1,
             lon2_or_azi1=lon2,
-            lat2_or_nan=lat2,
+            lat2=lat2,
             npts=npts,
             del_s=del_s,
             radians=radians,
@@ -699,6 +700,7 @@ def inv_intermediate(
             out_lats=out_lats,
             out_azis=out_azis,
             return_back_azimuth=return_back_azimuth,
+            is_fwd=False,
         )
 
     def fwd_intermediate(
@@ -835,7 +837,7 @@ def fwd_intermediate(
             lon1=lon1,
             lat1=lat1,
             lon2_or_azi1=azi1,
-            lat2_or_nan=math.nan,
+            lat2=math.nan,
             npts=npts,
             del_s=del_s,
             radians=radians,
@@ -846,6 +848,7 @@ def fwd_intermediate(
             out_lats=out_lats,
             out_azis=out_azis,
             return_back_azimuth=return_back_azimuth,
+            is_fwd=True,
         )
 
     def line_length(self, lons: Any, lats: Any, radians: bool = False) -> float:

test/test_geod.py

@@ -5,7 +5,7 @@
 
 import numpy
 import pytest
-from numpy.testing import assert_almost_equal
+from numpy.testing import assert_almost_equal, assert_array_equal
 
 from pyproj import Geod
 from pyproj.geod import GeodIntermediateFlag, reverse_azimuth
@@ -443,6 +443,21 @@ def test_geodesic_npts(include_initial, include_terminus):
     assert_almost_equal(lonlats, expected_lonlats, decimal=3)
 
 
+@pytest.mark.parametrize(
+    "input_data",
+    [
+        [1.0, 2.0, 3.0, math.nan],
+        [1.0, 2.0, math.nan, 4.0],
+        [1.0, math.nan, 3.0, 4.0],
+        [math.nan, 2.0, 3.0, 4.0],
+        [math.nan, math.nan, math.nan, math.nan],
+    ],
+)
+def test_geodesic_npts__nan(input_data):
+    geod = Geod(ellps="WGS84")
+    assert_array_equal(geod.npts(*input_data, npts=1), [(math.nan, math.nan)])
+
+
 @pytest.mark.parametrize(
     "ellipsoid,expected_azi12,expected_az21,expected_dist",
     [("clrk66", -66.531, 75.654, 4164192.708), ("WGS84", -66.530, 75.654, 4164074.239)],