From 527c0b5c94893d01e628226d674634efe5a6d31b Mon Sep 17 00:00:00 2001 From: Idan Miara Date: Sun, 9 May 2021 16:05:51 +0300 Subject: [PATCH] pyproj/[geod.py, _geod.pyi, _geod.pyx] - add `inv_intermediate_by_npts` functions that uses given buffers `out_lons`, `out_lats` (and `out_azis`) buffers instead of creating new buffers test/test_geod.py - add tests for `inv_intermediate_by_npts` --- pyproj/_geod.pyi | 14 +++++++++++++ pyproj/_geod.pyx | 53 +++++++++++++++++++++++++++++++++++++++++------ pyproj/geod.py | 32 ++++++++++++++++++++++++++++ test/test_geod.py | 43 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 136 insertions(+), 6 deletions(-) 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.