diff --git a/Examples/plot_geodetic2ecef.py b/Examples/plot_geodetic2ecef.py index 16c9487..a95876d 100644 --- a/Examples/plot_geodetic2ecef.py +++ b/Examples/plot_geodetic2ecef.py @@ -1,5 +1,7 @@ #!/usr/bin/env python3 +from __future__ import annotations +import typing import argparse import matplotlib.pyplot as mpl @@ -15,7 +17,7 @@ x, y, z = pm.geodetic2ecef(lat, lon, args.alt_m) -def panel(ax, val, name: str, cmap: str = None): +def panel(ax, val, name: str, cmap: str | None = None): hi = ax.pcolormesh(lon, lat, val, cmap=cmap) ax.set_title(name) fg.colorbar(hi, ax=ax).set_label(name + " [m]") @@ -23,7 +25,7 @@ def panel(ax, val, name: str, cmap: str = None): fg = mpl.figure(figsize=(16, 5)) -axs = fg.subplots(1, 3, sharey=True) +axs: typing.Any = fg.subplots(1, 3, sharey=True) fg.suptitle("geodetic2ecef") panel(axs[0], x, "x", "bwr") diff --git a/pyproject.toml b/pyproject.toml index eaffd35..f6e1012 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -42,8 +42,6 @@ known_third_party = ["pymap3d"] files = ["src", "Examples", "scripts"] ignore_missing_imports = true -strict_optional = false -show_column_numbers = true [tool.coverage.run] branch = true diff --git a/src/pymap3d/__init__.py b/src/pymap3d/__init__.py index d49250a..6f4575b 100644 --- a/src/pymap3d/__init__.py +++ b/src/pymap3d/__init__.py @@ -29,7 +29,7 @@ * Fortran: [Maptran3D](https://github.com/geospace-code/maptran3d) """ -__version__ = "3.0.1" +__version__ = "3.1.0" from .aer import aer2ecef, aer2geodetic, ecef2aer, geodetic2aer from .ecef import ( diff --git a/src/pymap3d/aer.py b/src/pymap3d/aer.py index 1c8128d..631345b 100644 --- a/src/pymap3d/aer.py +++ b/src/pymap3d/aer.py @@ -15,6 +15,8 @@ __all__ = ["aer2ecef", "ecef2aer", "geodetic2aer", "aer2geodetic", "eci2aer", "aer2eci"] +ELL = Ellipsoid.from_name("wgs84") + def ecef2aer( x, @@ -23,7 +25,7 @@ def ecef2aer( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -72,7 +74,7 @@ def geodetic2aer( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -120,7 +122,7 @@ def aer2geodetic( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -213,7 +215,7 @@ def aer2eci( lon0, h0, t: datetime, - ell=None, + ell: Ellipsoid = ELL, *, deg: bool = True, ) -> tuple: @@ -269,7 +271,7 @@ def aer2ecef( lat0, lon0, alt0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ diff --git a/src/pymap3d/ecef.py b/src/pymap3d/ecef.py index 4c7d2c4..88423a0 100644 --- a/src/pymap3d/ecef.py +++ b/src/pymap3d/ecef.py @@ -16,7 +16,6 @@ from .ellipsoid import Ellipsoid from .mathfun import atan, atan2, cos, degrees, hypot, isclose, radians, sin, sqrt, tan -from .utils import sanitize __all__ = [ "geodetic2ecef", @@ -30,12 +29,14 @@ "enu2ecef", ] +ELL = Ellipsoid.from_name("wgs84") + def geodetic2ecef( lat, lon, alt, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -68,8 +69,9 @@ def geodetic2ecef( z target z ECEF coordinate (meters) """ - lat, ell = sanitize(lat, ell, deg) + if deg: + lat = radians(lat) lon = radians(lon) # radius of curvature of the prime vertical section @@ -88,7 +90,7 @@ def ecef2geodetic( x, y, z, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -121,9 +123,6 @@ def ecef2geodetic( Journal of Surveying Engineering. doi: 10.1061/(ASCE)0733-9453 """ - if ell is None: - ell = Ellipsoid.from_name("wgs84") - try: x = asarray(x) y = asarray(y) @@ -167,10 +166,7 @@ def ecef2geodetic( # eqn. 13 Beta += ( (ell.semiminor_axis * u - ell.semimajor_axis * huE + E**2) * sin(Beta) - ) / ( - ell.semimajor_axis * huE * 1 / cos(Beta) - - E**2 * cos(Beta) - ) + ) / (ell.semimajor_axis * huE * 1 / cos(Beta) - E**2 * cos(Beta)) except (ArithmeticError, RuntimeWarning): if isclose(z, 0): Beta = 0 @@ -281,7 +277,7 @@ def ecef2enu( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -407,7 +403,7 @@ def uvw2enu(u, v, w, lat0, lon0, deg: bool = True) -> tuple: return East, North, Up -def eci2geodetic(x, y, z, t: datetime, ell: Ellipsoid = None, *, deg: bool = True) -> tuple: +def eci2geodetic(x, y, z, t: datetime, ell: Ellipsoid = ELL, *, deg: bool = True) -> tuple: """ convert Earth Centered Internal ECI to geodetic coordinates @@ -448,7 +444,7 @@ def eci2geodetic(x, y, z, t: datetime, ell: Ellipsoid = None, *, deg: bool = Tru return ecef2geodetic(xecef, yecef, zecef, ell, deg) -def geodetic2eci(lat, lon, alt, t: datetime, ell: Ellipsoid = None, *, deg: bool = True) -> tuple: +def geodetic2eci(lat, lon, alt, t: datetime, ell: Ellipsoid = ELL, *, deg: bool = True) -> tuple: """ convert geodetic coordinates to Earth Centered Internal ECI @@ -496,7 +492,7 @@ def enu2ecef( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -532,6 +528,7 @@ def enu2ecef( z target z ECEF coordinate (meters) """ + x0, y0, z0 = geodetic2ecef(lat0, lon0, h0, ell, deg=deg) dx, dy, dz = enu2uvw(e1, n1, u1, lat0, lon0, deg=deg) diff --git a/src/pymap3d/ellipsoid.py b/src/pymap3d/ellipsoid.py index fa0b9f6..72ea9a6 100644 --- a/src/pymap3d/ellipsoid.py +++ b/src/pymap3d/ellipsoid.py @@ -151,7 +151,7 @@ def __init__( } @classmethod - def from_name(cls, name: str) -> Ellipsoid | None: + def from_name(cls, name: str) -> Ellipsoid: """Create an Ellipsoid from a name.""" return cls( diff --git a/src/pymap3d/enu.py b/src/pymap3d/enu.py index c445ad2..1e90086 100644 --- a/src/pymap3d/enu.py +++ b/src/pymap3d/enu.py @@ -14,6 +14,8 @@ __all__ = ["enu2aer", "aer2enu", "enu2geodetic", "geodetic2enu"] +ELL = Ellipsoid.from_name("wgs84") + def enu2aer(e, n, u, deg: bool = True) -> tuple: """ @@ -115,7 +117,7 @@ def enu2geodetic( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -163,7 +165,7 @@ def geodetic2enu( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ diff --git a/src/pymap3d/latitude.py b/src/pymap3d/latitude.py index adb11cc..e58658a 100644 --- a/src/pymap3d/latitude.py +++ b/src/pymap3d/latitude.py @@ -22,8 +22,8 @@ sqrt, tan, ) -from .utils import sanitize +ELL = Ellipsoid.from_name("wgs84") COS_EPS = 1e-9 # tolerance for angles near abs([90, 270]) __all__ = [ @@ -47,7 +47,7 @@ def geoc2geod( geocentric_lat, geocentric_distance, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ): """ @@ -78,7 +78,9 @@ def geoc2geod( and Geodetic Coordinates. Celestial Mechanics (12), 2, p. 225-230 (1975) doi: 10.1007/BF01230214" """ - geocentric_lat, ell = sanitize(geocentric_lat, ell, deg) + + if deg: + geocentric_lat = radians(geocentric_lat) r = geocentric_distance / ell.semimajor_axis @@ -91,7 +93,7 @@ def geoc2geod( return degrees(geodetic_lat) if deg else geodetic_lat -def geodetic2geocentric(geodetic_lat, alt_m, ell: Ellipsoid = None, deg: bool = True): +def geodetic2geocentric(geodetic_lat, alt_m, ell: Ellipsoid = ELL, deg: bool = True): """ convert geodetic latitude to geocentric latitude on spheroid surface @@ -120,7 +122,10 @@ def geodetic2geocentric(geodetic_lat, alt_m, ell: Ellipsoid = None, deg: bool = US Geological Survey Professional Paper 1395, US Government Printing Office, Washington, DC, 1987, pp. 13-18. """ - geodetic_lat, ell = sanitize(geodetic_lat, ell, deg) + + if deg: + geodetic_lat = radians(geodetic_lat) + r = rcurve.transverse(geodetic_lat, ell, deg=False) geocentric_lat = atan((1 - ell.eccentricity**2 * (r / (r + alt_m))) * tan(geodetic_lat)) @@ -130,7 +135,7 @@ def geodetic2geocentric(geodetic_lat, alt_m, ell: Ellipsoid = None, deg: bool = geod2geoc = geodetic2geocentric -def geocentric2geodetic(geocentric_lat, alt_m, ell: Ellipsoid = None, deg: bool = True): +def geocentric2geodetic(geocentric_lat, alt_m, ell: Ellipsoid = ELL, deg: bool = True): """ converts from geocentric latitude to geodetic latitude @@ -159,14 +164,17 @@ def geocentric2geodetic(geocentric_lat, alt_m, ell: Ellipsoid = None, deg: bool US Geological Survey Professional Paper 1395, US Government Printing Office, Washington, DC, 1987, pp. 13-18. """ - geocentric_lat, ell = sanitize(geocentric_lat, ell, deg) + + if deg: + geocentric_lat = radians(geocentric_lat) + r = rcurve.transverse(geocentric_lat, ell, deg=False) geodetic_lat = atan(tan(geocentric_lat) / (1 - ell.eccentricity**2 * (r / (r + alt_m)))) return degrees(geodetic_lat) if deg else geodetic_lat -def geodetic2isometric(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): +def geodetic2isometric(geodetic_lat, ell: Ellipsoid = ELL, deg: bool = True): """ computes isometric latitude on an ellipsoid @@ -196,7 +204,8 @@ def geodetic2isometric(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): January 2010 """ - geodetic_lat, ell = sanitize(geodetic_lat, ell, deg) + if deg: + geodetic_lat = radians(geodetic_lat) e = ell.eccentricity @@ -226,7 +235,7 @@ def geodetic2isometric(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): return isometric_lat -def isometric2geodetic(isometric_lat, ell: Ellipsoid = None, deg: bool = True): +def isometric2geodetic(isometric_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from isometric latitude to geodetic latitude @@ -252,7 +261,7 @@ def isometric2geodetic(isometric_lat, ell: Ellipsoid = None, deg: bool = True): US Geological Survey Professional Paper 1395, US Government Printing Office, Washington, DC, 1987, pp. 13-18. """ - # NOT sanitize for isometric2geo + if deg: isometric_lat = radians(isometric_lat) @@ -262,7 +271,7 @@ def isometric2geodetic(isometric_lat, ell: Ellipsoid = None, deg: bool = True): return degrees(geodetic_lat) if deg else geodetic_lat -def conformal2geodetic(conformal_lat, ell: Ellipsoid = None, deg: bool = True): +def conformal2geodetic(conformal_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from conformal latitude to geodetic latitude @@ -288,7 +297,9 @@ def conformal2geodetic(conformal_lat, ell: Ellipsoid = None, deg: bool = True): US Geological Survey Professional Paper 1395, US Government Printing Office, Washington, DC, 1987, pp. 13-18. """ - conformal_lat, ell = sanitize(conformal_lat, ell, deg) + + if deg: + conformal_lat = radians(conformal_lat) e = ell.eccentricity f1 = e**2 / 2 + 5 * e**4 / 24 + e**6 / 12 + 13 * e**8 / 360 @@ -307,7 +318,7 @@ def conformal2geodetic(conformal_lat, ell: Ellipsoid = None, deg: bool = True): return degrees(geodetic_lat) if deg else geodetic_lat -def geodetic2conformal(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): +def geodetic2conformal(geodetic_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from geodetic latitude to conformal latitude @@ -335,7 +346,8 @@ def geodetic2conformal(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): """ - geodetic_lat, ell = sanitize(geodetic_lat, ell, deg) + if deg: + geodetic_lat = radians(geodetic_lat) e = ell.eccentricity f1 = 1 - e * sin(geodetic_lat) @@ -355,7 +367,7 @@ def geodetic2conformal(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): # %% rectifying -def geodetic2rectifying(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): +def geodetic2rectifying(geodetic_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from geodetic latitude to rectifying latitude @@ -382,7 +394,9 @@ def geodetic2rectifying(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): Office, Washington, DC, 1987, pp. 13-18. """ - geodetic_lat, ell = sanitize(geodetic_lat, ell, deg) + + if deg: + geodetic_lat = radians(geodetic_lat) n = ell.thirdflattening f1 = 3 * n / 2 - 9 * n**3 / 16 @@ -401,7 +415,7 @@ def geodetic2rectifying(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): return degrees(rectifying_lat) if deg else rectifying_lat -def rectifying2geodetic(rectifying_lat, ell: Ellipsoid = None, deg: bool = True): +def rectifying2geodetic(rectifying_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from rectifying latitude to geodetic latitude @@ -427,7 +441,9 @@ def rectifying2geodetic(rectifying_lat, ell: Ellipsoid = None, deg: bool = True) US Geological Survey Professional Paper 1395, US Government Printing Office, Washington, DC, 1987, pp. 13-18. """ - rectifying_lat, ell = sanitize(rectifying_lat, ell, deg) + + if deg: + rectifying_lat = radians(rectifying_lat) n = ell.thirdflattening f1 = 3 * n / 2 - 27 * n**3 / 32 @@ -447,7 +463,7 @@ def rectifying2geodetic(rectifying_lat, ell: Ellipsoid = None, deg: bool = True) # %% authalic -def geodetic2authalic(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): +def geodetic2authalic(geodetic_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from geodetic latitude to authalic latitude @@ -474,7 +490,9 @@ def geodetic2authalic(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): Office, Washington, DC, 1987, pp. 13-18. """ - geodetic_lat, ell = sanitize(geodetic_lat, ell, deg) + + if deg: + geodetic_lat = radians(geodetic_lat) e = ell.eccentricity f1 = e**2 / 3 + 31 * e**4 / 180 + 59 * e**6 / 560 @@ -491,7 +509,7 @@ def geodetic2authalic(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): return degrees(authalic_lat) if deg else authalic_lat -def authalic2geodetic(authalic_lat, ell: Ellipsoid = None, deg: bool = True): +def authalic2geodetic(authalic_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from authalic latitude to geodetic latitude @@ -517,7 +535,10 @@ def authalic2geodetic(authalic_lat, ell: Ellipsoid = None, deg: bool = True): US Geological Survey Professional Paper 1395, US Government Printing Office, Washington, DC, 1987, pp. 13-18. """ - authalic_lat, ell = sanitize(authalic_lat, ell, deg) + + if deg: + authalic_lat = radians(authalic_lat) + e = ell.eccentricity f1 = e**2 / 3 + 31 * e**4 / 180 + 517 * e**6 / 5040 f2 = 23 * e**4 / 360 + 251 * e**6 / 3780 @@ -534,7 +555,7 @@ def authalic2geodetic(authalic_lat, ell: Ellipsoid = None, deg: bool = True): # %% parametric -def geodetic2parametric(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): +def geodetic2parametric(geodetic_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from geodetic latitude to parametric latitude @@ -561,14 +582,16 @@ def geodetic2parametric(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): Office, Washington, DC, 1987, pp. 13-18. """ - geodetic_lat, ell = sanitize(geodetic_lat, ell, deg) + + if deg: + geodetic_lat = radians(geodetic_lat) parametric_lat = atan(sqrt(1 - (ell.eccentricity) ** 2) * tan(geodetic_lat)) return degrees(parametric_lat) if deg else parametric_lat -def parametric2geodetic(parametric_lat, ell: Ellipsoid = None, deg: bool = True): +def parametric2geodetic(parametric_lat, ell: Ellipsoid = ELL, deg: bool = True): """ converts from parametric latitude to geodetic latitude @@ -594,7 +617,9 @@ def parametric2geodetic(parametric_lat, ell: Ellipsoid = None, deg: bool = True) US Geological Survey Professional Paper 1395, US Government Printing Office, Washington, DC, 1987, pp. 13-18. """ - parametric_lat, ell = sanitize(parametric_lat, ell, deg) + + if deg: + parametric_lat = radians(parametric_lat) geodetic_lat = atan(tan(parametric_lat) / sqrt(1 - (ell.eccentricity) ** 2)) diff --git a/src/pymap3d/los.py b/src/pymap3d/los.py index d4e05c7..46ba9cb 100644 --- a/src/pymap3d/los.py +++ b/src/pymap3d/los.py @@ -16,6 +16,8 @@ __all__ = ["lookAtSpheroid"] +ELL = Ellipsoid.from_name("wgs84") + def lookAtSpheroid( lat0, @@ -23,7 +25,7 @@ def lookAtSpheroid( h0, az, tilt, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -63,7 +65,7 @@ def lookAtSpheroid( """ if ell is None: - ell = Ellipsoid.from_name("wgs84") + ell = ELL try: lat0 = asarray(lat0) diff --git a/src/pymap3d/lox.py b/src/pymap3d/lox.py index 2f578b1..6699d20 100644 --- a/src/pymap3d/lox.py +++ b/src/pymap3d/lox.py @@ -30,11 +30,11 @@ "meanm", ] - +ELL = Ellipsoid.from_name("wgs84") COS_EPS = 1e-9 -def meridian_dist(lat, ell: Ellipsoid = None, deg: bool = True) -> float: +def meridian_dist(lat, ell: Ellipsoid = ELL, deg: bool = True) -> float: """ Computes the ground distance on an ellipsoid from the equator to the input latitude. @@ -55,7 +55,7 @@ def meridian_dist(lat, ell: Ellipsoid = None, deg: bool = True) -> float: return meridian_arc(0.0, lat, ell, deg) -def meridian_arc(lat1, lat2, ell: Ellipsoid = None, deg: bool = True) -> float: +def meridian_arc(lat1, lat2, ell: Ellipsoid = ELL, deg: bool = True) -> float: """ Computes the ground distance on an ellipsoid between two latitudes. @@ -88,7 +88,7 @@ def loxodrome_inverse( lon1, lat2, lon2, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple[float, float]: """ @@ -180,7 +180,7 @@ def loxodrome_direct( lon1, rng, a12, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -261,7 +261,7 @@ def loxodrome_direct( return lat2, lon2 -def departure(lon1, lon2, lat, ell: Ellipsoid = None, deg: bool = True) -> float: +def departure(lon1, lon2, lat, ell: Ellipsoid = ELL, deg: bool = True) -> float: """ Computes the distance along a specific parallel between two meridians. @@ -289,7 +289,7 @@ def departure(lon1, lon2, lat, ell: Ellipsoid = None, deg: bool = True) -> float return rcurve.parallel(lat, ell=ell, deg=False) * (abs(lon2 - lon1) % pi) -def meanm(lat, lon, ell: Ellipsoid = None, deg: bool = True) -> tuple: +def meanm(lat, lon, ell: Ellipsoid = ELL, deg: bool = True) -> tuple: """ Computes geographic mean for geographic points on an ellipsoid diff --git a/src/pymap3d/ned.py b/src/pymap3d/ned.py index c888576..7c90fef 100644 --- a/src/pymap3d/ned.py +++ b/src/pymap3d/ned.py @@ -6,6 +6,18 @@ from .ellipsoid import Ellipsoid from .enu import aer2enu, enu2aer, geodetic2enu +__all__ = [ + "aer2ned", + "ned2aer", + "ned2geodetic", + "ned2ecef", + "ecef2ned", + "geodetic2ned", + "ecef2nedv", +] + +ELL = Ellipsoid.from_name("wgs84") + def aer2ned(az, elev, slantRange, deg: bool = True) -> tuple: """ @@ -73,7 +85,7 @@ def ned2geodetic( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -122,7 +134,7 @@ def ned2ecef( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -168,7 +180,7 @@ def ecef2ned( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -217,7 +229,7 @@ def geodetic2ned( lat0, lon0, h0, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ diff --git a/src/pymap3d/rcurve.py b/src/pymap3d/rcurve.py index d3ec1cd..bb62f1a 100644 --- a/src/pymap3d/rcurve.py +++ b/src/pymap3d/rcurve.py @@ -3,19 +3,22 @@ from __future__ import annotations from .ellipsoid import Ellipsoid -from .mathfun import cos, sin, sqrt -from .utils import sanitize +from .mathfun import cos, sin, sqrt, radians __all__ = ["parallel", "meridian", "transverse", "geocentric_radius"] +ELL = Ellipsoid.from_name("wgs84") -def geocentric_radius(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): + +def geocentric_radius(geodetic_lat, ell: Ellipsoid = ELL, deg: bool = True): """ compute geocentric radius at geodetic latitude https://en.wikipedia.org/wiki/Earth_radius#Geocentric_radius """ - geodetic_lat, ell = sanitize(geodetic_lat, ell, deg) + + if deg: + geodetic_lat = radians(geodetic_lat) return sqrt( ( @@ -29,7 +32,7 @@ def geocentric_radius(geodetic_lat, ell: Ellipsoid = None, deg: bool = True): ) -def parallel(lat, ell: Ellipsoid = None, deg: bool = True) -> float: +def parallel(lat, ell: Ellipsoid = ELL, deg: bool = True) -> float: """ computes the radius of the small circle encompassing the globe at the specified latitude @@ -50,12 +53,13 @@ def parallel(lat, ell: Ellipsoid = None, deg: bool = True) -> float: radius of ellipsoid (meters) """ - lat, ell = sanitize(lat, ell, deg) + if deg: + lat = radians(lat) return cos(lat) * transverse(lat, ell, deg=False) -def meridian(lat, ell: Ellipsoid = None, deg: bool = True): +def meridian(lat, ell: Ellipsoid = ELL, deg: bool = True): """computes the meridional radius of curvature for the ellipsoid like Matlab rcurve('meridian', ...) @@ -75,14 +79,15 @@ def meridian(lat, ell: Ellipsoid = None, deg: bool = True): radius of ellipsoid """ - lat, ell = sanitize(lat, ell, deg) + if deg: + lat = radians(lat) f1 = ell.semimajor_axis * (1 - ell.eccentricity**2) f2 = 1 - (ell.eccentricity * sin(lat)) ** 2 return f1 / sqrt(f2**3) -def transverse(lat, ell: Ellipsoid = None, deg: bool = True): +def transverse(lat, ell: Ellipsoid = ELL, deg: bool = True): """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 @@ -104,6 +109,7 @@ def transverse(lat, ell: Ellipsoid = None, deg: bool = True): radius of ellipsoid (meters) """ - lat, ell = sanitize(lat, ell, deg) + if deg: + lat = radians(lat) return ell.semimajor_axis / sqrt(1 - (ell.eccentricity * sin(lat)) ** 2) diff --git a/src/pymap3d/rsphere.py b/src/pymap3d/rsphere.py index 3f610c3..00424df 100644 --- a/src/pymap3d/rsphere.py +++ b/src/pymap3d/rsphere.py @@ -22,8 +22,10 @@ "biaxial", ] +ELL = Ellipsoid.from_name("wgs84") -def eqavol(ell: Ellipsoid = None) -> float: + +def eqavol(ell: Ellipsoid = ELL) -> float: """computes the radius of the sphere with equal volume as the ellipsoid Parameters @@ -36,15 +38,13 @@ def eqavol(ell: Ellipsoid = None) -> float: radius: float radius of sphere """ - if ell is None: - ell = Ellipsoid.from_name("wgs84") f = ell.flattening return ell.semimajor_axis * (1 - f / 3 - f**2 / 9) -def authalic(ell: Ellipsoid = None) -> float: +def authalic(ell: Ellipsoid = ELL) -> float: """computes the radius of the sphere with equal surface area as the ellipsoid Parameters @@ -57,8 +57,6 @@ def authalic(ell: Ellipsoid = None) -> float: radius: float radius of sphere """ - if ell is None: - ell = Ellipsoid.from_name("wgs84") e = ell.eccentricity @@ -71,7 +69,7 @@ def authalic(ell: Ellipsoid = None) -> float: return ell.semimajor_axis -def rectifying(ell: Ellipsoid = None) -> float: +def rectifying(ell: Ellipsoid = ELL) -> float: """computes the radius of the sphere with equal meridional distances as the ellipsoid Parameters @@ -84,8 +82,7 @@ def rectifying(ell: Ellipsoid = None) -> float: radius: float radius of sphere """ - if ell is None: - ell = Ellipsoid.from_name("wgs84") + return ((ell.semimajor_axis ** (3 / 2) + ell.semiminor_axis ** (3 / 2)) / 2) ** (2 / 3) @@ -94,7 +91,7 @@ def euler( lon1, lat2, lon2, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ): """computes the Euler radii of curvature at the midpoint of the @@ -140,7 +137,7 @@ def euler( return rho * nu / den -def curve(lat, ell: Ellipsoid = None, deg: bool = True, method: str = "mean"): +def curve(lat, ell: Ellipsoid = ELL, deg: bool = True, method: str = "mean"): """computes the arithmetic average of the transverse and meridional radii of curvature at a specified latitude point @@ -175,7 +172,7 @@ def curve(lat, ell: Ellipsoid = None, deg: bool = True, method: str = "mean"): raise ValueError("method must be mean or norm") -def triaxial(ell: Ellipsoid = None, method: str = "mean") -> float: +def triaxial(ell: Ellipsoid = ELL, method: str = "mean") -> float: """computes triaxial average of the semimajor and semiminor axes of the ellipsoid Parameters @@ -191,9 +188,6 @@ def triaxial(ell: Ellipsoid = None, method: str = "mean") -> float: radius of sphere """ - if ell is None: - ell = Ellipsoid.from_name("wgs84") - if method == "mean": return (2 * ell.semimajor_axis + ell.semiminor_axis) / 3 elif method == "norm": @@ -202,7 +196,7 @@ def triaxial(ell: Ellipsoid = None, method: str = "mean") -> float: raise ValueError("method must be mean or norm") -def biaxial(ell: Ellipsoid = None, method: str = "mean") -> float: +def biaxial(ell: Ellipsoid = ELL, method: str = "mean") -> float: """computes biaxial average of the semimajor and semiminor axes of the ellipsoid Parameters @@ -218,9 +212,6 @@ def biaxial(ell: Ellipsoid = None, method: str = "mean") -> float: radius of sphere """ - if ell is None: - ell = Ellipsoid.from_name("wgs84") - if method == "mean": return (ell.semimajor_axis + ell.semiminor_axis) / 2 elif method == "norm": diff --git a/src/pymap3d/spherical.py b/src/pymap3d/spherical.py index f6662c0..ecfc4ed 100644 --- a/src/pymap3d/spherical.py +++ b/src/pymap3d/spherical.py @@ -7,19 +7,20 @@ from .ellipsoid import Ellipsoid from .mathfun import asin, atan2, cbrt, degrees, hypot, power, radians, sin, sqrt -from .utils import sanitize __all__ = [ "geodetic2spherical", "spherical2geodetic", ] +ELL = Ellipsoid.from_name("wgs84") + def geodetic2spherical( lat, lon, alt, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -58,8 +59,9 @@ def geodetic2spherical( geodetic coordinates. Journal of Geodesy. 76. 451-454. doi:10.1007/s00190-002-0273-6 """ - lat, ell = sanitize(lat, ell, deg) + if deg: + lat = radians(lat) lon = radians(lon) # Pre-compute to avoid repeated trigonometric functions @@ -90,7 +92,7 @@ def spherical2geodetic( lat, lon, radius, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, deg: bool = True, ) -> tuple: """ @@ -124,8 +126,9 @@ def spherical2geodetic( geodetic coordinates. Journal of Geodesy. 76. 451-454. doi:10.1007/s00190-002-0273-6 """ - lat, ell = sanitize(lat, ell, deg) + if deg: + lat = radians(lat) lon = radians(lon) # Pre-compute to avoid repeated trigonometric functions diff --git a/src/pymap3d/tests/test_geodetic.py b/src/pymap3d/tests/test_geodetic.py index 4b9e456..38dc574 100755 --- a/src/pymap3d/tests/test_geodetic.py +++ b/src/pymap3d/tests/test_geodetic.py @@ -185,9 +185,6 @@ def test_ecef(): assert y == approx(xyz[1]) assert z == approx(xyz[2]) - with pytest.raises(ValueError): - pm.geodetic2ecef(-100, lla0[1], lla0[2]) - assert pm.ecef2geodetic(*xyz) == approx(lla0) assert pm.ecef2geodetic(*xyz, deg=False) == approx(rlla0) diff --git a/src/pymap3d/tests/test_latitude.py b/src/pymap3d/tests/test_latitude.py index fc6e502..058c8b4 100644 --- a/src/pymap3d/tests/test_latitude.py +++ b/src/pymap3d/tests/test_latitude.py @@ -167,28 +167,3 @@ def test_numpy_geodetic_parametric(): pytest.importorskip("numpy") assert latitude.geodetic2parametric([45, 0]) == approx([44.9037878, 0]) assert latitude.parametric2geodetic([44.9037878, 0]) == approx([45, 0]) - - -@pytest.mark.parametrize("lat", [91, -91]) -def test_badvals(lat): - # geodetic_isometric is not included on purpose - with pytest.raises(ValueError): - latitude.geodetic2geocentric(lat, 0) - with pytest.raises(ValueError): - latitude.geocentric2geodetic(lat, 0) - with pytest.raises(ValueError): - latitude.geodetic2conformal(lat) - with pytest.raises(ValueError): - latitude.conformal2geodetic(lat) - with pytest.raises(ValueError): - latitude.geodetic2rectifying(lat) - with pytest.raises(ValueError): - latitude.rectifying2geodetic(lat) - with pytest.raises(ValueError): - latitude.geodetic2authalic(lat) - with pytest.raises(ValueError): - latitude.authalic2geodetic(lat) - with pytest.raises(ValueError): - latitude.geodetic2parametric(lat) - with pytest.raises(ValueError): - latitude.parametric2geodetic(lat) diff --git a/src/pymap3d/tests/test_rsphere.py b/src/pymap3d/tests/test_rsphere.py index 2ad4b62..461b56e 100644 --- a/src/pymap3d/tests/test_rsphere.py +++ b/src/pymap3d/tests/test_rsphere.py @@ -15,12 +15,6 @@ def test_geocentric_radius(): assert rcurve.geocentric_radius(30) == approx(6372824.0) -@pytest.mark.parametrize("bad_lat", [-91, 91]) -def test_geocentric_radius_badval(bad_lat): - with pytest.raises(ValueError): - rcurve.geocentric_radius(bad_lat) - - def test_rsphere_eqavol(): assert rsphere.eqavol() == approx(6371000.8049) diff --git a/src/pymap3d/utils.py b/src/pymap3d/utils.py index b70c4e3..334cef2 100644 --- a/src/pymap3d/utils.py +++ b/src/pymap3d/utils.py @@ -4,15 +4,7 @@ from __future__ import annotations -from math import pi - -try: - from numpy import asarray -except ImportError: - pass - -from .ellipsoid import Ellipsoid -from .mathfun import atan2, cos, hypot, radians, sin +from .mathfun import atan2, cos, hypot, sin __all__ = ["cart2pol", "pol2cart", "cart2sph", "sph2cart"] @@ -43,25 +35,3 @@ def sph2cart(az, el, r) -> tuple: y = rcos_theta * sin(az) z = r * sin(el) return x, y, z - - -def sanitize(lat, ell: Ellipsoid | None, deg: bool) -> tuple: - if ell is None: - ell = Ellipsoid.from_name("wgs84") - - try: - lat = asarray(lat) - except NameError: - pass - - if deg: - lat = radians(lat) - - try: - if (abs(lat) > pi / 2).any(): - raise ValueError("-pi/2 <= latitude <= pi/2") - except AttributeError: - if abs(lat) > pi / 2: - raise ValueError("-pi/2 <= latitude <= pi/2") - - return lat, ell diff --git a/src/pymap3d/vincenty.py b/src/pymap3d/vincenty.py index 15b76f2..d4b8fd1 100644 --- a/src/pymap3d/vincenty.py +++ b/src/pymap3d/vincenty.py @@ -33,8 +33,10 @@ __all__ = ["vdist", "vreckon", "track2"] +ELL = Ellipsoid.from_name("wgs84") -def vdist(Lat1, Lon1, Lat2, Lon2, ell: Ellipsoid = None, deg: bool = True) -> tuple: + +def vdist(Lat1, Lon1, Lat2, Lon2, ell: Ellipsoid = ELL, deg: bool = True) -> tuple: """ Using the reference ellipsoid, compute the distance between two points within a few millimeters of accuracy, compute forward azimuth, @@ -106,9 +108,6 @@ def vdist(Lat1, Lon1, Lat2, Lon2, ell: Ellipsoid = None, deg: bool = True) -> tu 12. No warranties; use at your own risk. """ - if ell is None: - ell = Ellipsoid.from_name("wgs84") - # %% Input check: try: Lat1 = atleast_1d(Lat1) @@ -284,7 +283,7 @@ def vdist(Lat1, Lon1, Lat2, Lon2, ell: Ellipsoid = None, deg: bool = True) -> tu return dist_m, a12 -def vreckon(Lat1, Lon1, Rng, Azim, ell: Ellipsoid = None, deg: bool = True) -> tuple: +def vreckon(Lat1, Lon1, Rng, Azim, ell: Ellipsoid = ELL, deg: bool = True) -> tuple: """ This is the Vincenty "forward" solution. @@ -355,14 +354,9 @@ def vreckon(Lat1, Lon1, Rng, Azim, ell: Ellipsoid = None, deg: bool = True) -> t if Rng < 0.0: raise ValueError("Ground distance must be positive") - if ell is not None: - a = ell.semimajor_axis - b = ell.semiminor_axis - f = ell.flattening - else: # Supply WGS84 earth ellipsoid axis lengths in meters: - a = 6378137 # semimajor axis - b = 6356752.31424518 # WGS84 earth flattening coefficient definition - f = (a - b) / a + a = ell.semimajor_axis + b = ell.semiminor_axis + f = ell.flattening if deg: Lat1 = radians(Lat1) # intial latitude in radians @@ -480,7 +474,7 @@ def track2( lon1, lat2, lon2, - ell: Ellipsoid = None, + ell: Ellipsoid = ELL, npts: int = 100, deg: bool = True, ) -> tuple[list, list]: @@ -516,9 +510,6 @@ def track2( Based on code posted to the GMT mailing list in Dec 1999 by Jim Levens and by Jeff Whitaker """ - if ell is None: - ell = Ellipsoid.from_name("wgs84") - if npts < 2: raise ValueError("npts must be greater than 1") if npts == 2: