diff --git a/pyproj/_geod.pyi b/pyproj/_geod.pyi
index 412ba57a1..ed0196761 100644
--- a/pyproj/_geod.pyi
+++ b/pyproj/_geod.pyi
@@ -45,6 +45,19 @@ class Geod:
         include_initial: bool = False,
         include_terminus: bool = False,
     ) -> None: ...
+    def _fwd_npts(
+        self,
+        out_lons: Any,
+        out_lats: Any,
+        out_azis: Any,
+        lon1: float,
+        lat1: float,
+        az1: float,
+        npts: int,
+        radians: bool = False,
+        include_initial: bool = False,
+        include_terminus: bool = False,
+    ) -> None: ...
     def _line_length(self, lons: Any, lats: Any, radians: bool = False) -> float: ...
     def _polygon_area_perimeter(
         self, lons: Any, lats: Any, radians: bool = False
diff --git a/pyproj/_geod.pyx b/pyproj/_geod.pyx
index d99520f1c..7e02b2611 100644
--- a/pyproj/_geod.pyx
+++ b/pyproj/_geod.pyx
@@ -245,6 +245,66 @@ cdef class Geod:
                 if store_az:
                     azis_buff.data[iii] = pazi2
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def _fwd_npts(
+        self,
+        object out_lons,
+        object out_lats,
+        object out_azis,
+        double lon1,
+        double lat1,
+        double az1,
+        int npts,
+        bint radians=False,
+        bint include_initial=False,
+        bint include_terminus=False,
+    ):
+        """
+        given initial and terminus lat/lon, find npts intermediate points.
+        using given lons, lats buffers
+        """
+        cdef Py_ssize_t iii
+        cdef double del_s
+        cdef double pazi2
+        cdef double s12
+        cdef double plon2
+        cdef double plat2
+        cdef geod_geodesicline line
+
+        cdef PyBuffWriteManager lons_buff = PyBuffWriteManager(out_lons)
+        cdef PyBuffWriteManager lats_buff = PyBuffWriteManager(out_lats)
+
+        if lons_buff.len < npts or lats_buff.len < npts:
+            raise GeodError("lons or lats arrays are not long enough.")
+
+        cdef PyBuffWriteManager azis_buff
+        cdef bint store_az = out_azis is not None
+
+        if store_az:
+            azis_buff = PyBuffWriteManager(out_azis)
+            if azis_buff.len < npts:
+                raise GeodError("az array is not long enough.")
+
+        with nogil:
+            if radians:
+                lon1 *= _RAD2DG
+                lat1 *= _RAD2DG
+            # do fwd computation to set azimuths, distance.
+            geod_lineinit(&line, &self._geod_geodesic, lat1, lon1, azi1, 0u)
+            # distance increment.
+            del_s = line.s13 / (npts - include_initial - include_terminus + 1)
+            # loop over intermediate points, compute lat/lons.
+            for iii in range(0, npts):
+                s12 = (iii + (not include_initial)) * del_s
+                geod_position(&line, s12, &plat2, &plon2, &pazi2)
+                if radians:
+                    plat2 *= _DG2RAD
+                    plon2 *= _DG2RAD
+                lats_buff.data[iii] = plat2
+                lons_buff.data[iii] = plon2
+                if store_az:
+                    azis_buff.data[iii] = pazi2
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
diff --git a/pyproj/geod.py b/pyproj/geod.py
index 965d86092..6d3a19e97 100644
--- a/pyproj/geod.py
+++ b/pyproj/geod.py
@@ -417,6 +417,35 @@ def inv_npts(
             include_terminus=include_terminus,
         )
 
+    def fwd_npts(
+        self,
+        out_lons: Any,
+        out_lats: Any,
+        out_azis: Any,
+        lon1: float,
+        lat1: float,
+        az1: float,
+        npts: int,
+        radians: bool = False,
+        include_initial: bool = False,
+        include_terminus: bool = False,
+    ) -> None:
+        """
+        Similar to npts_inv(), but using a forward computation to calculate geod_line
+        """
+        super()._fwd_npts(
+            out_lons=out_lons,
+            out_lats=out_lats,
+            out_azis=out_azis,
+            lon1=lon1,
+            lat1=lat1,
+            az1=az1,
+            npts=npts,
+            radians=radians,
+            include_initial=include_initial,
+            include_terminus=include_terminus,
+        )
+
     def line_length(self, lons: Any, lats: Any, radians: bool = False) -> float:
         """
         .. versionadded:: 2.3.0