[Git][debian-gis-team/pymap3d][upstream] New upstream version 2.6.0
Antonio Valentino
gitlab at salsa.debian.org
Sun Mar 21 11:32:25 GMT 2021
Antonio Valentino pushed to branch upstream at Debian GIS Project / pymap3d
Commits:
486084e0 by Antonio Valentino at 2021-03-21T11:27:56+00:00
New upstream version 2.6.0
- - - - -
22 changed files:
- .github/workflows/ci.yml
- .github/workflows/ci_stdlib_only.yml
- README.md
- pyproject.toml
- + scripts/benchmark_ecef2geo.py
- scripts/benchmark_vincenty.py
- setup.cfg
- src/pymap3d/aer.py
- src/pymap3d/ecef.py
- src/pymap3d/enu.py
- src/pymap3d/latitude.py
- src/pymap3d/los.py
- src/pymap3d/lox.py
- src/pymap3d/ned.py
- src/pymap3d/rcurve.py
- src/pymap3d/rsphere.py
- + src/pymap3d/tests/__init__.py
- src/pymap3d/tests/test_aer.py
- src/pymap3d/tests/test_geodetic.py
- src/pymap3d/tests/test_vincenty_vreckon.py
- src/pymap3d/utils.py
- src/pymap3d/vincenty.py
Changes:
=====================================
.github/workflows/ci.yml
=====================================
@@ -4,6 +4,7 @@ on:
push:
paths:
- "**.py"
+ - .github/workflows/ci.yml
jobs:
=====================================
.github/workflows/ci_stdlib_only.yml
=====================================
@@ -4,6 +4,7 @@ on:
push:
paths:
- "**.py"
+ - .github/workflows/ci_stdlib_only.yml
jobs:
=====================================
README.md
=====================================
@@ -2,7 +2,6 @@
[](https://doi.org/10.5281/zenodo.213676)
[](https://doi.org/10.21105/joss.00580)
-[](https://github.com/Ullaakut/astronomer/)
[](https://lgtm.com/projects/g/geospace-code/pymap3d/context:python)
=====================================
pyproject.toml
=====================================
@@ -1,8 +1,44 @@
-[build-system]
-requires = ["setuptools", "wheel"]
+[project]
+name = "pymap3d"
+version = "2.5.1"
+description = "pure Python (no prereqs) coordinate conversions, following convention of several popular Matlab routines."
+readme = "README.md"
+requires-python = ">=3.7"
+license = {file = "LICENSE.txt"}
+keywords = ["geodesy"]
+classifiers = [
+ 'Development Status :: 5 - Production/Stable',
+ 'Environment :: Console',
+ 'Intended Audience :: Science/Research',
+ 'Operating System :: OS Independent',
+ 'Programming Language :: Python :: 3',
+ 'Topic :: Scientific/Engineering :: GIS'
+]
+authors = [{name = "Michael Hirsch"}]
+
+[project.urls]
+homepage = "https://github.com/geospace-code/pymap3d"
+
+[project.optional-dependencies]
+test = [
+ "pytest"
+]
+lint = [
+ "flake8",
+ "flake8-bugbear",
+ "flake8-builtins",
+ "flake8-blind-except",
+ "mypy >= 0.800"
+]
+full = [
+ "python-dateutil",
+ "numpy >= 1.10.0",
+ "astropy",
+ "xarray"
+]
[tool.black]
line-length = 100
[tool.pytest.ini_options]
-addopts = "-ra -v"
+addopts = "-ra"
=====================================
scripts/benchmark_ecef2geo.py
=====================================
@@ -0,0 +1,31 @@
+#!/usr/bin/env python3
+"""
+benchmark ecef2geodetic
+"""
+import time
+from pymap3d.ecef import ecef2geodetic
+import numpy as np
+import argparse
+
+ll0 = (42.0, 82.0)
+
+
+def bench(N: int) -> float:
+
+ x = np.random.random(N)
+ y = np.random.random(N)
+ z = np.random.random(N)
+
+ tic = time.monotonic()
+ lat, lon, alt = ecef2geodetic(x, y, z)
+
+ return time.monotonic() - tic
+
+
+if __name__ == "__main__":
+ p = argparse.ArgumentParser()
+ p.add_argument("N", type=int)
+ p = p.parse_args()
+ N = p.N
+
+ print(f"ecef2geodetic: {bench(N):.3f} seconds")
=====================================
scripts/benchmark_vincenty.py
=====================================
@@ -31,7 +31,7 @@ def bench_vreckon(N: int) -> float:
az = np.random.random(N)
tic = time.monotonic()
- a, b = vreckon(*ll0, sr, az)
+ a, b = vreckon(ll0[0], ll0[1], sr, az)
return time.monotonic() - tic
@@ -41,7 +41,7 @@ def bench_vdist(N: int) -> float:
lon = np.random.random(N)
tic = time.monotonic()
- asr, aaz = vdist(*ll0, lat, lon)
+ asr, aaz = vdist(ll0[0], ll0[1], lat, lon)
return time.monotonic() - tic
=====================================
setup.cfg
=====================================
@@ -1,9 +1,11 @@
[metadata]
name = pymap3d
-version = 2.5.0
+version = 2.6.0
author = Michael Hirsch, Ph.D.
author_email = scivision at users.noreply.github.com
description = pure Python (no prereqs) coordinate conversions, following convention of several popular Matlab routines.
+long_description = file: README.md
+long_description_content_type = text/markdown
url = https://github.com/geospace-code/pymap3d
keywords =
coordinate conversion
=====================================
src/pymap3d/aer.py
=====================================
@@ -20,15 +20,15 @@ __all__ = ["aer2ecef", "ecef2aer", "geodetic2aer", "aer2geodetic", "eci2aer", "a
def ecef2aer(
- x: float | ndarray,
- y: float | ndarray,
- z: float | ndarray,
- lat0: float,
- lon0: float,
- h0: float,
+ x: ndarray,
+ y: ndarray,
+ z: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
compute azimuth, elevation and slant range from an Observer to a Point with ECEF coordinates.
@@ -69,15 +69,15 @@ def ecef2aer(
def geodetic2aer(
- lat: float,
- lon: float,
- h: float,
- lat0: float,
- lon0: float,
- h0: float,
+ lat: ndarray,
+ lon: ndarray,
+ h: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
gives azimuth, elevation and slant range from an Observer to a Point with geodetic coordinates.
@@ -117,15 +117,15 @@ def geodetic2aer(
def aer2geodetic(
- az: float,
- el: float,
- srange: float,
- lat0: float,
- lon0: float,
- h0: float,
+ az: ndarray,
+ el: ndarray,
+ srange: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
gives geodetic coordinates of a point with az, el, range
from an observer at lat0, lon0, h0
@@ -167,17 +167,17 @@ def aer2geodetic(
def eci2aer(
- x: float,
- y: float,
- z: float,
- lat0: float,
- lon0: float,
- h0: float,
+ x: ndarray,
+ y: ndarray,
+ z: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
t: datetime,
*,
deg: bool = True,
use_astropy: bool = True
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
takes Earth Centered Inertial x,y,z ECI coordinates of point and gives az, el, slant range from Observer
@@ -221,18 +221,18 @@ def eci2aer(
def aer2eci(
- az: float,
- el: float,
- srange: float,
- lat0: float,
- lon0: float,
- h0: float,
+ az: ndarray,
+ el: ndarray,
+ srange: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
t: datetime,
ell=None,
*,
deg: bool = True,
use_astropy: bool = True
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
gives ECI of a point from an observer at az, el, slant range
@@ -280,15 +280,15 @@ def aer2eci(
def aer2ecef(
- az: float,
- el: float,
- srange: float,
- lat0: float,
- lon0: float,
- alt0: float,
+ az: ndarray,
+ el: ndarray,
+ srange: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ alt0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
converts target azimuth, elevation, range from observer at lat0,lon0,alt0 to ECEF coordinates.
=====================================
src/pymap3d/ecef.py
=====================================
@@ -15,12 +15,10 @@ try:
arctan2 as atan2,
sqrt,
pi,
- vectorize,
)
except ImportError:
from math import radians, sin, cos, tan, atan, hypot, degrees, atan2, sqrt, pi # type: ignore
- vectorize = None
ndarray = typing.Any # type: ignore
from datetime import datetime
@@ -51,12 +49,12 @@ __all__ = [
def geodetic2ecef(
- lat: float | ndarray,
- lon: float | ndarray,
- alt: float | ndarray,
+ lat: ndarray,
+ lon: ndarray,
+ alt: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
point transformation from Geodetic of specified ellipsoid (default WGS-84) to ECEF
@@ -105,27 +103,12 @@ def geodetic2ecef(
def ecef2geodetic(
- x: float | ndarray,
- y: float | ndarray,
- z: float | ndarray,
- ell: Ellipsoid = None,
- deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
- if vectorize is not None:
- fun = vectorize(ecef2geodetic_point)
- lat, lon, alt = fun(x, y, z, ell, deg)
- return lat[()], lon[()], alt[()]
- else:
- return ecef2geodetic_point(x, y, z, ell, deg)
-
-
-def ecef2geodetic_point(
- x: float | ndarray,
- y: float | ndarray,
- z: float | ndarray,
+ x: ndarray,
+ y: ndarray,
+ z: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
convert ECEF (meters) to geodetic coordinates
@@ -155,6 +138,7 @@ def ecef2geodetic_point(
You, Rey-Jer. (2000). Transformation of Cartesian to Geodetic Coordinates without Iterations.
Journal of Surveying Engineering. doi: 10.1061/(ASCE)0733-9453
"""
+
if ell is None:
ell = Ellipsoid()
@@ -199,8 +183,14 @@ def ecef2geodetic_point(
+ z ** 2 / ell.semiminor_axis ** 2
< 1
)
- if inside:
- alt = -alt
+
+ try:
+ if inside.any(): # type: ignore
+ # avoid all false assignment bug
+ alt[inside] = -alt
+ except (TypeError, AttributeError):
+ if inside:
+ alt = -alt
if deg:
lat = degrees(lat)
@@ -255,15 +245,15 @@ def ecef2enuv(
def ecef2enu(
- x: float | ndarray,
- y: float | ndarray,
- z: float | ndarray,
- lat0: float,
- lon0: float,
- h0: float,
+ x: ndarray,
+ y: ndarray,
+ z: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
from observer to target, ECEF => ENU
@@ -302,13 +292,13 @@ def ecef2enu(
def enu2uvw(
- east: float,
- north: float,
- up: float,
- lat0: float,
- lon0: float,
+ east: ndarray,
+ north: ndarray,
+ up: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
Parameters
----------
@@ -342,8 +332,8 @@ def enu2uvw(
def uvw2enu(
- u: float, v: float, w: float, lat0: float, lon0: float, deg: bool = True
-) -> tuple[float, float, float]:
+ u: ndarray, v: ndarray, w: ndarray, lat0: ndarray, lon0: ndarray, deg: bool = True
+) -> tuple[ndarray, ndarray, ndarray]:
"""
Parameters
----------
@@ -384,7 +374,7 @@ def eci2geodetic(
*,
deg: bool = True,
use_astropy: bool = True
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
convert Earth Centered Internal ECI to geodetic coordinates
@@ -427,9 +417,9 @@ def eci2geodetic(
def geodetic2eci(
- lat: float | ndarray,
- lon: float | ndarray,
- alt: float | ndarray,
+ lat: ndarray,
+ lon: ndarray,
+ alt: ndarray,
t: datetime,
ell: Ellipsoid = None,
*,
@@ -478,15 +468,15 @@ def geodetic2eci(
def enu2ecef(
- e1: float,
- n1: float,
- u1: float,
- lat0: float,
- lon0: float,
- h0: float,
+ e1: ndarray,
+ n1: ndarray,
+ u1: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
ENU to ECEF
=====================================
src/pymap3d/enu.py
=====================================
@@ -3,11 +3,10 @@ from __future__ import annotations
import typing
try:
- from numpy import radians, sin, cos, hypot, arctan2 as atan2, degrees, pi, vectorize, ndarray
+ from numpy import asarray, radians, sin, cos, hypot, arctan2 as atan2, degrees, pi, ndarray
except ImportError:
from math import radians, sin, cos, hypot, atan2, degrees, pi # type: ignore
- vectorize = None
ndarray = typing.Any # type: ignore
from .ecef import geodetic2ecef, ecef2geodetic, enu2ecef, uvw2enu
@@ -20,8 +19,8 @@ __all__ = ["enu2aer", "aer2enu", "enu2geodetic", "geodetic2enu"]
def enu2aer(
- e: float | ndarray, n: float | ndarray, u: float | ndarray, deg: bool = True
-) -> tuple[float, float, float]:
+ e: ndarray, n: ndarray, u: ndarray, deg: bool = True
+) -> tuple[ndarray, ndarray, ndarray]:
"""
ENU to Azimuth, Elevation, Range
@@ -48,26 +47,19 @@ def enu2aer(
slant range [meters]
"""
- if vectorize is not None:
- fun = vectorize(enu2aer_point)
- az, el, rng = fun(e, n, u, deg)
- return az[()], el[()], rng[()]
- else:
- return enu2aer_point(e, n, u, deg)
-
-
-def enu2aer_point(
- e: float | ndarray, n: float | ndarray, u: float | ndarray, deg: bool = True
-) -> tuple[float, float, float]:
-
# 1 millimeter precision for singularity
- if abs(e) < 1e-3:
- e = 0.0
- if abs(n) < 1e-3:
- n = 0.0
- if abs(u) < 1e-3:
- u = 0.0
+ try:
+ e[abs(e) < 1e-3] = 0.0
+ n[abs(n) < 1e-3] = 0.0
+ u[abs(u) < 1e-3] = 0.0
+ except TypeError:
+ if abs(e) < 1e-3:
+ e = 0.0 # type: ignore
+ if abs(n) < 1e-3:
+ n = 0.0 # type: ignore
+ if abs(u) < 1e-3:
+ u = 0.0 # type: ignore
r = hypot(e, n)
slantRange = hypot(r, u)
@@ -81,18 +73,9 @@ def enu2aer_point(
return az, elev, slantRange
-def aer2enu(az: float, el: float, srange: float, deg: bool = True) -> tuple[float, float, float]:
- if vectorize is not None:
- fun = vectorize(aer2enu_point)
- e, n, u = fun(az, el, srange, deg)
- return e[()], n[()], u[()]
- else:
- return aer2enu_point(az, el, srange, deg)
-
-
-def aer2enu_point(
- az: float, el: float, srange: float, deg: bool = True
-) -> tuple[float, float, float]:
+def aer2enu(
+ az: ndarray, el: ndarray, srange: float | ndarray, deg: bool = True
+) -> tuple[ndarray, ndarray, ndarray]:
"""
Azimuth, Elevation, Slant range to target to East, North, Up
@@ -120,8 +103,12 @@ def aer2enu_point(
el = radians(el)
az = radians(az)
- if srange < 0:
- raise ValueError("Slant range [0, Infinity)")
+ try:
+ if (asarray(srange) < 0).any():
+ raise ValueError("Slant range [0, Infinity)")
+ except NameError:
+ if srange < 0:
+ raise ValueError("Slant range [0, Infinity)")
r = srange * cos(el)
@@ -129,15 +116,15 @@ def aer2enu_point(
def enu2geodetic(
- e: float,
- n: float,
- u: float,
- lat0: float,
- lon0: float,
- h0: float,
+ e: ndarray,
+ n: ndarray,
+ u: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
East, North, Up to target to geodetic coordinates
@@ -177,15 +164,15 @@ def enu2geodetic(
def geodetic2enu(
- lat: float,
- lon: float,
- h: float,
- lat0: float,
- lon0: float,
- h0: float,
+ lat: ndarray,
+ lon: ndarray,
+ h: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
Parameters
----------
=====================================
src/pymap3d/latitude.py
=====================================
@@ -7,7 +7,7 @@ from .utils import sanitize
from .rcurve import rcurve_transverse
try:
- from numpy import radians, degrees, tan, sin, exp, pi, sqrt, inf, vectorize, ndarray
+ from numpy import radians, degrees, tan, sin, exp, pi, sqrt, inf, ndarray
from numpy import arctan as atan, arcsinh as asinh, arctanh as atanh # noqa: A001
use_numpy = True
@@ -37,8 +37,8 @@ __all__ = [
def geoc2geod(
- geocentric_lat: float,
- geocentric_distance: float,
+ geocentric_lat: ndarray,
+ geocentric_distance: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
) -> float:
@@ -162,28 +162,7 @@ def geocentric2geodetic(
return degrees(geodetic_lat) if deg else geodetic_lat
-def geodetic2isometric_point(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
- geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)
-
- e = ell.eccentricity
-
- if abs(geodetic_lat - pi / 2) <= 1e-9:
- isometric_lat = inf
- elif abs(-geodetic_lat - pi / 2) <= 1e-9:
- isometric_lat = -inf
- else:
- isometric_lat = asinh(tan(geodetic_lat)) - e * atanh(e * sin(geodetic_lat))
- # same results
- # a1 = e * sin(geodetic_lat)
- # y = (1 - a1) / (1 + a1)
- # a2 = pi / 4 + geodetic_lat / 2
- # isometric_lat = log(tan(a2) * (y ** (e / 2)))
- # isometric_lat = log(tan(a2)) + e/2 * log((1-e*sin(geodetic_lat)) / (1+e*sin(geodetic_lat)))
-
- return degrees(isometric_lat) if deg else isometric_lat
-
-
-def geodetic2isometric(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def geodetic2isometric(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
"""
computes isometric latitude on an ellipsoid
@@ -212,11 +191,29 @@ def geodetic2isometric(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = T
School of Mathematical and Geospatial Sciences, RMIT University,
January 2010
"""
- if use_numpy:
- fun = vectorize(geodetic2isometric_point)
- return fun(geodetic_lat, ell, deg)[()]
- else:
- return geodetic2isometric_point(geodetic_lat, ell, deg)
+
+ geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)
+
+ e = ell.eccentricity
+
+ isometric_lat = asinh(tan(geodetic_lat)) - e * atanh(e * sin(geodetic_lat))
+ # same results
+ # a1 = e * sin(geodetic_lat)
+ # y = (1 - a1) / (1 + a1)
+ # a2 = pi / 4 + geodetic_lat / 2
+ # isometric_lat = log(tan(a2) * (y ** (e / 2)))
+ # isometric_lat = log(tan(a2)) + e/2 * log((1-e*sin(geodetic_lat)) / (1+e*sin(geodetic_lat)))
+
+ try:
+ isometric_lat[abs(geodetic_lat - pi / 2) <= 1e-9] = inf
+ isometric_lat[abs(-geodetic_lat - pi / 2) <= 1e-9] = -inf
+ except TypeError:
+ if abs(geodetic_lat - pi / 2) <= 1e-9:
+ isometric_lat = inf
+ elif abs(-geodetic_lat - pi / 2) <= 1e-9:
+ isometric_lat = -inf
+
+ return degrees(isometric_lat) if deg else isometric_lat
def isometric2geodetic(isometric_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
@@ -327,14 +324,7 @@ def geodetic2conformal(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = T
Office, Washington, DC, 1987, pp. 13-18.
"""
- if use_numpy:
- fun = vectorize(geodetic2conformal_point)
- return fun(geodetic_lat, ell, deg)[()]
- else:
- return geodetic2conformal_point(geodetic_lat, ell, deg)
-
-def geodetic2conformal_point(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)
e = ell.eccentricity
@@ -353,7 +343,7 @@ def geodetic2conformal_point(geodetic_lat: float, ell: Ellipsoid = None, deg: bo
# %% rectifying
-def geodetic2rectifying(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def geodetic2rectifying(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
"""
converts from geodetic latitude to rectifying latitude
@@ -399,7 +389,9 @@ def geodetic2rectifying(geodetic_lat: float, ell: Ellipsoid = None, deg: bool =
return degrees(rectifying_lat) if deg else rectifying_lat
-def rectifying2geodetic(rectifying_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def rectifying2geodetic(
+ rectifying_lat: ndarray, ell: Ellipsoid = None, deg: bool = True
+) -> ndarray:
"""
converts from rectifying latitude to geodetic latitude
=====================================
src/pymap3d/los.py
=====================================
@@ -4,11 +4,10 @@ from __future__ import annotations
import typing
try:
- from numpy import pi, nan, sqrt, vectorize, ndarray
+ from numpy import pi, nan, sqrt, ndarray, atleast_1d
except ImportError:
from math import pi, nan, sqrt # type: ignore
- vectorize = None
ndarray = typing.Any # type: ignore
from .aer import aer2enu
@@ -19,14 +18,14 @@ __all__ = ["lookAtSpheroid"]
def lookAtSpheroid(
- lat0: float,
- lon0: float,
- h0: float,
- az: float,
- tilt: float,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
+ az: ndarray,
+ tilt: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
Calculates line-of-sight intersection with Earth (or other ellipsoid) surface from above surface / orbit
@@ -62,23 +61,6 @@ def lookAtSpheroid(
Algorithm based on https://medium.com/@stephenhartzell/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6 Stephen Hartzell
"""
- if vectorize is not None:
- fun = vectorize(lookAtSpheroid_point)
- lat, lon, d = fun(lat0, lon0, h0, az, tilt, ell, deg)
- return lat[()], lon[()], d[()]
- else:
- return lookAtSpheroid_point(lat0, lon0, h0, az, tilt, ell, deg)
-
-
-def lookAtSpheroid_point(
- lat0: float,
- lon0: float,
- h0: float,
- az: float,
- tilt: float,
- ell: Ellipsoid = None,
- deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
if h0 < 0:
raise ValueError("Intersection calculation requires altitude [0, Infinity)")
@@ -86,6 +68,13 @@ def lookAtSpheroid_point(
if ell is None:
ell = Ellipsoid()
+ try:
+ lat0 = atleast_1d(lat0)
+ lon0 = atleast_1d(lon0)
+ tilt = atleast_1d(tilt)
+ except NameError:
+ pass
+
a = ell.semimajor_axis
b = ell.semimajor_axis
c = ell.semiminor_axis
@@ -115,14 +104,22 @@ def lookAtSpheroid_point(
magnitude = a ** 2 * b ** 2 * w ** 2 + a ** 2 * c ** 2 * v ** 2 + b ** 2 * c ** 2 * u ** 2
# %% Return nan if radical < 0 or d < 0 because LOS vector does not point towards Earth
- if radical > 0:
+ try:
d = (value - a * b * c * sqrt(radical)) / magnitude
- else:
- d = nan
+ d[radical < 0] = nan
+ d[d < 0] = nan
+ except ValueError:
+ if radical < 0:
+ d = nan
+ if d < 0:
+ d = nan
+ except TypeError:
+ pass
- if d < 0:
- d = nan
# %% cartesian to ellipsodal
lat, lon, _ = ecef2geodetic(x + d * u, y + d * v, z + d * w, deg=deg)
- return lat, lon, d
+ try:
+ return lat.squeeze()[()], lon.squeeze()[()], d.squeeze()[()]
+ except AttributeError:
+ return lat, lon, d
=====================================
src/pymap3d/lox.py
=====================================
@@ -4,11 +4,10 @@ from __future__ import annotations
import typing
try:
- from numpy import radians, degrees, cos, arctan2 as atan2, tan, pi, vectorize, ndarray
+ from numpy import radians, degrees, cos, arctan2 as atan2, tan, pi, ndarray, atleast_1d
except ImportError:
from math import radians, degrees, cos, atan2, tan, pi # type: ignore
- vectorize = None
ndarray = typing.Any # type: ignore
import typing
@@ -34,7 +33,7 @@ __all__ = [
]
-def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def meridian_dist(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
"""
Computes the ground distance on an ellipsoid from the equator to the input latitude.
@@ -52,10 +51,10 @@ def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
dist : float
distance (meters)
"""
- return meridian_arc(0, lat, ell, deg)
+ return meridian_arc(0.0, lat, ell, deg)
-def meridian_arc(lat1: float, lat2: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def meridian_arc(lat1, lat2: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
"""
Computes the ground distance on an ellipsoid between two latitudes.
@@ -84,10 +83,10 @@ def meridian_arc(lat1: float, lat2: float, ell: Ellipsoid = None, deg: bool = Tr
def loxodrome_inverse(
- lat1: float,
- lon1: float,
- lat2: float,
- lon2: float,
+ lat1: ndarray,
+ lon1: ndarray,
+ lat2: ndarray,
+ lon2: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
) -> tuple[float, float]:
@@ -135,17 +134,7 @@ def loxodrome_inverse(
Survey, U.S. Department of Commerce, Washington, DC: U.S.
Government Printing Office, p. 66.
"""
- if vectorize is not None:
- fun = vectorize(loxodrome_inverse_point)
- lox, az = fun(lat1, lon1, lat2, lon2, ell, deg)
- return lox[()], az[()]
- else:
- return loxodrome_inverse_point(lat1, lon1, lat2, lon2, ell, deg)
-
-def loxodrome_inverse_point(
- lat1: float, lon1: float, lat2: float, lon2: float, ell: Ellipsoid = None, deg: bool = True
-) -> typing.Tuple[float, float]:
if deg:
lat1, lon1, lat2, lon2 = radians(lat1), radians(lon1), radians(lat2), radians(lon2)
@@ -160,10 +149,13 @@ def loxodrome_inverse_point(
cosaz12 = cos(az12)
# compute distance along loxodromic curve
- if abs(cosaz12) < 1e-9: # straight east or west
- dist = departure(lon2, lon1, lat1, ell, deg=False)
- else:
- dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / abs(cos(az12))
+ dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / abs(cos(az12))
+ try:
+ if (abs(cosaz12) < 1e-9).any():
+ dist[abs(cosaz12) < 1e-9] = departure(lon2, lon1, lat1, ell, deg=False)
+ except (AttributeError, TypeError):
+ if abs(cosaz12) < 1e-9: # straight east or west
+ dist = departure(lon2, lon1, lat1, ell, deg=False)
if deg:
az12 = degrees(az12) % 360.0
@@ -172,13 +164,13 @@ def loxodrome_inverse_point(
def loxodrome_direct(
- lat1: float,
- lon1: float,
- rng: float,
+ lat1: ndarray,
+ lon1: ndarray,
+ rng: ndarray,
a12: float,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float, float]:
+) -> tuple[ndarray, ndarray]:
"""
Given starting lat, lon with arclength and azimuth, compute final lat, lon
@@ -206,25 +198,22 @@ def loxodrome_direct(
lon2 : float
final geodetic longitude (degrees)
"""
- if vectorize is not None:
- fun = vectorize(loxodrome_direct_point)
- lat2, lon2 = fun(lat1, lon1, rng, a12, ell, deg)
- return lat2[()], lon2[()]
- else:
- return loxodrome_direct_point(lat1, lon1, rng, a12, ell, deg)
-
-
-def loxodrome_direct_point(
- lat1: float, lon1: float, rng: float, a12: float, ell: Ellipsoid = None, deg: bool = True
-) -> typing.Tuple[float, float]:
if deg:
lat1, lon1, a12 = radians(lat1), radians(lon1), radians(a12)
- if abs(lat1) > pi / 2:
- raise ValueError("-90 <= latitude <= 90")
- if rng < 0:
- raise ValueError("ground distance must be >= 0")
+ try:
+ lat1 = atleast_1d(lat1)
+ rng = atleast_1d(rng)
+ if (abs(lat1) > pi / 2).any():
+ raise ValueError("-90 <= latitude <= 90")
+ if (rng < 0).any():
+ raise ValueError("ground distance must be >= 0")
+ except NameError:
+ if abs(lat1) > pi / 2:
+ raise ValueError("-90 <= latitude <= 90")
+ if rng < 0:
+ raise ValueError("ground distance must be >= 0")
# compute rectifying sphere latitude and radius
reclat = geodetic2rectifying(lat1, ell, deg=False)
@@ -243,11 +232,14 @@ def loxodrome_direct_point(
if deg:
lat2, lon2 = degrees(lat2), degrees(lon2)
- return lat2, lon2
+ try:
+ return lat2.squeeze()[()], lon2.squeeze()[()]
+ except AttributeError:
+ return lat2, lon2
def departure(
- lon1: float, lon2: float, lat: float, ell: Ellipsoid = None, deg: bool = True
+ lon1: ndarray, lon2: ndarray, lat: ndarray, ell: Ellipsoid = None, deg: bool = True
) -> float:
"""
Computes the distance along a specific parallel between two meridians.
=====================================
src/pymap3d/ned.py
=====================================
@@ -14,8 +14,8 @@ from .ellipsoid import Ellipsoid
def aer2ned(
- az: float, elev: float, slantRange: float, deg: bool = True
-) -> tuple[float, float, float]:
+ az: ndarray, elev: ndarray, slantRange: ndarray, deg: bool = True
+) -> tuple[ndarray, ndarray, ndarray]:
"""
converts azimuth, elevation, range to target from observer to North, East, Down
@@ -45,7 +45,9 @@ def aer2ned(
return n, e, -u
-def ned2aer(n: float, e: float, d: float, deg: bool = True) -> tuple[float, float, float]:
+def ned2aer(
+ n: ndarray, e: ndarray, d: ndarray, deg: bool = True
+) -> tuple[ndarray, ndarray, ndarray]:
"""
converts North, East, Down to azimuth, elevation, range
@@ -75,15 +77,15 @@ def ned2aer(n: float, e: float, d: float, deg: bool = True) -> tuple[float, floa
def ned2geodetic(
- n: float,
- e: float,
- d: float,
- lat0: float,
- lon0: float,
- h0: float,
+ n: ndarray,
+ e: ndarray,
+ d: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
Converts North, East, Down to target latitude, longitude, altitude
@@ -124,15 +126,15 @@ def ned2geodetic(
def ned2ecef(
- n: float,
- e: float,
- d: float,
- lat0: float,
- lon0: float,
- h0: float,
+ n: ndarray,
+ e: ndarray,
+ d: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float, float, float]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
North, East, Down to target ECEF coordinates
@@ -170,15 +172,15 @@ def ned2ecef(
def ecef2ned(
- x: float,
- y: float,
- z: float,
- lat0: float,
- lon0: float,
- h0: float,
+ x: ndarray,
+ y: ndarray,
+ z: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
Convert ECEF x,y,z to North, East, Down
@@ -219,15 +221,15 @@ def ecef2ned(
def geodetic2ned(
- lat: float,
- lon: float,
- h: float,
- lat0: float,
- lon0: float,
- h0: float,
+ lat: ndarray,
+ lon: ndarray,
+ h: ndarray,
+ lat0: ndarray,
+ lon0: ndarray,
+ h0: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray, ndarray]:
"""
convert latitude, longitude, altitude of target to North, East, Down from observer
=====================================
src/pymap3d/rcurve.py
=====================================
@@ -36,7 +36,7 @@ def geocentric_radius(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = Tr
)
-def rcurve_parallel(lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def rcurve_parallel(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
"""
computes the radius of the small circle encompassing the globe at the specified latitude
=====================================
src/pymap3d/rsphere.py
=====================================
@@ -94,13 +94,13 @@ def rsphere_rectifying(ell: Ellipsoid = None) -> float:
def rsphere_euler(
- lat1: float | ndarray,
- lon1: float,
- lat2: float | ndarray,
- lon2: float,
+ lat1: ndarray,
+ lon1: ndarray,
+ lat2: ndarray,
+ lon2: ndarray,
ell: Ellipsoid = None,
deg: bool = True,
-) -> float:
+) -> ndarray:
"""computes the Euler radii of curvature at the midpoint of the
great circle arc defined by the endpoints (lat1,lon1) and (lat2,lon2)
=====================================
src/pymap3d/tests/__init__.py
=====================================
=====================================
src/pymap3d/tests/test_aer.py
=====================================
@@ -16,6 +16,9 @@ def test_aer2ecef(aer, lla, xyz):
assert x == approx(xyz[0])
assert y == approx(xyz[1])
assert z == approx(xyz[2])
+ assert isinstance(x, float)
+ assert isinstance(y, float)
+ assert isinstance(z, float)
raer = (radians(aer[0]), radians(aer[1]), aer[2])
rlla = (radians(lla[0]), radians(lla[1]), lla[2])
=====================================
src/pymap3d/tests/test_geodetic.py
=====================================
@@ -37,17 +37,11 @@ def test_3d_geodetic2ecef():
assert (x0, y0, z0) == approx(xyz0)
- at pytest.mark.parametrize(
- "xyz", [(xyz0[0], xyz0[1], xyz0[2]), ([xyz0[0]], [xyz0[1]], [xyz0[2]])], ids=("scalar", "list")
-)
-def test_scalar_ecef2geodetic(xyz):
+def test_scalar_ecef2geodetic():
"""
verify we can handle the wide variety of input data type users might use
"""
- if isinstance(xyz[0], list):
- pytest.importorskip("numpy")
-
- lat, lon, alt = pm.ecef2geodetic(*xyz)
+ lat, lon, alt = pm.ecef2geodetic(xyz0[0], xyz0[1], xyz0[2])
assert [lat, lon, alt] == approx(lla0, rel=1e-4)
=====================================
src/pymap3d/tests/test_vincenty_vreckon.py
=====================================
@@ -30,7 +30,10 @@ def test_unit(lat, lon, srange, az, lato, lono):
lat1, lon1 = vincenty.vreckon(lat, lon, srange, az)
assert lat1 == approx(lato)
+ assert isinstance(lat1, float)
+
assert lon1 == approx(lono, rel=0.001)
+ assert isinstance(lon1, float)
def test_az_vector():
=====================================
src/pymap3d/utils.py
=====================================
@@ -8,26 +8,26 @@ import typing
from .ellipsoid import Ellipsoid
try:
- from numpy import hypot, cos, sin, arctan2 as atan2, radians, pi, asarray, ndarray
+ from numpy import hypot, cos, sin, arctan2 as atan2, radians, pi, asarray, ndarray, sign
except ImportError:
from math import atan2, hypot, cos, sin, radians, pi # type: ignore
asarray = None
ndarray = typing.Any # type: ignore
-__all__ = ["cart2pol", "pol2cart", "cart2sph", "sph2cart", "sign"]
+ def sign(x: float) -> float: # type: ignore
+ """ signum function """
+ if x < 0:
+ y = -1.0
+ elif x > 0:
+ y = 1.0
+ else:
+ y = 0.0
+ return y
-def sign(x: float) -> float:
- """ signum function """
- if x < 0:
- y = -1.0
- elif x > 0:
- y = 1.0
- else:
- y = 0.0
- return y
+__all__ = ["cart2pol", "pol2cart", "cart2sph", "sph2cart", "sign"]
def cart2pol(x: float, y: float) -> tuple[float, float]:
=====================================
src/pymap3d/vincenty.py
=====================================
@@ -5,13 +5,27 @@ Vincenty's methods for computing ground distance and reckoning
from __future__ import annotations
import typing
import logging
-from math import sqrt, tan, sin, cos, isnan, atan, atan2, asin, nan, pi, radians, degrees
+from math import nan, pi
from copy import copy
try:
- from numpy import ndarray, vectorize
+ from numpy import (
+ ndarray,
+ atleast_1d,
+ sqrt,
+ tan,
+ sin,
+ cos,
+ isnan,
+ arctan as atan,
+ arctan2 as atan2,
+ arcsin as asin,
+ radians,
+ degrees,
+ )
except ImportError:
- vectorize = None
+ from math import sqrt, tan, sin, cos, isnan, atan, atan2, asin, radians, degrees # type: ignore
+
ndarray = typing.Any # type: ignore
from .ellipsoid import Ellipsoid
@@ -21,8 +35,12 @@ __all__ = ["vdist", "vreckon", "track2"]
def vdist(
- Lat1: float | ndarray, Lon1: float, Lat2: float | ndarray, Lon2: float, ell: Ellipsoid = None
-) -> tuple[float, float]:
+ Lat1: float | ndarray,
+ Lon1: float | ndarray,
+ Lat2: float | ndarray,
+ Lon2: float | ndarray,
+ ell: Ellipsoid = None,
+) -> tuple[ndarray, ndarray]:
"""
Using the reference ellipsoid, compute the distance between two points
within a few millimeters of accuracy, compute forward azimuth,
@@ -93,22 +111,20 @@ def vdist(
11. Azimuths agree with the Mapping Toolbox, version 2.2 (R14SP3) to within about a hundred-thousandth of a degree, except when traversing to or from a pole, where the convention for this function is described in (10), and except in the cases noted above in (9).
12. No warranties; use at your own risk.
"""
- if vectorize is not None:
- fun = vectorize(vdist_point)
- dist, az = fun(Lat1, Lon1, Lat2, Lon2, ell)
- return dist[()], az[()]
- else:
- return vdist_point(Lat1, Lon1, Lat2, Lon2, ell)
-
-def vdist_point(
- Lat1: float | ndarray, Lon1: float, Lat2: float | ndarray, Lon2: float, ell: Ellipsoid = None
-) -> tuple[float, float]:
if ell is None:
ell = Ellipsoid()
# %% Input check:
- if (abs(Lat1) > 90) | (abs(Lat2) > 90):
- raise ValueError("Input latitudes must be in [-90, 90] degrees.")
+ try:
+ Lat1 = atleast_1d(Lat1)
+ Lon1 = atleast_1d(Lon1)
+ Lat2 = atleast_1d(Lat2)
+ Lon2 = atleast_1d(Lon2)
+ if (abs(Lat1) > 90).any() | (abs(Lat2) > 90).any():
+ raise ValueError("Input latitudes must be in [-90, 90] degrees.")
+ except NameError:
+ if (abs(Lat1) > 90) | (abs(Lat2) > 90):
+ raise ValueError("Input latitudes must be in [-90, 90] degrees.")
# %% Supply WGS84 earth ellipsoid axis lengths in meters:
a = ell.semimajor_axis
b = ell.semiminor_axis
@@ -119,19 +135,30 @@ def vdist_point(
lat2 = radians(Lat2)
lon2 = radians(Lon2)
# %% correct for errors at exact poles by adjusting 0.6 millimeters:
- if abs(pi / 2 - abs(lat1)) < 1e-10:
- lat1 = sign(lat1) * (pi / 2 - 1e-10)
+ try:
+ i = abs(pi / 2 - abs(lat1)) < 1e-10
+ lat1[i] = sign(lat1[i]) * (pi / 2 - 1e-10)
- if abs(pi / 2 - abs(lat2)) < 1e-10:
- lat2 = sign(lat2) * (pi / 2 - 1e-10)
+ i = abs(pi / 2 - abs(lat2)) < 1e-10
+ lat2[i] = sign(lat2[i]) * (pi / 2 - 1e-10)
+ except TypeError:
+ if abs(pi / 2 - abs(lat1)) < 1e-10:
+ lat1 = sign(lat1) * (pi / 2 - 1e-10)
+
+ if abs(pi / 2 - abs(lat2)) < 1e-10:
+ lat2 = sign(lat2) * (pi / 2 - 1e-10)
U1 = atan((1 - f) * tan(lat1))
U2 = atan((1 - f) * tan(lat2))
lon1 = lon1 % (2 * pi)
lon2 = lon2 % (2 * pi)
L = abs(lon2 - lon1)
- if L > pi:
- L = 2 * pi - L
+
+ try:
+ L[L > pi] = 2 * pi - L[L > pi]
+ except TypeError:
+ if L > pi:
+ L = 2 * pi - L
lamb = copy(L) # NOTE: program will fail without copy!
itercount = 0
@@ -161,15 +188,23 @@ def vdist_point(
try:
sinAlpha = cos(U1) * cos(U2) * sin(lamb) / sin(sigma)
- except ZeroDivisionError:
- sinAlpha = 0.0
-
- if isnan(sinAlpha):
- alpha = 0.0
- elif sinAlpha > 1 or abs(sinAlpha - 1) < 1e-16:
- alpha = pi / 2
- else:
+
alpha = asin(sinAlpha)
+ alpha[isnan(sinAlpha)] = 0
+ alpha[(sinAlpha > 1) | (abs(sinAlpha - 1) < 1e-16)] = pi / 2
+
+ except (ZeroDivisionError, TypeError, ValueError):
+ try:
+ sinAlpha = cos(U1) * cos(U2) * sin(lamb) / sin(sigma)
+ except ZeroDivisionError:
+ sinAlpha = 0.0
+
+ if isnan(sinAlpha):
+ alpha = 0.0
+ elif sinAlpha > 1 or abs(sinAlpha - 1) < 1e-16:
+ alpha = pi / 2
+ else:
+ alpha = asin(sinAlpha)
cos2sigmam = cos(sigma) - 2 * sin(U1) * sin(U2) / cos(alpha) ** 2
@@ -180,7 +215,12 @@ def vdist_point(
)
# print(f'then, lambda(21752) = {lamb[21752],20})
# correct for convergence failure for essentially antipodal points
- if lamb > pi:
+ try:
+ i = (lamb > pi).any()
+ except AttributeError:
+ i = lamb > pi
+
+ if i:
logging.warning(
"Essentially antipodal points encountered. Precision may be reduced slightly."
)
@@ -188,7 +228,10 @@ def vdist_point(
lambdaold = pi
lamb = pi
- notdone = abs(lamb - lambdaold) > 1e-12
+ try:
+ notdone = (abs(lamb - lambdaold) > 1e-12).any()
+ except AttributeError:
+ notdone = abs(lamb - lambdaold) > 1e-12
u2 = cos(alpha) ** 2 * (a ** 2 - b ** 2) / b ** 2
A = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)))
@@ -212,8 +255,14 @@ def vdist_point(
# %% From point #1 to point #2
# correct sign of lambda for azimuth calcs:
lamb = abs(lamb)
- if sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0:
- lamb = -lamb
+
+ try:
+ i = sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0
+ lamb[i] = -lamb[i]
+ except TypeError:
+ if sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0:
+ lamb = -lamb
+
numer = cos(U2) * sin(lamb)
denom = cos(U1) * sin(U2) - sin(U1) * cos(U2) * cos(lamb)
a12 = atan2(numer, denom)
@@ -221,12 +270,15 @@ def vdist_point(
az = degrees(a12)
- return dist_m, az
+ try:
+ return dist_m.squeeze()[()], az.squeeze()[()]
+ except AttributeError:
+ return dist_m, az
def vreckon(
- Lat1: float, Lon1: float, Rng: float, Azim: float, ell: Ellipsoid = None
-) -> tuple[float, float]:
+ Lat1: float | ndarray, Lon1: float | ndarray, Rng: ndarray, Azim: ndarray, ell: Ellipsoid = None
+) -> tuple[ndarray, ndarray]:
"""
This is the Vincenty "forward" solution.
@@ -283,20 +335,21 @@ def vreckon(
Added ellipsoid and vectorized whenever possible. Also, lon2 is always converted to the [-180 180] interval.
Joaquim Luis
"""
- if vectorize is not None:
- fun = vectorize(vreckon_point)
- return fun(Lat1, Lon1, Rng, Azim, ell)
- else:
- return vreckon_point(Lat1, Lon1, Rng, Azim, ell)
-
-def vreckon_point(
- Lat1: float, Lon1: float, Rng: float, Azim: float, ell: Ellipsoid = None
-) -> tuple[float, float]:
- if abs(Lat1) > 90:
- raise ValueError("Input lat. must be between -90 and 90 deg., inclusive.")
- if Rng < 0:
- raise ValueError("Ground distance must be positive")
+ try:
+ Lat1 = atleast_1d(Lat1)
+ Lon1 = atleast_1d(Lon1)
+ Rng = atleast_1d(Rng)
+ Azim = atleast_1d(Azim)
+ if (abs(Lat1) > 90).any():
+ raise ValueError("Input lat. must be between -90 and 90 deg., inclusive.")
+ if (Rng < 0).any():
+ raise ValueError("Ground distance must be positive")
+ except NameError:
+ if abs(Lat1) > 90:
+ raise ValueError("Input lat. must be between -90 and 90 deg., inclusive.")
+ if Rng < 0:
+ raise ValueError("Ground distance must be positive")
if ell is not None:
a = ell.semimajor_axis
@@ -311,8 +364,12 @@ def vreckon_point(
lon1 = radians(Lon1) # intial longitude in radians
# correct for errors at exact poles by adjusting 0.6 millimeters:
- if abs(pi / 2 - abs(lat1)) < 1e-10:
- lat1 = sign(lat1) * (pi / 2 - (1e-10))
+ try:
+ i = abs(pi / 2 - abs(lat1)) < 1e-10
+ lat1[i] = sign(lat1[i]) * (pi / 2 - (1e-10))
+ except TypeError:
+ if abs(pi / 2 - abs(lat1)) < 1e-10:
+ lat1 = sign(lat1) * (pi / 2 - (1e-10))
alpha1 = radians(Azim) # inital azimuth in radians
sinAlpha1 = sin(alpha1)
@@ -334,7 +391,13 @@ def vreckon_point(
sinSigma = nan
cosSigma = nan
cos2SigmaM = nan
- while abs(sigma - sigmaP) > 1e-12:
+
+ try:
+ i = (abs(sigma - sigmaP) > 1e-12).any()
+ except AttributeError:
+ i = abs(sigma - sigmaP) > 1e-12
+
+ while i:
cos2SigmaM = cos(2 * sigma1 + sigma)
sinSigma = sin(sigma)
cosSigma = cos(sigma)
@@ -357,6 +420,10 @@ def vreckon_point(
)
sigmaP = sigma
sigma = Rng / (b * A) + deltaSigma
+ try:
+ i = (abs(sigma - sigmaP) > 1e-12).any()
+ except AttributeError:
+ i = abs(sigma - sigmaP) > 1e-12
tmp = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1
lat2 = atan2(
@@ -383,18 +450,21 @@ def vreckon_point(
lon2 = lon2 % 360 # follow [0, 360) convention
- return degrees(lat2), lon2
+ try:
+ return degrees(lat2).squeeze()[()], lon2.squeeze()[()]
+ except AttributeError:
+ return degrees(lat2), lon2
def track2(
- lat1: float,
- lon1: float,
- lat2: float,
- lon2: float,
+ lat1: ndarray,
+ lon1: ndarray,
+ lat2: ndarray,
+ lon2: ndarray,
ell: Ellipsoid = None,
npts: int = 100,
deg: bool = True,
-) -> tuple[list[float], list[float]]:
+) -> tuple[list[ndarray], list[ndarray]]:
"""
computes great circle tracks starting at the point lat1, lon1 and ending at lat2, lon2
View it on GitLab: https://salsa.debian.org/debian-gis-team/pymap3d/-/commit/486084e0fba575d50c2a1c56e23acbde439e2a5a
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pymap3d/-/commit/486084e0fba575d50c2a1c56e23acbde439e2a5a
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20210321/8017ec90/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list