diff --git a/pyproj/_geod.pyi b/pyproj/_geod.pyi
index 214aea51d..f34c2d7b1 100644
--- a/pyproj/_geod.pyi
+++ b/pyproj/_geod.pyi
@@ -31,6 +31,20 @@ class Geod:
         initial_idx: int = 1,
         terminus_idx: int = 1,
     ) -> Tuple[Tuple[float], Tuple[float]]: ...
+    def _inv_intermediate_by_npts(
+        self,
+        out_lons: Any,
+        out_lats: Any,
+        out_azis: Any,
+        lon1: float,
+        lat1: float,
+        lon2: float,
+        lat2: float,
+        npts: int,
+        radians: bool = False,
+        initial_idx: int = 1,
+        terminus_idx: int = 1,
+    ) -> int: ...
     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 19b889b81..33e8420e3 100644
--- a/pyproj/_geod.pyx
+++ b/pyproj/_geod.pyx
@@ -168,6 +168,38 @@ cdef class Geod:
         given initial and terminus lat/lon, find npts intermediate points.
         returns lons, lats buffers
         """
+        cdef array.array array_template = array.array("d", [])
+        cdef array.array out_lons = array.clone(array_template, npts, zero=False)
+        cdef array.array out_lats = array.clone(array_template, npts, zero=False)
+
+        self._inv_intermediate_by_npts(
+            out_lons=out_lons, out_lats=out_lats, out_azis=None,
+            lon1=lon1, lat1=lat1, lon2=lon2, lat2=lat2,
+            npts=npts, initial_idx=initial_idx, terminus_idx=terminus_idx,
+            radians=radians)
+
+        return out_lons, out_lats
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def _inv_intermediate_by_npts(
+        self,
+        object out_lons,
+        object out_lats,
+        object out_azis,
+        double lon1,
+        double lat1,
+        double lon2,
+        double lat2,
+        int npts,
+        bint radians=False,
+        int initial_idx=1,
+        int terminus_idx=1,
+    ) -> int:
+        """
+        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
@@ -176,12 +208,19 @@ cdef class Geod:
         cdef double plat2
         cdef geod_geodesicline line
 
-        cdef array.array array_template = array.array("d", [])
-        cdef array.array out_lons = array.clone(array_template, npts, zero=False)
-        cdef array.array out_lats = array.clone(array_template, npts, zero=False)
-
-        cdef PyBuffWriteManager lats_buff = PyBuffWriteManager(out_lats)
         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:
@@ -204,8 +243,10 @@ cdef class Geod:
                     plon2 *= _DG2RAD
                 lats_buff.data[iii] = plat2
                 lons_buff.data[iii] = plon2
+                if store_az:
+                    azis_buff.data[iii] = pazi2
 
-        return out_lons, out_lats
+        return npts
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
diff --git a/pyproj/geod.py b/pyproj/geod.py
index 11af537c8..dfe61f8e5 100644
--- a/pyproj/geod.py
+++ b/pyproj/geod.py
@@ -375,6 +375,38 @@ def npts(
         )
         return list(zip(lons, lats)) if return_points else (lons, lats)
 
+    def inv_intermediate_by_npts(
+        self,
+        out_lons: Any,
+        out_lats: Any,
+        out_azis: Any,
+        lon1: float,
+        lat1: float,
+        lon2: float,
+        lat2: float,
+        npts: int,
+        initial_idx: int = 1,
+        terminus_idx: int = 1,
+        radians: bool = False,
+    ) -> None:
+        """
+        Similar to npts(), but using the given buffers for the output.
+        out_lons, out_lats, and out_az are the output buffers (out_azis is optional)
+        """
+        super()._inv_intermediate_by_npts(
+            out_lons=out_lons,
+            out_lats=out_lats,
+            out_azis=out_azis,
+            lon1=lon1,
+            lat1=lat1,
+            lon2=lon2,
+            lat2=lat2,
+            npts=npts,
+            radians=radians,
+            initial_idx=initial_idx,
+            terminus_idx=terminus_idx,
+        )
+
     def line_length(self, lons: Any, lats: Any, radians: bool = False) -> float:
         """
         .. versionadded:: 2.3.0
diff --git a/test/test_geod.py b/test/test_geod.py
index 3dcaa511f..a17753637 100644
--- a/test/test_geod.py
+++ b/test/test_geod.py
@@ -7,6 +7,7 @@
 from contextlib import contextmanager
 from itertools import permutations
 
+import numpy as np
 import pytest
 from numpy.testing import assert_almost_equal
 
@@ -132,6 +133,48 @@ def test_geod_inverse_transform():
                 (lat2pt, lon2pt, o_dist), (45.517, -123.683, 832838.542), decimal=3
             )
 
+        if include_initial and include_terminus:
+            lons, lats, azis12, azis21, dists = np.hstack(
+                (lonlats, res_az12_az21_dists)
+            ).transpose()
+
+    lons_a = np.empty(point_count)
+    lats_a = np.empty(point_count)
+    azis_a = np.empty(point_count)
+
+    gg.inv_intermediate_by_npts(
+        out_lons=lons_a,
+        out_lats=lats_a,
+        out_azis=azis_a,
+        lon1=lon1pt,
+        lat1=lat1pt,
+        lon2=lon2pt,
+        lat2=lat2pt,
+        npts=point_count,
+        initial_idx=0,
+        terminus_idx=0,
+    )
+
+    assert_almost_equal(lons_a, lons)
+    assert_almost_equal(lats_a, lats)
+    assert_almost_equal(lats_a, azis12)
+
+    gg.inv_intermediate_by_npts(
+        out_lons=lons_a,
+        out_lats=lats_a,
+        out_azis=None,
+        lon1=lon1pt,
+        lat1=lat1pt,
+        lon2=lon2pt,
+        lat2=lat2pt,
+        npts=point_count,
+        initial_idx=0,
+        terminus_idx=0,
+    )
+
+    assert_almost_equal(lons_a, lons)
+    assert_almost_equal(lats_a, lats)
+
 
 def test_geod_cities():
     # specify the lat/lons of some cities.