[Git][debian-gis-team/pymap3d][upstream] New upstream version 2.7.2
Antonio Valentino (@antonio.valentino)
gitlab at salsa.debian.org
Wed Oct 20 14:03:33 BST 2021
Antonio Valentino pushed to branch upstream at Debian GIS Project / pymap3d
Commits:
6e091f6a by Antonio Valentino at 2021-10-20T12:14:38+00:00
New upstream version 2.7.2
- - - - -
30 changed files:
- .github/workflows/ci.yml
- .github/workflows/ci_stdlib_only.yml
- − .github/workflows/codeql-analysis.yml
- Examples/angle_distance.py
- + Examples/has_map_toolbox.m
- + Examples/lox_stability.py
- Examples/vdist_poi.py
- Examples/vdist_stability.py
- README.md
- pyproject.toml
- setup.cfg
- src/pymap3d/aer.py
- src/pymap3d/ecef.py
- src/pymap3d/eci.py
- src/pymap3d/ellipsoid.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/sidereal.py
- src/pymap3d/tests/test_latitude.py
- src/pymap3d/tests/test_rcurve.py
- src/pymap3d/tests/test_rhumb.py
- src/pymap3d/timeconv.py
- src/pymap3d/utils.py
- src/pymap3d/vincenty.py
- src/pymap3d/vreckon/__main__.py
Changes:
=====================================
.github/workflows/ci.yml
=====================================
@@ -5,7 +5,9 @@ on:
paths:
- "**.py"
- .github/workflows/ci.yml
-
+ pull-request:
+ - "**.py"
+ - .github/workflows/ci_stdlib_only.yml
jobs:
=====================================
.github/workflows/ci_stdlib_only.yml
=====================================
@@ -5,17 +5,25 @@ on:
paths:
- "**.py"
- .github/workflows/ci_stdlib_only.yml
-
+ pull-request:
+ - "**.py"
+ - .github/workflows/ci_stdlib_only.yml
jobs:
stdlib_only:
runs-on: ubuntu-latest
+
+ name: Python ${{ matrix.python-version }}
+ strategy:
+ matrix:
+ python-version: [ '3.7', '3.10' ]
+
steps:
- uses: actions/checkout at v2
- uses: actions/setup-python at v2
with:
- python-version: 3.7
+ python-version: ${{ matrix.python-version }}
- run: pip install .[tests]
=====================================
.github/workflows/codeql-analysis.yml deleted
=====================================
@@ -1,31 +0,0 @@
-name: "CodeQL"
-
-on:
- schedule:
- - cron: '0 10 * * 6'
-
-jobs:
- analyse:
- name: Analyse
- runs-on: ubuntu-latest
-
- steps:
- - name: Checkout repository
- uses: actions/checkout at v2
- with:
- # We must fetch at least the immediate parents so that if this is
- # a pull request then we can checkout the head.
- fetch-depth: 2
-
- # If this run was triggered by a pull request event, then checkout
- # the head of the pull request instead of the merge commit.
- - run: git checkout HEAD^2
- if: ${{ github.event_name == 'pull_request' }}
-
- - name: Initialize CodeQL
- uses: github/codeql-action/init at v1
- with:
- languages: python
-
- - name: Perform CodeQL Analysis
- uses: github/codeql-action/analyze at v1
=====================================
Examples/angle_distance.py
=====================================
@@ -15,7 +15,7 @@ def main():
dist_deg = anglesep_meeus(a.r0, a.d0, a.r1, a.d1)
dist_deg_astropy = anglesep(a.r0, a.d0, a.r1, a.d1)
- print("{:.6f} deg sep".format(dist_deg))
+ print(f"{dist_deg:.6f} deg sep")
assert dist_deg == approx(dist_deg_astropy)
=====================================
Examples/has_map_toolbox.m
=====================================
@@ -0,0 +1,3 @@
+function has_map = has_map_toolbox()
+ has_map = ~isempty(ver("map"));
+end
=====================================
Examples/lox_stability.py
=====================================
@@ -0,0 +1,38 @@
+#!/usr/bin/env python3
+from __future__ import annotations
+import logging
+from pathlib import Path
+from math import isclose
+
+from pymap3d.lox import loxodrome_direct
+
+import matlab.engine
+
+cwd = Path(__file__).parent
+eng = None # don't start Matlab engine over and over when script is interactive
+
+if eng is None:
+ eng = matlab.engine.start_matlab("-nojvm")
+ eng.addpath(eng.genpath(str(cwd)), nargout=0)
+
+if not eng.has_map_toolbox():
+ raise EnvironmentError("Matlab does not have Mapping Toolbox")
+
+
+def matlab_func(lat1: float, lon1: float, rng: float, az: float) -> tuple[float, float]:
+ """Using Matlab Engine to do same thing as Pymap3d"""
+ return eng.reckon("rh", lat1, lon1, rng, az, eng.wgs84Ellipsoid(), nargout=2) # type: ignore
+
+
+clat, clon, rng = 35.0, 140.0, 50000.0 # arbitrary
+
+for i in range(20):
+ for azi in (90 + 10.0 ** (-i), -90 + 10.0 ** (-i), 270 + 10.0 ** (-i), -270 + 10.0 ** (-i)):
+ lat, lon = loxodrome_direct(clat, clon, rng, azi)
+
+ lat_matlab, lon_matlab = matlab_func(clat, clon, rng, azi)
+ rstr = f"azimuth: {azi} lat,lon: Python: {lat}, {lon} Matlab: {lat_matlab}, {lon_matlab}"
+ if not (
+ isclose(lat_matlab, lat, rel_tol=0.005) and isclose(lon_matlab, lon, rel_tol=0.001)
+ ):
+ logging.error(rstr)
=====================================
Examples/vdist_poi.py
=====================================
@@ -32,13 +32,13 @@ def get_place_coords(
keyfn = Path(keyfn).expanduser()
key = keyfn.read_text()
- stub = URL + "location={},{}".format(latitude, longitude)
+ stub = URL + f"location={latitude},{longitude}"
- stub += "&radius={}".format(search_radius_km * 1000)
+ stub += f"&radius={search_radius_km * 1000}"
- stub += "&types={}".format(place_type)
+ stub += f"&types={place_type}"
- stub += "&key={}".format(key)
+ stub += f"&key={key}"
r = requests.get(stub)
r.raise_for_status()
=====================================
Examples/vdist_stability.py
=====================================
@@ -1,24 +1,26 @@
-#!/usr/bin/env python
+#!/usr/bin/env python3
+from __future__ import annotations
from pymap3d.vincenty import vdist
import sys
+from pathlib import Path
from math import isclose, nan
-import typing
-eng = None # don't start engine over and over when script is interactive
-try:
- import matlab.engine
+import matlab.engine
- if eng is None:
- eng = matlab.engine.start_matlab("-nojvm")
-except ImportError as exc:
- print(exc, file=sys.stderr)
- eng = None
+cwd = Path(__file__).parent
+eng = None # don't start Matlab engine over and over when script is interactive
+if eng is None:
+ eng = matlab.engine.start_matlab("-nojvm")
+ eng.addpath(eng.genpath(str(cwd)), nargout=0)
-def matlab_func(lat1: float, lon1: float, lat2: float, lon2: float) -> typing.Tuple[float, float]:
+if not eng.has_map_toolbox():
+ raise EnvironmentError("Matlab does not have Mapping Toolbox")
+
+
+def matlab_func(lat1, lon1, lat2, lon2) -> tuple[float, float]:
"""Using Matlab Engine to do same thing as Pymap3d"""
- ell = eng.wgs84Ellipsoid()
- return eng.distance(lat1, lon1, lat2, lon2, ell, nargout=2)
+ return eng.distance(lat1, lon1, lat2, lon2, eng.wgs84Ellipsoid(), nargout=2) # type: ignore
dlast, alast = nan, nan
=====================================
README.md
=====================================
@@ -5,8 +5,6 @@
[](https://lgtm.com/projects/g/geospace-code/pymap3d/context:python)
-
-[](https://codecov.io/gh/geospace-code/pymap3d)
[](https://pypi.python.org/pypi/pymap3d)
[](http://pepy.tech/project/pymap3d)
@@ -16,12 +14,12 @@ PyMap3D is intended for non-interactive use on massively parallel (HPC) and embe
[API docs](https://geospace-code.github.io/pymap3d/)
-Thanks to our [contributors](./contributors.md).
+Thanks to our [contributors](./.github/contributors.md).
## Similar toolboxes in other code languages
* [Matlab, GNU Octave](https://github.com/geospace-code/matmap3d)
-* [Fortran](https://github.com/geospace-code/maptran)
+* [Fortran](https://github.com/geospace-code/maptran3d)
* [Rust](https://github.com/gberrante/map_3d)
## Prerequisites
@@ -46,7 +44,7 @@ pip install -e pymap3d
One can verify Python functionality after installation by:
```sh
-pytest pymap3d -r a -v
+pytest pymap3d
```
## Usage
@@ -63,7 +61,7 @@ az,el,range = pm.geodetic2aer(lat, lon, alt, observer_lat, observer_lon, 0)
```
[Python](https://www.python.org/dev/peps/pep-0448/)
-[argument unpacking](https://docs.python.org/3.6/tutorial/controlflow.html#unpacking-argument-lists)
+[argument unpacking](https://docs.python.org/3/tutorial/controlflow.html#unpacking-argument-lists)
can be used for compact function arguments with scalars or arbitrarily
shaped N-D arrays:
@@ -76,7 +74,7 @@ lla = pm.aer2geodetic(*aer,*obslla)
where tuple `lla` is comprised of scalar or N-D arrays `(lat,lon,alt)`.
-Example scripts are in the [examples](./examples) directory.
+Example scripts are in the [examples](./Examples) directory.
### Functions
=====================================
pyproject.toml
=====================================
@@ -39,6 +39,3 @@ full = [
[tool.black]
line-length = 100
-
-[tool.pytest.ini_options]
-addopts = "-ra"
=====================================
setup.cfg
=====================================
@@ -1,6 +1,6 @@
[metadata]
name = pymap3d
-version = 2.7.0
+version = 2.7.2
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.
@@ -38,6 +38,8 @@ lint =
flake8-builtins
flake8-blind-except
mypy >= 0.800
+ types-python-dateutil
+ types-requests
full =
python-dateutil
numpy >= 1.10.0
=====================================
src/pymap3d/aer.py
=====================================
@@ -10,11 +10,11 @@ from .ellipsoid import Ellipsoid
try:
from .eci import eci2ecef, ecef2eci
- from numpy import ndarray
except ImportError:
- eci2ecef = ecef2eci = None
- ndarray = typing.Any # type: ignore
+ eci2ecef = ecef2eci = None # type: ignore
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
__all__ = ["aer2ecef", "ecef2aer", "geodetic2aer", "aer2geodetic", "eci2aer", "aer2eci"]
=====================================
src/pymap3d/ecef.py
=====================================
@@ -3,24 +3,11 @@ from __future__ import annotations
import typing
try:
- from numpy import (
- ndarray,
- radians,
- sin,
- cos,
- tan,
- arctan as atan,
- hypot,
- degrees,
- arctan2 as atan2,
- sqrt,
- pi,
- )
+ from numpy import radians, sin, cos, tan, arctan as atan, hypot, degrees, arctan2 as atan2, sqrt
except ImportError:
- from math import radians, sin, cos, tan, atan, hypot, degrees, atan2, sqrt, pi # type: ignore
-
- ndarray = typing.Any # type: ignore
+ from math import radians, sin, cos, tan, atan, hypot, degrees, atan2, sqrt # type: ignore
+from math import pi
from datetime import datetime
from .ellipsoid import Ellipsoid
@@ -29,10 +16,10 @@ from .utils import sanitize
try:
from .eci import eci2ecef, ecef2eci
except ImportError:
- eci2ecef = ecef2eci = None
+ eci2ecef = ecef2eci = None # type: ignore
-# py < 3.6 compatible
-tau = 2 * pi
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
__all__ = [
=====================================
src/pymap3d/eci.py
=====================================
@@ -2,8 +2,9 @@
from __future__ import annotations
+import typing
from datetime import datetime
-from numpy import ndarray, array, sin, cos, column_stack, empty, atleast_1d
+from numpy import array, sin, cos, column_stack, empty, atleast_1d
try:
from astropy.coordinates import GCRS, ITRS, EarthLocation, CartesianRepresentation
@@ -14,16 +15,14 @@ except ImportError:
from .sidereal import greenwichsrt, juliandate
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
__all__ = ["eci2ecef", "ecef2eci"]
def eci2ecef(
- x: float | ndarray,
- y: float | ndarray,
- z: float | ndarray,
- time: datetime,
- *,
- use_astropy: bool = True
+ x: ndarray, y: ndarray, z: ndarray, time: datetime, *, use_astropy: bool = True
) -> tuple[ndarray, ndarray, ndarray]:
"""
Observer => Point ECI => ECEF
@@ -81,12 +80,7 @@ def eci2ecef(
def ecef2eci(
- x: float | ndarray,
- y: float | ndarray,
- z: float | ndarray,
- time: datetime,
- *,
- use_astropy: bool = True
+ x: ndarray, y: ndarray, z: ndarray, time: datetime, *, use_astropy: bool = True
) -> tuple[ndarray, ndarray, ndarray]:
"""
Point => Point ECEF => ECI
=====================================
src/pymap3d/ellipsoid.py
=====================================
@@ -22,6 +22,7 @@ class Ellipsoid:
model : str
name of ellipsoid
"""
+
if model == "wgs84":
"""https://en.wikipedia.org/wiki/World_Geodetic_System#WGS84"""
self.semimajor_axis = 6378137.0
@@ -62,9 +63,7 @@ class Ellipsoid:
self.semiminor_axis = self.semimajor_axis
else:
raise NotImplementedError(
- "{} model not implemented, let us know and we will add it (or make a pull request)".format(
- model
- )
+ f"{model} model not implemented, let us know and we will add it (or make a pull request)"
)
self.flattening = (self.semimajor_axis - self.semiminor_axis) / self.semimajor_axis
=====================================
src/pymap3d/enu.py
=====================================
@@ -5,15 +5,16 @@ import typing
from math import tau
try:
- from numpy import asarray, radians, sin, cos, hypot, arctan2 as atan2, degrees, ndarray
+ from numpy import asarray, radians, sin, cos, hypot, arctan2 as atan2, degrees
except ImportError:
from math import radians, sin, cos, hypot, atan2, degrees # type: ignore
- ndarray = typing.Any # type: ignore
-
from .ecef import geodetic2ecef, ecef2geodetic, enu2ecef, uvw2enu
from .ellipsoid import Ellipsoid
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
__all__ = ["enu2aer", "aer2enu", "enu2geodetic", "geodetic2enu"]
@@ -46,7 +47,7 @@ def enu2aer(
slant range [meters]
"""
- # 1 millimeter precision for singularity
+ # 1 millimeter precision for singularity stability
try:
e[abs(e) < 1e-3] = 0.0
=====================================
src/pymap3d/latitude.py
=====================================
@@ -4,20 +4,23 @@ from __future__ import annotations
import typing
from .ellipsoid import Ellipsoid
-from .utils import sanitize
+from .utils import sanitize, sign
from . import rcurve
try:
- from numpy import radians, degrees, tan, sin, exp, pi, sqrt, inf, ndarray
+ from numpy import radians, degrees, tan, sin, cos, exp, pi, sqrt, inf
from numpy import arctan as atan, arcsinh as asinh, arctanh as atanh # noqa: A001
use_numpy = True
except ImportError:
- from math import atan, radians, degrees, tan, sin, asinh, atanh, exp, pi, sqrt, inf # type: ignore
+ from math import atan, radians, degrees, tan, sin, cos, asinh, atanh, exp, pi, sqrt, inf # type: ignore
use_numpy = False
- ndarray = typing.Any # type: ignore
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
+COS_EPS = 1e-9 # tolerance for angles near abs([90, 270])
__all__ = [
"geodetic2isometric",
@@ -85,8 +88,8 @@ def geoc2geod(
def geodetic2geocentric(
- geodetic_lat: float, alt_m: float, ell: Ellipsoid = None, deg: bool = True
-) -> float:
+ geodetic_lat: ndarray, alt_m: float, ell: Ellipsoid = None, deg: bool = True
+) -> ndarray:
"""
convert geodetic latitude to geocentric latitude on spheroid surface
@@ -126,8 +129,8 @@ geod2geoc = geodetic2geocentric
def geocentric2geodetic(
- geocentric_lat: float, alt_m: float, ell: Ellipsoid = None, deg: bool = True
-) -> float:
+ geocentric_lat: ndarray, alt_m: float, ell: Ellipsoid = None, deg: bool = True
+) -> ndarray:
"""
converts from geocentric latitude to geodetic latitude
@@ -163,7 +166,9 @@ def geocentric2geodetic(
return degrees(geodetic_lat) if deg else geodetic_lat
-def geodetic2isometric(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
+def geodetic2isometric(
+ geodetic_lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True
+) -> float:
"""
computes isometric latitude on an ellipsoid
@@ -205,19 +210,25 @@ def geodetic2isometric(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool =
# 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)))
+ coslat = cos(geodetic_lat)
+ i = abs(coslat) <= COS_EPS
+
try:
- isometric_lat[abs(geodetic_lat - pi / 2) <= 1e-9] = inf
- isometric_lat[abs(-geodetic_lat - pi / 2) <= 1e-9] = -inf
+ isometric_lat[i] = sign(geodetic_lat[i]) * inf # type: ignore
except TypeError:
- if abs(geodetic_lat - pi / 2) <= 1e-9:
- isometric_lat = inf
- elif abs(-geodetic_lat - pi / 2) <= 1e-9:
- isometric_lat = -inf
+ if i:
+ isometric_lat = sign(geodetic_lat) * inf
- return degrees(isometric_lat) if deg else isometric_lat
+ if deg:
+ isometric_lat = degrees(isometric_lat)
+
+ try:
+ return isometric_lat.squeeze()[()]
+ except AttributeError:
+ return isometric_lat
-def isometric2geodetic(isometric_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def isometric2geodetic(isometric_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""
converts from isometric latitude to geodetic latitude
@@ -253,7 +264,7 @@ def isometric2geodetic(isometric_lat: float, ell: Ellipsoid = None, deg: bool =
return degrees(geodetic_lat) if deg else geodetic_lat
-def conformal2geodetic(conformal_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def conformal2geodetic(conformal_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""
converts from conformal latitude to geodetic latitude
@@ -298,7 +309,7 @@ def conformal2geodetic(conformal_lat: float, ell: Ellipsoid = None, deg: bool =
return degrees(geodetic_lat) if deg else geodetic_lat
-def geodetic2conformal(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def geodetic2conformal(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""
converts from geodetic latitude to conformal latitude
@@ -344,7 +355,9 @@ def geodetic2conformal(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = T
# %% rectifying
-def geodetic2rectifying(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
+def geodetic2rectifying(
+ geodetic_lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True
+) -> float:
"""
converts from geodetic latitude to rectifying latitude
@@ -438,9 +451,7 @@ def rectifying2geodetic(
# %% authalic
-def geodetic2authalic(
- geodetic_lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True
-) -> float:
+def geodetic2authalic(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""
converts from geodetic latitude to authalic latitude
@@ -484,9 +495,7 @@ def geodetic2authalic(
return degrees(authalic_lat) if deg else authalic_lat
-def authalic2geodetic(
- authalic_lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True
-) -> float:
+def authalic2geodetic(authalic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""
converts from authalic latitude to geodetic latitude
@@ -529,7 +538,7 @@ def authalic2geodetic(
# %% parametric
-def geodetic2parametric(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def geodetic2parametric(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""
converts from geodetic latitude to parametric latitude
@@ -563,7 +572,9 @@ def geodetic2parametric(geodetic_lat: float, ell: Ellipsoid = None, deg: bool =
return degrees(parametric_lat) if deg else parametric_lat
-def parametric2geodetic(parametric_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def parametric2geodetic(
+ parametric_lat: ndarray, ell: Ellipsoid = None, deg: bool = True
+) -> ndarray:
"""
converts from parametric latitude to geodetic latitude
=====================================
src/pymap3d/los.py
=====================================
@@ -4,16 +4,17 @@ from __future__ import annotations
import typing
try:
- from numpy import pi, nan, sqrt, ndarray, atleast_1d
+ from numpy import pi, nan, sqrt, atleast_1d
except ImportError:
from math import pi, nan, sqrt # type: ignore
- ndarray = typing.Any # type: ignore
-
from .aer import aer2enu
from .ecef import enu2uvw, geodetic2ecef, ecef2geodetic
from .ellipsoid import Ellipsoid
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
__all__ = ["lookAtSpheroid"]
=====================================
src/pymap3d/lox.py
=====================================
@@ -4,14 +4,14 @@ from __future__ import annotations
import typing
try:
- from numpy import radians, degrees, cos, arctan2 as atan2, tan, pi, ndarray, atleast_1d
+ from numpy import radians, degrees, cos, arctan2 as atan2, tan, array, broadcast_arrays
except ImportError:
- from math import radians, degrees, cos, atan2, tan, pi # type: ignore
+ from math import radians, degrees, cos, atan2, tan # type: ignore
- ndarray = typing.Any # type: ignore
+from math import pi, tau
-import typing
from .ellipsoid import Ellipsoid
+from .utils import sign
from . import rcurve
from . import rsphere
from .latitude import (
@@ -23,6 +23,9 @@ from .latitude import (
)
from .utils import sph2cart, cart2sph
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
__all__ = [
"loxodrome_inverse",
"loxodrome_direct",
@@ -33,6 +36,9 @@ __all__ = [
]
+COS_EPS = 1e-9
+
+
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.
@@ -138,6 +144,12 @@ def loxodrome_inverse(
if deg:
lat1, lon1, lat2, lon2 = radians(lat1), radians(lon1), radians(lat2), radians(lon2)
+ try:
+ lat1, lon1, lat2, lon2 = broadcast_arrays(lat1, lon1, lat2, lon2)
+
+ except NameError:
+ pass
+
# compute changes in isometric latitude and longitude between points
disolat = geodetic2isometric(lat2, deg=False, ell=ell) - geodetic2isometric(
lat1, deg=False, ell=ell
@@ -146,35 +158,41 @@ def loxodrome_inverse(
# compute azimuth
az12 = atan2(dlon, disolat)
- cosaz12 = cos(az12)
+ aux = abs(cos(az12))
# compute distance along loxodromic curve
- dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / abs(cos(az12))
+ dist = meridian_arc(lat2, lat1, deg=False, ell=ell) / aux
+
+ # straight east or west
+ i = aux < COS_EPS
try:
- if (abs(cosaz12) < 1e-9).any():
- dist[abs(cosaz12) < 1e-9] = departure(lon2, lon1, lat1, ell, deg=False)
+ dist[i] = departure(lon2[i], lon1[i], lat1[i], ell, deg=False)
except (AttributeError, TypeError):
- if abs(cosaz12) < 1e-9: # straight east or west
+ if i:
dist = departure(lon2, lon1, lat1, ell, deg=False)
if deg:
az12 = degrees(az12) % 360.0
- return dist, az12
+ try:
+ return dist.squeeze()[()], az12.squeeze()[()]
+ except AttributeError:
+ return dist, az12
def loxodrome_direct(
- lat1: ndarray,
- lon1: ndarray,
- rng: ndarray,
+ lat1: float | ndarray,
+ lon1: float | ndarray,
+ rng: float | ndarray,
a12: float,
ell: Ellipsoid = None,
deg: bool = True,
-) -> tuple[ndarray, ndarray]:
+) -> tuple[float | ndarray, float | ndarray]:
"""
Given starting lat, lon with arclength and azimuth, compute final lat, lon
like Matlab reckon('rh', ...)
+ except that "rng" in meters instead of "arclen" degrees of arc
Parameters
----------
@@ -202,15 +220,16 @@ def loxodrome_direct(
if deg:
lat1, lon1, a12 = radians(lat1), radians(lon1), radians(a12)
+ a12 = a12 % tau
+
try:
- lat1 = atleast_1d(lat1)
- rng = atleast_1d(rng)
- if (abs(lat1) > pi / 2).any():
+ lat1, rng, a12 = broadcast_arrays(lat1, rng, a12)
+ if (abs(lat1) > pi / 2).any(): # type: ignore
raise ValueError("-90 <= latitude <= 90")
- if (rng < 0).any():
+ if (rng < 0).any(): # type: ignore
raise ValueError("ground distance must be >= 0")
except NameError:
- if abs(lat1) > pi / 2:
+ if abs(lat1) > pi / 2: # type: ignore
raise ValueError("-90 <= latitude <= 90")
if rng < 0:
raise ValueError("ground distance must be >= 0")
@@ -226,14 +245,23 @@ def loxodrome_direct(
newiso = geodetic2isometric(lat2, ell, deg=False)
iso = geodetic2isometric(lat1, ell, deg=False)
+ # stability near singularities
+ i = abs(cos(a12)) < COS_EPS
dlon = tan(a12) * (newiso - iso)
+
+ try:
+ dlon[i] = sign(pi - a12[i]) * rng[i] / rcurve.parallel(lat1[i], ell=ell, deg=False) # type: ignore
+ except (AttributeError, TypeError):
+ if i: # straight east or west
+ dlon = sign(pi - a12) * rng / rcurve.parallel(lat1, ell=ell, deg=False)
+
lon2 = lon1 + dlon
if deg:
lat2, lon2 = degrees(lat2), degrees(lon2)
try:
- return lat2.squeeze()[()], lon2.squeeze()[()]
+ return lat2.squeeze()[()], lon2.squeeze()[()] # type: ignore
except AttributeError:
return lat2, lon2
@@ -265,12 +293,12 @@ def departure(
if deg:
lon1, lon2, lat = radians(lon1), radians(lon2), radians(lat)
- return rcurve.parallel(lat, ell, deg=False) * ((lon2 - lon1) % pi)
+ return rcurve.parallel(lat, ell=ell, deg=False) * ((lon2 - lon1) % pi)
def meanm(
lat: ndarray, lon: ndarray, ell: Ellipsoid = None, deg: bool = True
-) -> tuple[float | ndarray, float | ndarray]:
+) -> tuple[ndarray, ndarray]:
"""
Computes geographic mean for geographic points on an ellipsoid
@@ -297,8 +325,8 @@ def meanm(
lat, lon = radians(lat), radians(lon)
lat = geodetic2authalic(lat, ell, deg=False)
- assert isinstance(lat, ndarray)
- x, y, z = sph2cart(lon, lat, 1)
+
+ x, y, z = sph2cart(lon, lat, array(1.0))
lonbar, latbar, _ = cart2sph(x.sum(), y.sum(), z.sum())
latbar = authalic2geodetic(latbar, ell, deg=False)
=====================================
src/pymap3d/ned.py
=====================================
@@ -3,15 +3,13 @@
from __future__ import annotations
import typing
-try:
- from numpy import ndarray
-except ImportError:
- ndarray = typing.Any # type: ignore
-
from .enu import geodetic2enu, aer2enu, enu2aer
from .ecef import ecef2geodetic, ecef2enuv, ecef2enu, enu2ecef
from .ellipsoid import Ellipsoid
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
def aer2ned(
az: ndarray, elev: ndarray, slantRange: ndarray, deg: bool = True
=====================================
src/pymap3d/rcurve.py
=====================================
@@ -4,19 +4,20 @@ from __future__ import annotations
import typing
try:
- from numpy import radians, sin, cos, sqrt, ndarray
+ from numpy import sin, cos, sqrt
except ImportError:
- from math import radians, sin, cos, sqrt # type: ignore
-
- ndarray = typing.Any # type: ignore
+ from math import sin, cos, sqrt # type: ignore
from .ellipsoid import Ellipsoid
from .utils import sanitize
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
__all__ = ["parallel", "meridian", "transverse", "geocentric_radius"]
-def geocentric_radius(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def geocentric_radius(geodetic_lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""
compute geocentric radius at geodetic latitude
@@ -36,10 +37,12 @@ def geocentric_radius(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = Tr
)
-def parallel(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
+def parallel(lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
"""
computes the radius of the small circle encompassing the globe at the specified latitude
+ like Matlab rcurve('parallel', ...)
+
Parameters
----------
lat : float
@@ -52,16 +55,15 @@ def parallel(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
Returns
-------
radius: float
- radius of ellipsoid
+ radius of ellipsoid (meters)
"""
- if deg:
- lat = radians(lat)
+ lat, ell = sanitize(lat, ell, deg)
return cos(lat) * transverse(lat, ell, deg=False)
-def meridian(lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+def meridian(lat: ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""computes the meridional radius of curvature for the ellipsoid
like Matlab rcurve('meridian', ...)
@@ -81,17 +83,14 @@ def meridian(lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
radius of ellipsoid
"""
- if ell is None:
- ell = Ellipsoid()
- if deg:
- lat = radians(lat)
+ lat, ell = sanitize(lat, ell, deg)
f1 = ell.semimajor_axis * (1 - ell.eccentricity ** 2)
f2 = 1 - (ell.eccentricity * sin(lat)) ** 2
return f1 / sqrt(f2 ** 3)
-def transverse(lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True) -> float:
+def transverse(lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True) -> ndarray:
"""computes the radius of the curve formed by a plane
intersecting the ellipsoid at the latitude which is
normal to the surface of the ellipsoid
@@ -110,7 +109,7 @@ def transverse(lat: float | ndarray, ell: Ellipsoid = None, deg: bool = True) ->
Returns
-------
radius: float
- radius of ellipsoid
+ radius of ellipsoid (meters)
"""
lat, ell = sanitize(lat, ell, deg)
=====================================
src/pymap3d/rsphere.py
=====================================
@@ -4,17 +4,18 @@ from __future__ import annotations
import typing
try:
- from numpy import radians, sin, cos, log, sqrt, degrees, asarray, ndarray
+ from numpy import radians, sin, cos, log, sqrt, degrees, asarray
except ImportError:
from math import radians, sin, cos, log, sqrt, degrees # type: ignore
- asarray = None
- ndarray = typing.Any # type: ignore
+ asarray = None # type: ignore
from .ellipsoid import Ellipsoid
from . import rcurve
from .vincenty import vdist
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
__all__ = [
"eqavol",
@@ -141,7 +142,7 @@ def euler(
return rho * nu / den
-def curve(lat: float, ell: Ellipsoid = None, deg: bool = True, method: str = "mean") -> float:
+def curve(lat: ndarray, ell: Ellipsoid = None, deg: bool = True, method: str = "mean") -> ndarray:
"""computes the arithmetic average of the transverse and meridional
radii of curvature at a specified latitude point
=====================================
src/pymap3d/sidereal.py
=====================================
@@ -1,6 +1,6 @@
# Copyright (c) 2014-2018 Michael Hirsch, Ph.D.
""" manipulations of sidereal time """
-from math import pi
+from math import tau
from datetime import datetime
from .timeconv import str2dt
@@ -128,4 +128,4 @@ def greenwichsrt(Jdate: float) -> float:
)
# 1/86400 and %(2*pi) implied by units of radians
- return gmst_sec * (2 * pi) / 86400.0 % (2 * pi)
+ return gmst_sec * tau / 86400.0 % tau
=====================================
src/pymap3d/tests/test_latitude.py
=====================================
@@ -49,7 +49,14 @@ def test_numpy_geodetic_geocentric():
@pytest.mark.parametrize(
"geodetic_lat, isometric_lat",
- [(0, 0), (90, inf), (-90, -inf), (45, 50.227466), (-45, -50.227466), (89, 271.275)],
+ [
+ (0, 0),
+ (90, inf),
+ (-90, -inf),
+ (45, 50.227466),
+ (-45, -50.227466),
+ (89, 271.275),
+ ],
)
def test_geodetic_isometric(geodetic_lat, isometric_lat):
isolat = latitude.geodetic2isometric(geodetic_lat)
=====================================
src/pymap3d/tests/test_rcurve.py
=====================================
@@ -41,4 +41,4 @@ def test_numpy_meridian():
def test_numpy_transverse():
pytest.importorskip("numpy")
- assert rcurve.transverse([0, 90]) == approx([A, 6399593.626])
+ assert rcurve.transverse([-90, 0, 90]) == approx([6399593.6258, A, 6399593.6258])
=====================================
src/pymap3d/tests/test_rhumb.py
=====================================
@@ -69,15 +69,36 @@ def test_loxodrome_inverse(lat1, lon1, lat2, lon2, arclen, az):
def test_numpy_loxodrome_inverse():
pytest.importorskip("numpy")
d, a = lox.loxodrome_inverse([40, 40], [-80, -80], 65, -148)
+ assert d == approx(5248666.209)
+ assert a == approx(302.00567)
+
+ d, a = lox.loxodrome_inverse([40, 40], [-80, -80], [65, 65], -148)
+ d, a = lox.loxodrome_inverse([40, 40], [-80, -80], 65, [-148, -148])
+
+def test_numpy_2d_loxodrome_inverse():
+ pytest.importorskip("numpy")
+ d, a = lox.loxodrome_inverse([[40, 40], [40, 40]], [[-80, -80], [-80, -80]], 65, -148)
assert d == approx(5248666.209)
assert a == approx(302.00567)
+ d, a = lox.loxodrome_inverse(
+ [[40, 40], [40, 40]], [[-80, -80], [-80, -80]], [[65, 65], [65, 65]], -148
+ )
+ d, a = lox.loxodrome_inverse(
+ [[40, 40], [40, 40]], [[-80, -80], [-80, -80]], 65, [[-148, -148], [-148, -148]]
+ )
+ d, a = lox.loxodrome_inverse(40, -80, [[65, 65], [65, 65]], [[-148, -148], [-148, -148]])
+
@pytest.mark.parametrize(
"lat0,lon0,rng,az,lat1,lon1",
[
- (40, -80, 10000, 30, 40.077995, -79.9414144),
+ (40, -80, 10000, 30, 40.0000779959676, -79.9999414477481),
+ (35, 140, 50000, 90, 35, 140.548934481815),
+ (35, 140, 50000, -270, 35, 140.548934481815),
+ (35, 140, 50000, 270, 35, 139.451065518185),
+ (35, 140, 50000, -90, 35, 139.451065518185),
(0, 0, 0, 0, 0, 0),
(0, 0, 10018754.17, 90, 0, 90),
(0, 0, 10018754.17, -90, 0, -90),
@@ -86,18 +107,25 @@ def test_numpy_loxodrome_inverse():
],
)
def test_loxodrome_direct(lat0, lon0, rng, az, lat1, lon1):
+ """
+ reckon('rh', 40, -80, 10, 30, wgs84Ellipsoid)
+ """
lat2, lon2 = lox.loxodrome_direct(lat0, lon0, rng, az)
- assert lat2 == approx(lat1, abs=1e-6)
- assert lon2 == approx(lon1)
+ assert lat2 == approx(lat1, rel=0.005, abs=1e-6)
+ assert lon2 == approx(lon1, rel=0.001)
assert isinstance(lat2, float)
assert isinstance(lon2, float)
def test_numpy_loxodrome_direct():
pytest.importorskip("numpy")
- lat, lon = lox.loxodrome_direct([40, 40], [-80, -80], [10000, 10000], [30, 30])
- assert lat == approx(40.077995)
- assert lon == approx(-79.941414)
+ lat, lon = lox.loxodrome_direct([40, 40], [-80, -80], [10, 10], [30, 30])
+ assert lat == approx(40.000078)
+ assert lon == approx(-79.99994145)
+
+ lat, lon = lox.loxodrome_direct([40, 40], [-80, -80], 10, 30)
+ lat, lon = lox.loxodrome_direct([40, 40], [-80, -80], [10, 10], 30)
+ lat, lon = lox.loxodrome_direct([40, 40], [-80, -80], 10, [30, 30])
@pytest.mark.parametrize("lat,lon", [([0, 45, 90], [0, 45, 90])])
=====================================
src/pymap3d/timeconv.py
=====================================
@@ -5,11 +5,11 @@ from datetime import datetime
try:
from dateutil.parser import parse
except ImportError:
- parse = None
+ parse = None # type: ignore
try:
import numpy as np
except ImportError:
- np = None
+ np = None # type: ignore
__all__ = ["str2dt"]
=====================================
src/pymap3d/utils.py
=====================================
@@ -8,12 +8,11 @@ import typing
from .ellipsoid import Ellipsoid
try:
- from numpy import hypot, cos, sin, arctan2 as atan2, radians, pi, asarray, ndarray, sign
+ from numpy import hypot, cos, sin, arctan2 as atan2, radians, pi, asarray, sign
except ImportError:
from math import atan2, hypot, cos, sin, radians, pi # type: ignore
- asarray = None
- ndarray = typing.Any # type: ignore
+ asarray = None # type: ignore
def sign(x: float) -> float: # type: ignore
"""signum function"""
@@ -27,6 +26,9 @@ except ImportError:
return y
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
__all__ = ["cart2pol", "pol2cart", "cart2sph", "sph2cart", "sign"]
@@ -40,9 +42,7 @@ def pol2cart(theta: float, rho: float) -> tuple[float, float]:
return rho * cos(theta), rho * sin(theta)
-def cart2sph(
- x: float | ndarray, y: float | ndarray, z: float | ndarray
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+def cart2sph(x: ndarray, y: ndarray, z: ndarray) -> tuple[ndarray, ndarray, ndarray]:
"""Transform Cartesian to spherical coordinates"""
hxy = hypot(x, y)
r = hypot(hxy, z)
@@ -51,9 +51,7 @@ def cart2sph(
return az, el, r
-def sph2cart(
- az: float | ndarray, el: float | ndarray, r: float | ndarray
-) -> tuple[float | ndarray, float | ndarray, float | ndarray]:
+def sph2cart(az: ndarray, el: ndarray, r: ndarray) -> tuple[ndarray, ndarray, ndarray]:
"""Transform spherical to Cartesian coordinates"""
rcos_theta = r * cos(el)
x = rcos_theta * cos(az)
@@ -62,7 +60,9 @@ def sph2cart(
return x, y, z
-def sanitize(lat: float | ndarray, ell: Ellipsoid, deg: bool) -> tuple[float | ndarray, Ellipsoid]:
+def sanitize(
+ lat: float | ndarray, ell: typing.Optional[Ellipsoid], deg: bool
+) -> tuple[float | ndarray, Ellipsoid]:
if ell is None:
ell = Ellipsoid()
if asarray is not None:
@@ -71,10 +71,10 @@ def sanitize(lat: float | ndarray, ell: Ellipsoid, deg: bool) -> tuple[float | n
lat = radians(lat)
if asarray is not None:
- if (abs(lat) > pi / 2).any():
+ if (abs(lat) > pi / 2).any(): # type: ignore
raise ValueError("-pi/2 <= latitude <= pi/2")
else:
- if abs(lat) > pi / 2:
+ if abs(lat) > pi / 2: # type: ignore
raise ValueError("-pi/2 <= latitude <= pi/2")
return lat, ell
=====================================
src/pymap3d/vincenty.py
=====================================
@@ -10,7 +10,6 @@ from copy import copy
try:
from numpy import (
- ndarray,
atleast_1d,
sqrt,
tan,
@@ -26,11 +25,12 @@ try:
except ImportError:
from math import sqrt, tan, sin, cos, isnan, atan, atan2, asin, radians, degrees # type: ignore
- ndarray = typing.Any # type: ignore
-
from .ellipsoid import Ellipsoid
from .utils import sign
+if typing.TYPE_CHECKING:
+ from numpy import ndarray
+
__all__ = ["vdist", "vreckon", "track2"]
@@ -123,7 +123,7 @@ def vdist(
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):
+ if (abs(Lat1) > 90) | (abs(Lat2) > 90): # type: ignore
raise ValueError("Input latitudes must be in [-90, 90] degrees.")
# %% Supply WGS84 earth ellipsoid axis lengths in meters:
a = ell.semimajor_axis
@@ -158,7 +158,7 @@ def vdist(
L[L > pi] = 2 * pi - L[L > pi]
except TypeError:
if L > pi:
- L = 2 * pi - L
+ L = 2 * pi - L # type: ignore
lamb = copy(L) # NOTE: program will fail without copy!
itercount = 0
@@ -170,7 +170,7 @@ def vdist(
if not warninggiven:
logging.warning("Essentially antipodal points--precision may be reduced slightly.")
- lamb = pi
+ lamb = pi # type: ignore
break
lambdaold = copy(lamb)
@@ -216,7 +216,7 @@ def vdist(
# print(f'then, lambda(21752) = {lamb[21752],20})
# correct for convergence failure for essentially antipodal points
try:
- i = (lamb > pi).any()
+ i = (lamb > pi).any() # type: ignore
except AttributeError:
i = lamb > pi
@@ -225,11 +225,11 @@ def vdist(
"Essentially antipodal points encountered. Precision may be reduced slightly."
)
warninggiven = True
- lambdaold = pi
- lamb = pi
+ lambdaold = pi # type: ignore
+ lamb = pi # type: ignore
try:
- notdone = (abs(lamb - lambdaold) > 1e-12).any()
+ notdone = (abs(lamb - lambdaold) > 1e-12).any() # type: ignore
except AttributeError:
notdone = abs(lamb - lambdaold) > 1e-12
@@ -341,14 +341,14 @@ def vreckon(
Lon1 = atleast_1d(Lon1)
Rng = atleast_1d(Rng)
Azim = atleast_1d(Azim)
- if (abs(Lat1) > 90).any():
+ if (abs(Lat1) > 90.0).any():
raise ValueError("Input lat. must be between -90 and 90 deg., inclusive.")
- if (Rng < 0).any():
+ if (Rng < 0.0).any():
raise ValueError("Ground distance must be positive")
except NameError:
- if abs(Lat1) > 90:
+ if abs(Lat1) > 90.0: # type: ignore
raise ValueError("Input lat. must be between -90 and 90 deg., inclusive.")
- if Rng < 0:
+ if Rng < 0.0:
raise ValueError("Ground distance must be positive")
if ell is not None:
=====================================
src/pymap3d/vreckon/__main__.py
=====================================
@@ -13,4 +13,4 @@ P = p.parse_args()
lat2, lon2 = vincenty.vreckon(P.lat, P.lon, P.range, P.azimuth)
-print("lat, lon = ({:.4f}, {:.4f})".format(lat2.item(), lon2.item()))
+print(f"lat, lon = ({lat2.item():.4f}, {lon2.item():.4f})")
View it on GitLab: https://salsa.debian.org/debian-gis-team/pymap3d/-/commit/6e091f6a78b215283bec1d611c6593006b5ff35f
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/pymap3d/-/commit/6e091f6a78b215283bec1d611c6593006b5ff35f
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/20211020/c32d903d/attachment-0001.htm>
More information about the Pkg-grass-devel
mailing list