diff --git a/.all-contributorsrc b/.all-contributorsrc index 329801756..ec25766fd 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -430,6 +430,18 @@ "contributions": [ "platform" ] + }, + { + "login": "idanmiara", + "name": "Idan Miara", + "avatar_url": "https://avatars.githubusercontent.com/u/26349741?v=4", + "profile": "https://github.com/idanmiara", + "contributions": [ + "code", + "doc", + "example", + "test" + ] } ], "contributorsPerLine": 7 diff --git a/README.md b/README.md index a7b2a53bc..a728b29d3 100644 --- a/README.md +++ b/README.md @@ -88,6 +88,7 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d
Poruri Sai Rahul
⚠️
Yann-Sebastien Tremblay-Johnston
📖
odidev
📦 +
Idan Miara
💻 📖 💡 ⚠️ diff --git a/docs/api/enums.rst b/docs/api/enums.rst index ef80bc94b..1af37cda2 100644 --- a/docs/api/enums.rst +++ b/docs/api/enums.rst @@ -15,3 +15,6 @@ Enumerations .. autoclass:: pyproj.enums.PJType :members: + +.. autoclass:: pyproj.enums.GeodIntermediateFlag + :members: diff --git a/docs/api/geod.rst b/docs/api/geod.rst index 9ade0d159..472ac65ac 100644 --- a/docs/api/geod.rst +++ b/docs/api/geod.rst @@ -9,3 +9,6 @@ pyproj.Geod :show-inheritance: :inherited-members: :special-members: __init__ + +.. autoclass:: pyproj.geod.GeodIntermediateReturn + :members: diff --git a/pyproj/_geod.pxd b/pyproj/_geod.pxd index d19788cf7..6ebbd2e8c 100644 --- a/pyproj/_geod.pxd +++ b/pyproj/_geod.pxd @@ -32,6 +32,13 @@ cdef extern from "geodesic.h": double* ps12, double* pazi1, double* pazi2) nogil + void geod_lineinit( + geod_geodesicline* l, + const geod_geodesic* g, + double lat1, + double lon1, + double azi1, + unsigned caps) nogil void geod_inverseline( geod_geodesicline* l, const geod_geodesic* g, diff --git a/pyproj/_geod.pyi b/pyproj/_geod.pyi index a996a0c01..219ead501 100644 --- a/pyproj/_geod.pyi +++ b/pyproj/_geod.pyi @@ -1,7 +1,15 @@ -from typing import Any, Tuple, Type +from typing import Any, NamedTuple, Tuple, Type geodesic_version_str: str +class GeodIntermediateReturn(NamedTuple): + npts: int + del_s: float + dist: float + lons: Any + lats: Any + azis: Any + class Geod: initstring: str a: float @@ -20,15 +28,22 @@ class Geod: def _inv( self, lons1: Any, lats1: Any, lons2: Any, lats2: Any, radians: bool = False ) -> None: ... - def _npts( + def _inv_or_fwd_intermediate( self, lon1: float, lat1: float, - lon2: float, - lat2: float, + lon2_or_azi1: float, + lat2_or_nan: float, npts: int, - radians: bool = False, - ) -> Tuple[Tuple[float], Tuple[float]]: ... + del_s: float, + radians: bool, + initial_idx: int, + terminus_idx: int, + flags: int, + out_lons: Any, + out_lats: Any, + out_azis: Any, + ) -> GeodIntermediateReturn: ... 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 e91573d1c..97655534d 100644 --- a/pyproj/_geod.pyx +++ b/pyproj/_geod.pyx @@ -2,16 +2,62 @@ include "base.pxi" cimport cython from cpython cimport array +from libc.math cimport ceil, isnan, round import array +from collections import namedtuple from pyproj.compat import cstrencode, pystrdecode +from pyproj.enums import GeodIntermediateFlag from pyproj.exceptions import GeodError geodesic_version_str = ( f"{GEODESIC_VERSION_MAJOR}.{GEODESIC_VERSION_MINOR}.{GEODESIC_VERSION_PATCH}" ) +GeodIntermediateReturn = namedtuple( + "GeodIntermediateReturn", ["npts", "del_s", "dist", "lons", "lats", "azis"] +) + +GeodIntermediateReturn.__doc__ = """ +.. versionadded:: 3.1.0 + +Geod Intermediate Return value (Named Tuple) + +Parameters +---------- + +npts: int + number of points +del_s: float + delimiter distance between two successive points +dist: float + distance between the initial and terminus points +out_lons: Any + array of the output lons +out_lats: Any + array of the output lats +out_azis: Any + array of the output azis +""" + + +cdef int GEOD_INTER_FLAG_DEFAULT = GeodIntermediateFlag.DEFAULT + +cdef int GEOD_INTER_FLAG_NPTS_MASK = GeodIntermediateFlag.NPTS_MASK +cdef int GEOD_INTER_FLAG_NPTS_ROUND = GeodIntermediateFlag.NPTS_ROUND +cdef int GEOD_INTER_FLAG_NPTS_CEIL = GeodIntermediateFlag.NPTS_CEIL +cdef int GEOD_INTER_FLAG_NPTS_TRUNC = GeodIntermediateFlag.NPTS_TRUNC + +cdef int GEOD_INTER_FLAG_DEL_S_MASK = GeodIntermediateFlag.DEL_S_MASK +cdef int GEOD_INTER_FLAG_DEL_S_RECALC = GeodIntermediateFlag.DEL_S_RECALC +cdef int GEOD_INTER_FLAG_DEL_S_NO_RECALC = GeodIntermediateFlag.DEL_S_NO_RECALC + +cdef int GEOD_INTER_FLAG_AZIS_MASK = GeodIntermediateFlag.AZIS_MASK +cdef int GEOD_INTER_FLAG_AZIS_DISCARD = GeodIntermediateFlag.AZIS_DISCARD +cdef int GEOD_INTER_FLAG_AZIS_KEEP = GeodIntermediateFlag.AZIS_KEEP + + cdef class Geod: def __init__(self, double a, double f, bint sphere, double b, double es): geod_init(&self._geod_geodesic, a, f) @@ -153,54 +199,117 @@ cdef class Geod: @cython.boundscheck(False) @cython.wraparound(False) - def _npts( + def _inv_or_fwd_intermediate( self, double lon1, double lat1, - double lon2, - double lat2, + double lon2_or_azi1, + double lat2_or_nan, int npts, - bint radians=False, - ): + double del_s, + bint radians, + int initial_idx, + int terminus_idx, + int flags, + object out_lons, + object out_lats, + object out_azis, + ) -> GeodIntermediateReturn: """ + .. versionadded:: 3.1.0 + 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 bint store_az = \ + out_azis is not None \ + or (flags & GEOD_INTER_FLAG_AZIS_MASK) == GEOD_INTER_FLAG_AZIS_KEEP + cdef array.array array_template = array.array("d", []) - cdef array.array lats = array.clone(array_template, npts, zero=False) - cdef array.array lons = array.clone(array_template, npts, zero=False) - cdef PyBuffWriteManager lats_buff = PyBuffWriteManager(lats) - cdef PyBuffWriteManager lons_buff = PyBuffWriteManager(lons) + cdef PyBuffWriteManager lons_buff + cdef PyBuffWriteManager lats_buff + cdef PyBuffWriteManager azis_buff + + cdef bint is_fwd = isnan(lat2_or_nan) + if not is_fwd and (del_s == 0) == (npts == 0): + raise GeodError("inv_intermediate: " + "npts and del_s are mutually exclusive, " + "only one of them must be != 0.") with nogil: if radians: lon1 *= _RAD2DG lat1 *= _RAD2DG - lon2 *= _RAD2DG - lat2 *= _RAD2DG - # do inverse computation to set azimuths, distance. - geod_inverseline( - &line, &self._geod_geodesic, lat1, lon1, lat2, lon2, 0u, - ) - # distance increment. - del_s = line.s13 / (npts + 1) + lon2_or_azi1 *= _RAD2DG + if not is_fwd: + lat2_or_nan *= _RAD2DG + + if is_fwd: + # do fwd computation to set azimuths, distance. + geod_lineinit(&line, &self._geod_geodesic, lat1, lon1, lon2_or_azi1, 0u) + line.s13 = del_s * (npts + initial_idx + terminus_idx - 1) + else: + # do inverse computation to set azimuths, distance. + geod_inverseline(&line, &self._geod_geodesic, lat1, lon1, + lat2_or_nan, lon2_or_azi1, 0u) + + if npts == 0: + # calc the number of required points by the distance increment + # s12 holds a temporary float value of npts (just reusing this var) + s12 = line.s13 / del_s - initial_idx - terminus_idx + 1 + if (flags & GEOD_INTER_FLAG_NPTS_MASK) == \ + GEOD_INTER_FLAG_NPTS_ROUND: + s12 = round(s12) + elif (flags & GEOD_INTER_FLAG_NPTS_MASK) == \ + GEOD_INTER_FLAG_NPTS_CEIL: + s12 = ceil(s12) + npts = int(s12) + if (flags & GEOD_INTER_FLAG_DEL_S_MASK) == GEOD_INTER_FLAG_DEL_S_RECALC: + # calc the distance increment by the number of required points + del_s = line.s13 / (npts + initial_idx + terminus_idx - 1) + + with gil: + if out_lons is None: + out_lons = array.clone(array_template, npts, zero=False) + if out_lats is None: + out_lats = array.clone(array_template, npts, zero=False) + if out_azis is None and store_az: + out_azis = array.clone(array_template, npts, zero=False) + + lons_buff = PyBuffWriteManager(out_lons) + lats_buff = PyBuffWriteManager(out_lats) + if store_az: + azis_buff = PyBuffWriteManager(out_azis) + + if lons_buff.len < npts \ + or lats_buff.len < npts \ + or (store_az and azis_buff.len < npts): + raise GeodError( + "Arrays are not long enough (" + f"{lons_buff.len}, {lats_buff.len}, " + f"{azis_buff.len if store_az else -1}) < {npts}.") + # loop over intermediate points, compute lat/lons. - for iii in range(1, npts + 1): - s12 = iii * del_s + for iii in range(0, npts): + s12 = (iii + initial_idx) * del_s geod_position(&line, s12, &plat2, &plon2, &pazi2) if radians: plat2 *= _DG2RAD plon2 *= _DG2RAD - lats_buff.data[iii - 1] = plat2 - lons_buff.data[iii - 1] = plon2 + lats_buff.data[iii] = plat2 + lons_buff.data[iii] = plon2 + if store_az: + azis_buff.data[iii] = pazi2 - return lons, lats + return GeodIntermediateReturn( + npts, del_s, line.s13, out_lons, out_lats, out_azis) @cython.boundscheck(False) @cython.wraparound(False) diff --git a/pyproj/enums.py b/pyproj/enums.py index 5caea79c4..d44b9bf5b 100644 --- a/pyproj/enums.py +++ b/pyproj/enums.py @@ -1,7 +1,7 @@ """ This module contains enumerations used in pyproj. """ -from enum import Enum +from enum import Enum, IntFlag class BaseEnum(Enum): @@ -136,3 +136,26 @@ class PJType(BaseEnum): TRANSFORMATION = "TRANSFORMATION" CONCATENATED_OPERATION = "CONCATENATED_OPERATION" OTHER_COORDINATE_OPERATION = "OTHER_COORDINATE_OPERATION" + + +class GeodIntermediateFlag(IntFlag): + """ + .. versionadded:: 3.1.0 + + Flags to be used in Geod.[inv|fwd]_intermediate() + """ + + DEFAULT = 0x0 + + NPTS_MASK = 0xF + NPTS_ROUND = 0x0 + NPTS_CEIL = 0x1 + NPTS_TRUNC = 0x2 + + DEL_S_MASK = 0xF0 + DEL_S_RECALC = 0x00 + DEL_S_NO_RECALC = 0x10 + + AZIS_MASK = 0xF00 + AZIS_DISCARD = 0x000 + AZIS_KEEP = 0x100 diff --git a/pyproj/geod.py b/pyproj/geod.py index 79b027bda..94d52471e 100644 --- a/pyproj/geod.py +++ b/pyproj/geod.py @@ -8,13 +8,20 @@ latitudes and longitudes of an initial and terminus point. """ -__all__ = ["Geod", "pj_ellps", "geodesic_version_str"] +__all__ = [ + "Geod", + "pj_ellps", + "geodesic_version_str", + "GeodIntermediateFlag", + "GeodIntermediateReturn", +] import math from typing import Any, Dict, List, Optional, Tuple, Union from pyproj._geod import Geod as _Geod -from pyproj._geod import geodesic_version_str +from pyproj._geod import GeodIntermediateReturn, geodesic_version_str +from pyproj.enums import GeodIntermediateFlag from pyproj.exceptions import GeodError from pyproj.list import get_ellps_map from pyproj.utils import _convertback, _copytobuffer @@ -273,13 +280,19 @@ def npts( lat2: float, npts: int, radians: bool = False, + initial_idx: int = 1, + terminus_idx: int = 1, ) -> List: """ + .. versionadded:: 3.1.0 initial_idx, terminus_idx + Given a single initial point and terminus point, returns a list of longitude/latitude pairs describing npts equally spaced intermediate points along the geodesic between the initial and terminus points. + Similar to inv_intermediate(), but with less options. + Example usage: >>> from pyproj import Geod @@ -336,17 +349,312 @@ def npts( Latitude of the terminus point npts: int Number of points to be returned + (including initial and/or terminus points, if required) radians: bool, optional If True, the input data is assumed to be in radians. - + initial_idx: int, optional + if initial_idx==0 then the initial point would be included in the output + (as the first point) + terminus_idx: int, optional + if terminus_idx==0 then the terminus point would be included in the output + (as the last point) Returns ------- list of tuples: list of (lon, lat) points along the geodesic between the initial and terminus points. """ - lons, lats = super()._npts(lon1, lat1, lon2, lat2, npts, radians=radians) - return list(zip(lons, lats)) + + res = self._inv_or_fwd_intermediate( + lon1=lon1, + lat1=lat1, + lon2_or_azi1=lon2, + lat2_or_nan=lat2, + npts=npts, + del_s=0, + radians=radians, + initial_idx=initial_idx, + terminus_idx=terminus_idx, + flags=GeodIntermediateFlag.AZIS_DISCARD, + out_lons=None, + out_lats=None, + out_azis=None, + ) + return list(zip(res.lons, res.lats)) + + def inv_intermediate( + self, + lon1: float, + lat1: float, + lon2: float, + lat2: float, + npts: int = 0, + del_s: float = 0, + initial_idx: int = 1, + terminus_idx: int = 1, + radians: bool = False, + flags: GeodIntermediateFlag = GeodIntermediateFlag.DEFAULT, + out_lons: Any = None, + out_lats: Any = None, + out_azis: Any = None, + ) -> GeodIntermediateReturn: + """ + .. versionadded:: 3.1.0 + + Given a single initial point and terminus point, + and the number of points, returns + a list of longitude/latitude pairs describing npts equally + spaced intermediate points along the geodesic between the + initial and terminus points. + + npts and del_s parameters are mutually exclusive: + + if npts != 0: + it calculates the distance between the points by + the distance between the initial point and the + terminus point divided by npts + (the number of intermediate points) + else: + it calculates the number of intermediate points by + dividing the distance between the initial and + terminus points by del_s + (delimiter distance between two successive points) + + Similar to npts(), but with more options. + + Example usage: + + >>> from pyproj import Geod + >>> g = Geod(ellps='clrk66') # Use Clarke 1866 ellipsoid. + >>> # specify the lat/lons of Boston and Portland. + >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.) + >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.) + >>> # find ten equally spaced points between Boston and Portland. + >>> r = g.inv_intermediate(boston_lon,boston_lat,portland_lon,portland_lat,10) + >>> for lon,lat in zip(r.lons, r.lats): f'{lat:.3f} {lon:.3f}' + '43.528 -75.414' + '44.637 -79.883' + '45.565 -84.512' + '46.299 -89.279' + '46.830 -94.156' + '47.149 -99.112' + '47.251 -104.106' + '47.136 -109.100' + '46.805 -114.051' + '46.262 -118.924' + >>> # test with radians=True (inputs/outputs in radians, not degrees) + >>> import math + >>> dg2rad = math.radians(1.) + >>> rad2dg = math.degrees(1.) + >>> r = g.inv_intermediate( + ... dg2rad*boston_lon, + ... dg2rad*boston_lat, + ... dg2rad*portland_lon, + ... dg2rad*portland_lat, + ... 10, + ... radians=True + ... ) + >>> for lon,lat in zip(r.lons, r.lats): f'{rad2dg*lat:.3f} {rad2dg*lon:.3f}' + '43.528 -75.414' + '44.637 -79.883' + '45.565 -84.512' + '46.299 -89.279' + '46.830 -94.156' + '47.149 -99.112' + '47.251 -104.106' + '47.136 -109.100' + '46.805 -114.051' + '46.262 -118.924' + + Parameters + ---------- + lon1: float + Longitude of the initial point + lat1: float + Latitude of the initial point + lon2: float + Longitude of the terminus point + lat2: float + Latitude of the terminus point + npts: int, optional + Number of points to be returned + npts == 0 iff del_s != 0 + del_s: float, optional + delimiter distance between two successive points + del_s == 0 iff npts != 0 + radians: bool, optional + If True, the input data is assumed to be in radians. + initial_idx: int, optional + if initial_idx==0 then the initial point would be included in the output + (as the first point) + terminus_idx: int, optional + if terminus_idx==0 then the terminus point would be included in the output + (as the last point) + flags: GeodIntermediateFlag, optional + 1st: round/ceil/trunc (see GeodIntermediateFlag.NPTS_*) + 2nd: update del_s to the new npts or not (see GeodIntermediateFlag.DEL_S_*) + 3rd: iff out_azis=None, indicates if to save or discard the azimuths + (see GeodIntermediateFlag.AZIS_*) + default - round npts, update del_s accordingly, discard azis + out_lons: array, :class:`numpy.ndarray`, optional + Longitude(s) of the intermediate point(s) + If None then buffers would be allocated internnaly + out_lats: array, :class:`numpy.ndarray`, optional + Latitudes(s) of the intermediate point(s) + If None then buffers would be allocated internnaly + out_azis: array, :class:`numpy.ndarray`, optional + az12(s) of the intermediate point(s) + If None then buffers would be allocated internnaly + unless requested otherwise by the flags + + Returns + ------- + GeodIntermediateReturn + number of points, distance and output arrays (GeodIntermediateReturn docs) + """ + return super()._inv_or_fwd_intermediate( + lon1=lon1, + lat1=lat1, + lon2_or_azi1=lon2, + lat2_or_nan=lat2, + npts=npts, + del_s=del_s, + radians=radians, + initial_idx=initial_idx, + terminus_idx=terminus_idx, + flags=int(flags), + out_lons=out_lons, + out_lats=out_lats, + out_azis=out_azis, + ) + + def fwd_intermediate( + self, + lon1: float, + lat1: float, + azi1: float, + npts: int, + del_s: float, + initial_idx: int = 1, + terminus_idx: int = 1, + radians: bool = False, + flags: GeodIntermediateFlag = GeodIntermediateFlag.DEFAULT, + out_lons: Any = None, + out_lats: Any = None, + out_azis: Any = None, + ) -> GeodIntermediateReturn: + """ + .. versionadded:: 3.1.0 + + Given a single initial point and azimuth, number of points (npts) + and delimiter distance between two successive points (del_s), returns + a list of longitude/latitude pairs describing npts equally + spaced intermediate points along the geodesic between the + initial and terminus points. + + Example usage: + + >>> from pyproj import Geod + >>> g = Geod(ellps='clrk66') # Use Clarke 1866 ellipsoid. + >>> # specify the lat/lons of Boston and Portland. + >>> boston_lat = 42.+(15./60.); boston_lon = -71.-(7./60.) + >>> portland_lat = 45.+(31./60.); portland_lon = -123.-(41./60.) + >>> az12,az21,dist = g.inv(boston_lon,boston_lat,portland_lon,portland_lat) + >>> # find ten equally spaced points between Boston and Portland. + >>> npts = 10 + >>> del_s = dist/(npts+1) + >>> r = g.fwd_intermediate(boston_lon,boston_lat,az12,npts=npts,del_s=del_s) + >>> for lon,lat in zip(r.lons, r.lats): f'{lat:.3f} {lon:.3f}' + '43.528 -75.414' + '44.637 -79.883' + '45.565 -84.512' + '46.299 -89.279' + '46.830 -94.156' + '47.149 -99.112' + '47.251 -104.106' + '47.136 -109.100' + '46.805 -114.051' + '46.262 -118.924' + >>> # test with radians=True (inputs/outputs in radians, not degrees) + >>> import math + >>> dg2rad = math.radians(1.) + >>> rad2dg = math.degrees(1.) + >>> r = g.fwd_intermediate( + ... dg2rad*boston_lon, + ... dg2rad*boston_lat, + ... dg2rad*az12, + ... npts=npts, + ... del_s=del_s, + ... radians=True + ... ) + >>> for lon,lat in zip(r.lons, r.lats): f'{rad2dg*lat:.3f} {rad2dg*lon:.3f}' + '43.528 -75.414' + '44.637 -79.883' + '45.565 -84.512' + '46.299 -89.279' + '46.830 -94.156' + '47.149 -99.112' + '47.251 -104.106' + '47.136 -109.100' + '46.805 -114.051' + '46.262 -118.924' + + Parameters + ---------- + lon1: float + Longitude of the initial point + lat1: float + Latitude of the initial point + azi1: float + Azimuth from the initial point towards the terminus point + npts: int + Number of points to be returned + (including initial and/or terminus points, if required) + del_s: float + delimiter distance between two successive points + radians: bool, optional + If True, the input data is assumed to be in radians. + initial_idx: int, optional + if initial_idx==0 then the initial point would be included in the output + (as the first point) + terminus_idx: int, optional + if terminus_idx==0 then the terminus point would be included in the output + (as the last point) + flags: GeodIntermediateFlag, optional + 3rd: iff out_azis=None, indicates if to save or discard the azimuths + (see GeodIntermediateFlag.AZIS_*) + default - discard azis + out_lons: array, :class:`numpy.ndarray`, optional + Longitude(s) of the intermediate point(s) + If None then buffers would be allocated internnaly + out_lats: array, :class:`numpy.ndarray`, optional + Latitudes(s) of the intermediate point(s) + If None then buffers would be allocated internnaly + out_azis: array, :class:`numpy.ndarray`, optional + az12(s) of the intermediate point(s) + If None then buffers would be allocated internnaly + unless requested otherwise by the flags + + Returns + ------- + GeodIntermediateReturn + number of points, distance and output arrays (GeodIntermediateReturn docs) + """ + return super()._inv_or_fwd_intermediate( + lon1=lon1, + lat1=lat1, + lon2_or_azi1=azi1, + lat2_or_nan=math.nan, + npts=npts, + del_s=del_s, + radians=radians, + initial_idx=initial_idx, + terminus_idx=terminus_idx, + flags=int(flags), + out_lons=out_lons, + out_lats=out_lats, + out_azis=out_azis, + ) def line_length(self, lons: Any, lats: Any, radians: bool = False) -> float: """ diff --git a/test/test_geod.py b/test/test_geod.py index 04ebd73e2..d3b409789 100644 --- a/test/test_geod.py +++ b/test/test_geod.py @@ -1,3 +1,4 @@ +import itertools import math import os import pickle @@ -6,10 +7,12 @@ from contextlib import contextmanager from itertools import permutations +import numpy as np import pytest from numpy.testing import assert_almost_equal from pyproj import Geod +from pyproj.geod import GeodIntermediateFlag try: from shapely.geometry import ( @@ -64,37 +67,216 @@ def test_geod_inverse_transform(): b'-66.531\t75.654\t4164192.708\n' """ + true_az12 = -66.5305947876623 + true_az21 = 75.65363415556968 print("from pyproj.Geod.inv:") az12, az21, dist = gg.inv(lon1pt, lat1pt, lon2pt, lat2pt) - assert_almost_equal((az12, az21, dist), (-66.531, 75.654, 4164192.708), decimal=3) + assert_almost_equal( + (az12, az21, dist), (true_az12, true_az21, 4164192.708), decimal=3 + ) print("forward transform") print("from proj.4 geod:") endlon, endlat, backaz = gg.fwd(lon1pt, lat1pt, az12, dist) - assert_almost_equal((endlon, endlat, backaz), (-123.683, 45.517, 75.654), decimal=3) - print("intermediate points:") - print("from geod with +lat_1,+lon_1,+lat_2,+lon_2,+n_S:") - npts = 4 - lonlats = gg.npts(lon1pt, lat1pt, lon2pt, lat2pt, npts) - lonprev = lon1pt - latprev = lat1pt - print(dist / (npts + 1)) - print("%6.3f %7.3f" % (lat1pt, lon1pt)) - result_dists = ( + assert_almost_equal( + (endlon, endlat, backaz), (lon2pt, lat2pt, true_az21), decimal=3 + ) + + inc_exc = ["excluding", "including"] + res_az12_az21_dists_all = [ + (180.0, 0.0, 0.0), (-66.53059478766238, 106.79071710136431, 832838.5416198927), (-73.20928289863558, 99.32289055927389, 832838.5416198935), (-80.67710944072617, 91.36325611787134, 832838.5416198947), (-88.63674388212858, 83.32809401477382, 832838.5416198922), + (-96.67190598522616, 75.65363415556973, 832838.5416198926), + ] + point_count = len(res_az12_az21_dists_all) + for include_initial, include_terminus in itertools.product( + (False, True), (False, True) + ): + initial_idx = int(not include_initial) + terminus_idx = int(not include_terminus) + + npts = point_count - initial_idx - terminus_idx + print("intermediate points:") + print("from geod with +lat_1,+lon_1,+lat_2,+lon_2,+n_S:") + print(f"{lat1pt:6.3f} {lon1pt:7.3f} {lat2pt:6.3f} {lon2pt:7.3f} {npts}") + + lonlats = gg.npts( + lon1pt, + lat1pt, + lon2pt, + lat2pt, + npts, + initial_idx=initial_idx, + terminus_idx=terminus_idx, + ) + assert len(lonlats) == npts + + npts1 = npts + initial_idx + terminus_idx - 1 + del_s = dist / npts1 + print( + f"Total distnace is {dist}, " + f"Points count: {npts}, " + f"{inc_exc[include_initial]} initial point, " + f"{inc_exc[include_terminus]} terminus point. " + f"The distance between successive points is {dist}/{npts1} = {del_s}" + ) + + from_idx = initial_idx + to_idx = point_count - terminus_idx + res_az12_az21_dists = res_az12_az21_dists_all[from_idx:to_idx] + + lonprev = lon1pt + latprev = lat1pt + for (lon, lat), (res12, res21, resdist) in zip(lonlats, res_az12_az21_dists): + o_az12, o_az21, o_dist = gg.inv(lonprev, latprev, lon, lat) + assert_almost_equal((o_az12, o_az21, o_dist), (res12, res21, resdist)) + latprev = lat + lonprev = lon + if not include_terminus: + o_az12, o_az21, o_dist = gg.inv(lonprev, latprev, lon2pt, lat2pt) + assert_almost_equal( + (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() + + del_s = dist / (point_count - 1) + lons_a = np.empty(point_count) + lats_a = np.empty(point_count) + azis_a = np.empty(point_count) + + print("test inv_intermediate (by npts) with azi output") + res = gg.inv_intermediate( + 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, ) - for (lon, lat), (res12, res21, resdist) in zip(lonlats, result_dists): - az12, az21, dist = gg.inv(lonprev, latprev, lon, lat) - assert_almost_equal((az12, az21, dist), (res12, res21, resdist)) - latprev = lat - lonprev = lon - az12, az21, dist = gg.inv(lonprev, latprev, lon2pt, lat2pt) - assert_almost_equal( - (lat2pt, lon2pt, dist), (45.517, -123.683, 832838.542), decimal=3 + assert res.npts == point_count + assert_almost_equal(res.del_s, del_s) + assert_almost_equal(res.dist, dist) + assert_almost_equal(res.lons, lons) + assert_almost_equal(res.lats, lats) + assert_almost_equal(res.azis[:-1], azis12[1:]) + assert res.lons is lons_a + assert res.lats is lats_a + assert res.azis is azis_a + + for flags in (GeodIntermediateFlag.AZIS_DISCARD, GeodIntermediateFlag.AZIS_KEEP): + print("test inv_intermediate (by npts) without azi output, no buffers") + res = gg.inv_intermediate( + lon1=lon1pt, + lat1=lat1pt, + lon2=lon2pt, + lat2=lat2pt, + npts=point_count, + initial_idx=0, + terminus_idx=0, + flags=flags, + ) + assert res.npts == point_count + assert_almost_equal(res.del_s, del_s) + assert_almost_equal(res.dist, dist) + assert_almost_equal(res.lons, lons_a) + assert_almost_equal(res.lats, lats_a) + if flags == GeodIntermediateFlag.AZIS_DISCARD: + assert res.azis is None + else: + assert_almost_equal(res.azis, azis_a) + + lons_b = np.empty(point_count) + lats_b = np.empty(point_count) + azis_b = np.empty(point_count) + + print("test inv_intermediate (by npts) without azi output") + res = gg.inv_intermediate( + out_lons=lons_b, + out_lats=lats_b, + out_azis=None, + lon1=lon1pt, + lat1=lat1pt, + lon2=lon2pt, + lat2=lat2pt, + npts=point_count, + initial_idx=0, + terminus_idx=0, + flags=flags, + ) + assert res.npts == point_count + assert_almost_equal(res.del_s, del_s) + assert_almost_equal(res.dist, dist) + assert_almost_equal(res.lons, lons_a) + assert_almost_equal(res.lats, lats_a) + assert res.lons is lons_b + assert res.lats is lats_b + if flags == GeodIntermediateFlag.AZIS_DISCARD: + assert res.azis is None + else: + assert_almost_equal(res.azis, azis_a) + + print("test fwd_intermediate") + res = gg.fwd_intermediate( + out_lons=lons_b, + out_lats=lats_b, + out_azis=azis_b, + lon1=lon1pt, + lat1=lat1pt, + azi1=true_az12, + npts=point_count, + del_s=del_s, + initial_idx=0, + terminus_idx=0, ) + assert res.npts == point_count + assert_almost_equal(res.del_s, del_s) + assert_almost_equal(res.dist, dist) + assert_almost_equal(res.lons, lons_a) + assert_almost_equal(res.lats, lats_a) + assert_almost_equal(res.azis, azis_a) + assert res.lons is lons_b + assert res.lats is lats_b + assert res.azis is azis_b + + print("test inv_intermediate (by del_s)") + for del_s_fact, flags in ( + (1, GeodIntermediateFlag.NPTS_ROUND), + ((point_count - 0.5) / point_count, GeodIntermediateFlag.NPTS_TRUNC), + ((point_count + 0.5) / point_count, GeodIntermediateFlag.NPTS_CEIL), + ): + res = gg.inv_intermediate( + out_lons=lons_b, + out_lats=lats_b, + out_azis=azis_b, + lon1=lon1pt, + lat1=lat1pt, + lon2=lon2pt, + lat2=lat2pt, + del_s=del_s * del_s_fact, + initial_idx=0, + terminus_idx=0, + flags=flags, + ) + assert res.npts == point_count + assert_almost_equal(res.del_s, del_s) + assert_almost_equal(res.dist, dist) + assert_almost_equal(res.lons, lons_a) + assert_almost_equal(res.lats, lats_a) + assert_almost_equal(res.azis, azis_a) + assert res.lons is lons_b + assert res.lats is lats_b + assert res.azis is azis_b def test_geod_cities():