From dcc2af80af068cfb40542150b316aea82b069da8 Mon Sep 17 00:00:00 2001 From: Greg Lucas Date: Fri, 16 Dec 2022 10:46:41 -0700 Subject: [PATCH] ENH: Add single point versions of transform, fwd, inv When transforming single points, this can add a significant speed benefit. --- pyproj/geod.py | 109 +++++++++++++++++++++++++++++++++++++++ pyproj/transformer.py | 103 ++++++++++++++++++++++++++++++++++++ test/test_geod.py | 17 ++++++ test/test_transformer.py | 11 ++++ 4 files changed, 240 insertions(+) diff --git a/pyproj/geod.py b/pyproj/geod.py index ed4b24fb0..b114d424f 100644 --- a/pyproj/geod.py +++ b/pyproj/geod.py @@ -17,6 +17,7 @@ ] import math +from array import array from typing import Any, Dict, List, Optional, Tuple, Union from pyproj._geod import Geod as _Geod @@ -310,6 +311,61 @@ def fwd( # pylint: disable=invalid-name outz = _convertback(z_data_type, inz) return outx, outy, outz + def fwd_point( # pylint: disable=invalid-name + self, + lon: float, + lat: float, + az: float, + dist: float, + radians: bool = False, + ) -> Tuple[float, float, float]: + """ + Forward transformation + + Determine longitude, latitude and back azimuth of terminus + point given longitude and latitude of initial point, + plus forward azimuth and distance. + + Accepted numeric scalar: + + - :class:`int` + - :class:`float` + - :class:`numpy.floating` + - :class:`numpy.integer` + + Parameters + ---------- + lon: scalar + Longitude of initial point + lat: scalar + Latitude of initial point + az: scalar + Forward azimuth + dist: scalar + Distance between initial and terminus point + in meters + radians: bool, default=False + If True, the input data is assumed to be in radians. + Otherwise, the data is assumed to be in degrees. + + Returns + ------- + scalar: + Longitude of terminus point + scalar: + Latitude of terminus point + scalar: + Back azimuth + """ + # process inputs, making copies that support buffer API. + lons = array("d", (float(lon),)) + lats = array("d", (float(lat),)) + azs = array("d", (float(az),)) + dists = array("d", (float(dist),)) + self._fwd(lons, lats, azs, dists, radians=radians) + # Output was computed inplace + return lons[0], lats[0], azs[0] + def inv( self, lons1: Any, @@ -380,6 +436,59 @@ def inv( outz = _convertback(z_data_type, inz) return outx, outy, outz + def inv_point( + self, + lon1: float, + lat1: float, + lon2: float, + lat2: float, + radians: bool = False, + ) -> Tuple[float, float, float]: + """ + Inverse transformation + + Determine forward and back azimuth, plus distance + between initial point and terminus point. + + Accepted numeric scalar or array: + + - :class:`int` + - :class:`float` + - :class:`numpy.floating` + - :class:`numpy.integer` + + Parameters + ---------- + lon1: scalar + Longitude of initial point + lat1: scalar + Latitude of initial point + lon2: scalar + Longitude of terminus point + lat2: scalar + Latitude of terminus point + radians: bool, default=False + If True, the input data is assumed to be in radians. + Otherwise, the data is assumed to be in degrees. + + Returns + ------- + scalar: + Forward azimuth + scalar: + Back azimuth + scalar: + Distance between initial and terminus point in meters + """ + # process inputs, making copies that support buffer API. + lon1s = array("d", (float(lon1),)) + lat1s = array("d", (float(lat1),)) + lon2s = array("d", (float(lon2),)) + lat2s = array("d", (float(lat2),)) + self._inv(lon1s, lat1s, lon2s, lat2s, radians=radians) + # Output was computed inplace + return lon1s[0], lat1s[0], lon2s[0] + def npts( self, lon1: float, diff --git a/pyproj/transformer.py b/pyproj/transformer.py index 0817746c8..bae787bb9 100644 --- a/pyproj/transformer.py +++ b/pyproj/transformer.py @@ -1,6 +1,7 @@ """ The transformer module is for performing cartographic transformations. """ +# pylint: disable=too-many-lines __all__ = [ "transform", "itransform", @@ -820,6 +821,108 @@ def transform( # pylint: disable=invalid-name return_data += (_convertback(t_data_type, intime),) return return_data + @overload + def transform_point( # pylint: disable=invalid-name + self, + x: float, + y: float, + radians: bool = False, + errcheck: bool = False, + direction: Union[TransformDirection, str] = TransformDirection.FORWARD, + ) -> Tuple[float, float]: + ... + + @overload + def transform_point( # pylint: disable=invalid-name + self, + x: float, + y: float, + z: float, + radians: bool = False, + errcheck: bool = False, + direction: Union[TransformDirection, str] = TransformDirection.FORWARD, + ) -> Tuple[float, float, float]: + ... + + @overload + def transform_point( # pylint: disable=invalid-name + self, + x: float, + y: float, + z: float, + t: float, + radians: bool = False, + errcheck: bool = False, + direction: Union[TransformDirection, str] = TransformDirection.FORWARD, + ) -> Tuple[float, float, float, float]: + ... + + def transform_point( # pylint: disable=invalid-name + self, + x, + y, + z=None, + t=None, + radians=False, + errcheck=False, + direction=TransformDirection.FORWARD, + ): + """ + Transform a single point between two coordinate systems. + + See: :c:func:`proj_trans_generic` + + Accepted numeric scalar or array: + + - :class:`int` + - :class:`float` + - :class:`numpy.floating` + - :class:`numpy.integer` + + Parameters + ---------- + x: scalar + Input x coordinate. + y: scalar + Input y coordinate. + z: scalar, optional + Input z coordinate. + t: scalar, optional + Input time coordinate. + radians: bool, default=False + If True, will expect input data to be in radians and will return radians + if the projection is geographic. Otherwise, it uses degrees. + errcheck: bool, default=False + If True, an exception is raised if the errors are found in the process. + If False, ``inf`` is returned for errors. + direction: pyproj.enums.TransformDirection, optional + The direction of the transform. + Default is :attr:`pyproj.enums.TransformDirection.FORWARD`. + """ + # process inputs, making copies that support buffer API. + x = array("d", (float(x),)) + y = array("d", (float(y),)) + if z is not None: + z = array("d", (float(z),)) + if t is not None: + t = array("d", (float(t),)) + # call pj_transform. xx, yy, zz buffers modified in place. + self._transformer._transform( + x, + y, + inz=z, + intime=t, + direction=direction, + radians=radians, + errcheck=errcheck, + ) + out = (x[0], y[0]) + if z is not None: + out += (z[0],) + if t is not None: + out += (t[0],) + return out + def itransform( self, points: Any, diff --git a/test/test_geod.py b/test/test_geod.py index 33cc1ad30..8e45a2df9 100644 --- a/test/test_geod.py +++ b/test/test_geod.py @@ -327,6 +327,23 @@ def test_geod_cities(): assert_almost_equal((az12, az21, dist), (-66.531, 75.654, 4164192.708), decimal=3) +@pytest.mark.parametrize("radians", [(True, False)]) +def test_geod_singlepoint(radians): + x1 = 0 + y1 = 0 + x2 = 1 + y2 = 1 + g1 = Geod(ellps="clrk66") + assert_almost_equal( + g1.fwd(x1, y1, x2, y2, radians=radians), + g1.fwd_point(x1, y1, x2, y2, radians=radians), + ) + assert_almost_equal( + g1.inv(x1, y1, x2, y2, radians=radians), + g1.inv_point(x1, y1, x2, y2, radians=radians), + ) + + def test_line_length__single_point(): geod = Geod(ellps="WGS84") assert geod.line_length(1, 1) == 0 diff --git a/test/test_transformer.py b/test/test_transformer.py index 8cc6d7d2b..3aa66b7eb 100644 --- a/test/test_transformer.py +++ b/test/test_transformer.py @@ -1505,6 +1505,17 @@ def test_4d_transform__inplace__numpy__int(): assert tarr[0] == 2019 +@pytest.mark.parametrize("z,t", [(None, None), (None, 0), (0, None)]) +def test_transform_single_point(z, t): + transformer = Transformer.from_crs(7789, 8401) + x = 0 + y = 0 + assert_almost_equal( + transformer.transform(xx=x, yy=y, zz=z, tt=t), + transformer.transform_point(x=x, y=y, z=z, t=t), + ) + + def test_transformer_source_target_crs(): transformer = Transformer.from_crs("EPSG:4326", "EPSG:4258") assert transformer.source_crs == "EPSG:4326"