From 01daf11740928e01510615f55f9024d992dc4219 Mon Sep 17 00:00:00 2001 From: rpavlick Date: Sun, 11 Aug 2019 13:59:34 -0700 Subject: [PATCH 1/5] Adding rcurve.py containing several functions that compute the radii of curvature for an ellipsoid. Adding rsphere.py containing several functions that compute the radii of auxiliary spheres. Adding latitude.py containing several functions that compute auxiliary latitudes and their inverse functions. Adding 'loxodrome_direct' to compute the end point of a loxodromic curve from a start lat/lon, azimuth, and distance. Fixing 'loxodrome_inverse', it previously failed on the degenerate case of east/west azimuth. Adding 'meanm' which computes the geographic mean of points on a spheroid. Adding 'departure' which computes the distance between two longitudes at a specified latitude. Adding 'meridianarc' which computes the distance between two latitudes along a specified longitude. Added thirdflattening and eccentricty attributes to ellipsoid definitions. Removed all instances '%%' in comments. --- pymap3d/__init__.py | 6 +- pymap3d/ecef.py | 2 +- pymap3d/eci.py | 8 +- pymap3d/ellipsoid.py | 22 +- pymap3d/latitude.py | 579 +++++++++++++++++++++++++++++++++++++++++++ pymap3d/los.py | 4 +- pymap3d/lox.py | 289 ++++++++++++--------- pymap3d/rcurve.py | 47 ++++ pymap3d/rsphere.py | 109 ++++++++ pymap3d/sidereal.py | 6 +- pymap3d/util.py | 42 ++++ pymap3d/vallado.py | 8 +- pymap3d/vincenty.py | 20 +- 13 files changed, 992 insertions(+), 150 deletions(-) create mode 100644 pymap3d/latitude.py create mode 100644 pymap3d/rcurve.py create mode 100644 pymap3d/rsphere.py create mode 100644 pymap3d/util.py diff --git a/pymap3d/__init__.py b/pymap3d/__init__.py index 2b733a1..3767edd 100644 --- a/pymap3d/__init__.py +++ b/pymap3d/__init__.py @@ -52,6 +52,10 @@ from .enu import enu2geodetic, geodetic2enu, aer2enu, enu2aer # noqa: F401 from .aer import ecef2aer, aer2ecef, geodetic2aer, aer2geodetic, eci2aer, aer2eci # noqa: F401 from .los import lookAtSpheroid # noqa: F401 - from .lox import isometric, meridian_dist, loxodrome_inverse # noqa: F401 + from .lox import loxodrome_direct, loxodrome_inverse, meridianarc, departure, meanm + from .latitude import geodetic2isometric, isometric2geodetic, geodetic2rectifying, rectifying2geodetic, geodetic2conformal, conformal2geodetic, geodetic2parametric, parametric2geodetic, geodetic2geocentric, geocentric2geodetic, geodetic2authalic, authalic2geodetic + from .rcurve import rcurve_meridian, rcurve_parallel, rcurve_transverse + from .rsphere import rsphere_authalic, rsphere_biaxial, rsphere_curve, rsphere_eqavol, rsphere_euler, rsphere_rectifying, rsphere_triaxial + from .util import wrapToPi, wrapTo2Pi, sph2cart, cart2sph, pol2cart, cart2pol else: # pure Python only from .math import * # type: ignore # noqa: F401, F403 diff --git a/pymap3d/ecef.py b/pymap3d/ecef.py index 7814d7f..f7961b7 100644 --- a/pymap3d/ecef.py +++ b/pymap3d/ecef.py @@ -126,7 +126,7 @@ def ecef2geodetic(x: float, y: float, z: float, eps = ((ell.semiminor_axis * u - ell.semimajor_axis * huE + E**2) * sin(Beta)) / (ell.semimajor_axis * huE * 1 / cos(Beta) - E**2 * cos(Beta)) Beta += eps -# %% final output +# final output lat = arctan(ell.semimajor_axis / ell.semiminor_axis * tan(Beta)) lon = arctan2(y, x) diff --git a/pymap3d/eci.py b/pymap3d/eci.py index 38026fc..f8ae26b 100644 --- a/pymap3d/eci.py +++ b/pymap3d/eci.py @@ -41,14 +41,14 @@ def eci2ecef(x: float, y: float, z: float = None, z : float target z ECEF coordinate """ -# %% +# # FIXME: temporary old API, which was a single N x 3 numpy.ndarray if z is None and isinstance(y, (str, datetime)): time = y z = x[:, 2] y = x[:, 1] x = x[:, 0] -# %% +# x = np.atleast_1d(x) y = np.atleast_1d(y) z = np.atleast_1d(z) @@ -112,14 +112,14 @@ def ecef2eci(x: float, y: float, z: float = None, z : float target z ECI coordinate """ -# %% +# # FIXME: temporary old API, which was a single N x 3 numpy.ndarray if z is None and isinstance(y, (str, datetime)): time = y z = x[:, 2] y = x[:, 1] x = x[:, 0] -# %% +# x = np.atleast_1d(x) y = np.atleast_1d(y) z = np.atleast_1d(z) diff --git a/pymap3d/ellipsoid.py b/pymap3d/ellipsoid.py index 1ed4939..8fe713f 100644 --- a/pymap3d/ellipsoid.py +++ b/pymap3d/ellipsoid.py @@ -1,4 +1,4 @@ - +from math import sqrt class Ellipsoid: """ @@ -26,19 +26,27 @@ def __init__(self, model='wgs84'): self.semimajor_axis = 6378137. self.flattening = 1 / 298.2572235630 self.semiminor_axis = self.semimajor_axis * (1 - self.flattening) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'wgs72': self.semimajor_axis = 6378135. self.flattening = 1 / 298.26 self.semiminor_axis = self.semimajor_axis * (1 - self.flattening) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'grs80': """https://en.wikipedia.org/wiki/GRS_80""" self.semimajor_axis = 6378137. self.flattening = 1 / 298.257222100882711243 self.semiminor_axis = self.semimajor_axis * (1 - self.flattening) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'clarke1866': self.semimajor_axis = 6378206.4 self.semiminor_axis = 6356583.8 self.flattening = -(self.semiminor_axis / self.semimajor_axis - 1) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'mars': """ https://tharsis.gsfc.nasa.gov/geodesy.html @@ -46,18 +54,26 @@ def __init__(self, model='wgs84'): self.semimajor_axis = 3396900. self.semiminor_axis = 3376097.80585952 self.flattening = 1 / 163.295274386012 + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'moon': self.semimajor_axis = 1738000. self.semiminor_axis = self.semimajor_axis self.flattening = 0. + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'venus': self.semimajor_axis = 6051000. self.semiminor_axis = self.semimajor_axis + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) self.flattening = 0. elif model == 'jupiter': self.semimajor_axis = 71492000. self.flattening = 1 / 15.415446277169725 self.flattening = -(self.semiminor_axis / self.semimajor_axis - 1) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'io': """ https://doi.org/10.1006/icar.1998.5987 @@ -65,9 +81,13 @@ def __init__(self, model='wgs84'): self.semimajor_axis = 1829.7 self.flattening = 1 / 131.633093525179 self.semiminor_axis = self.semimajor_axis * (1 - self.flattening) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'pluto': self.semimajor_axis = 1187000. self.semiminor_axis = self.semimajor_axis self.flattening = 0. + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) else: raise NotImplementedError('{} model not implemented, let us know and we will add it (or make a pull request)'.format(model)) diff --git a/pymap3d/latitude.py b/pymap3d/latitude.py new file mode 100644 index 0000000..a2deccc --- /dev/null +++ b/pymap3d/latitude.py @@ -0,0 +1,579 @@ +""" compute auxiliary latitudes and their inverse""" +from typing import Tuple +import numpy as np +from .ellipsoid import Ellipsoid +from numpy import radians, sin, cos, tan, exp, arctan, log, hypot, degrees, arctan2, sqrt, pi +from .rsphere import rsphere_rectifying +from .rcurve import rcurve_parallel +from .util import wrapToPi, sph2cart, cart2sph + +__all__ = ['geodetic2isometric', 'isometric2geodetic', 'geodetic2rectifying', 'rectifying2geodetic', 'geodetic2conformal', 'conformal2geodetic', 'geodetic2parametric', 'parametric2geodetic', 'geodetic2geocentric', 'geocentric2geodetic', 'geodetic2authalic', 'authalic2geodetic'] + +def geodetic2isometric(lat: float, ell: Ellipsoid = None, deg: bool = True): + """ + computes isometric latitude of a point on an ellipsoid + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + isolat : float or numpy.ndarray of float + isometric latiude + + Notes + ----- + + Isometric latitude is an auxiliary latitude proportional to the spacing + of parallels of latitude on an ellipsoidal mercator projection. + + Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, + School of Mathematical and Geospatial Sciences, RMIT University, + January 2010 + """ + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = np.deg2rad(lat) + + x = ell.eccentricity * np.sin(lat) + y = (1 - x)/(1 + x) + z = np.pi/4 + lat/2 + +# calculate the isometric latitude + isolat = np.log(np.tan(z) * (y**(ell.eccentricity /2))) + + if deg is True: + isolat = np.degrees(isolat) + + return isolat + +def isometric2geodetic(isolat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from isometric latitude to geodetic latitude + + Parameters + ---------- + + isolat : float or numpy.ndarray of float + isometric latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + + + if ell is None: + ell = Ellipsoid() + + n = ell.thirdflattening + f1 = 3*n/2 - 27*n**3/32 + f2 = 21*n**2/16 - 55*n**4/32 + f3 = 151*n**3/96 + f4 = 1097*n**4/512 + + if deg is True: + isolat = radians(isolat) + + cnflat = 2 * arctan(exp(isolat)) - (pi/2) + lat = conformal2geodetic(cnflat, ell, deg=False) + + if deg is True: + lat = degrees(lat) + + return lat + +def geodetic2rectifying(lat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from geodetic latitude to rectifying latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + reclat : float or numpy.ndarray of float + rectifying latiude + + Notes + ----- + + Rectifying latitude is an auxiliary latitude used to map an ellipsoid to a + sphere such that correct distances along meridians are preserved. + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + n = ell.thirdflattening + f1 = 3*n/2 - 9*n**3/16 + f2 = 15*n**2/16 - 15*n**4/32 + f3 = 35*n**3/48 + f4 = 315*n**4/512 + + if deg is True: + lat = radians(lat) + + reclat = lat - f1*sin(2*lat) + f2*sin(4*lat) - f3*sin(6*lat) + f4*sin(8*lat) + + if deg is True: + reclat = degrees(reclat) + + return reclat + + +def rectifying2geodetic(reclat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from rectifying latitude to geodetic latitude + + Parameters + ---------- + + reclat : float or numpy.ndarray of float + rectifying latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + n = ell.thirdflattening + f1 = 3*n/2 - 27*n**3/32 + f2 = 21*n**2/16 - 55*n**4/32 + f3 = 151*n**3/96 + f4 = 1097*n**4/512 + + if deg is True: + reclat = radians(reclat) + + lat = reclat + f1*sin(2*reclat) + f2*sin(4*reclat) + f3*sin(6*reclat) + f4*sin(8*reclat) + + if deg is True: + lat = degrees(lat) + + return lat + + +def geodetic2conformal(lat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from geodetic latitude to conformal latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + cnflat : float or numpy.ndarray of float + conformal latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + e = ell.eccentricity + f1 = 1 - e*sin(lat) + f2 = 1 + e*sin(lat) + f3 = 1 - sin(lat) + f4 = 1 + sin(lat) + + # compute conformal latitudes with correction for points at +90 + cnflat = np.where(f3 == 0, pi/2, 2 * arctan(sqrt((f4/f3) * ((f1/f2)**e))) - (pi/2)) + + if deg is True: + cnflat = degrees(cnflat) + + return cnflat + +def conformal2geodetic(cnflat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from conformal latitude to geodetic latitude + + Parameters + ---------- + + cnflat : float or numpy.ndarray of float + conformal latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + cnflat = radians(cnflat) + + e = ell.eccentricity + f1 = e**2/2 + 5*e**4/24 + e**6/12 + 13*e**8/360 + f2 = 7*e**4/48 + 29*e**6/240 + 811*e**8/11520 + f3 = 7*e**6/120 + 81*e**8/1120 + f4 = 4279*e**8/161280 + + lat = cnflat + f1*sin(2*cnflat) + f2*sin(4*cnflat) + f3*sin(6*cnflat) + f4*sin(8*cnflat) + + if deg is True: + lat = degrees(lat) + + return lat + +def geodetic2parametric(lat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from geodetic latitude to parametric latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + parlat : float or numpy.ndarray of float + parametric latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + epsilon = 1E-10 + + lat = np.where(lat == pi/2, pi/2 - epsilon, lat) + lat = np.where(lat == -pi/2, -pi/2 + epsilon, lat) + + parlat = arctan(sqrt(1-(ell.eccentricity)**2) * tan(lat)) + + if deg is True: + parlat = degrees(parlat) + + return parlat + +def parametric2geodetic(parlat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from parametric latitude to geodetic latitude + + Parameters + ---------- + + parlat : float or numpy.ndarray of float + parametric latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + parlat = radians(parlat) + + epsilon = 1E-10 + + parlat = np.where(parlat == pi/2, pi/2 - epsilon, parlat) + parlat = np.where(parlat == -pi/2, -pi/2 + epsilon, parlat) + + lat = arctan(tan(parlat)/ sqrt(1-(ell.eccentricity)**2)) + + if deg is True: + lat = degrees(lat) + + return lat + +def geodetic2geocentric(lat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from geodetic latitude to geocentric latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + cenlat : float or numpy.ndarray of float + geocentric latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + epsilon = 1E-10 + + lat = np.where(lat == pi/2, pi/2 - epsilon, lat) + lat = np.where(lat == -pi/2, -pi/2 + epsilon, lat) + + cenlat = arctan((1-(ell.eccentricity)**2) * tan(lat)) + + if deg is True: + cenlat = degrees(cenlat) + + return cenlat + +def geocentric2geodetic(cenlat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from geocentric latitude to geodetic latitude + + Parameters + ---------- + + cenlat : float or numpy.ndarray of float + geocentric latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + cenlat = radians(cenlat) + + epsilon = 1E-10 + + cenlat = np.where(cenlat == pi/2, pi/2 - epsilon, cenlat) + cenlat = np.where(cenlat == -pi/2, -pi/2 + epsilon, cenlat) + + lat = arctan(tan(cenlat)/ (1-(ell.eccentricity)**2)) + + if deg is True: + lat = degrees(lat) + + return lat + +def geodetic2authalic(lat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from geodetic latitude to authalic latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + autlat : float or numpy.ndarray of float + authalic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + e = ell.eccentricity + f1 = e**2/3 + 31*e**4/180 + 59*e**6/560 + f2 = 17*e**4/360 + 61*e**6/1260 + f3 = 383*e**6/45360 + + autlat = lat - f1*sin(2*lat) + f2*sin(4*lat) - f3*sin(6*lat) + + if deg is True: + autlat = degrees(autlat) + + return autlat + +def authalic2geodetic(autlat, ell: Ellipsoid = None, deg: bool = True): + ''' + converts from authalic latitude to geodetic latitude + + Parameters + ---------- + + autlat : float or numpy.ndarray of float + authalic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + autlat = radians(autlat) + + e = ell.eccentricity + f1 = e**2/3 + 31*e**4/180 + 517*e**6/5040 + f2 = 23*e**4/360 + 251*e**6/3780 + f3 = 761*e**6/45360 + + lat = autlat + f1*sin(2*autlat) + f2*sin(4*autlat) + f3*sin(6*autlat) + + if deg is True: + lat = degrees(lat) + + return lat \ No newline at end of file diff --git a/pymap3d/los.py b/pymap3d/los.py index ef175fb..bfe47a0 100644 --- a/pymap3d/los.py +++ b/pymap3d/los.py @@ -72,13 +72,13 @@ def lookAtSpheroid(lat0: float, lon0: float, h0: float, az: float, tilt: float, 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 +# Return nan if radical < 0 or d < 0 because LOS vector does not point towards Earth with np.errstate(invalid='ignore'): d = np.where(radical > 0, (value - a * b * c * np.sqrt(radical)) / magnitude, np.nan) d[d < 0] = np.nan -# %% cartesian to ellipsodal +# cartesian to ellipsodal lat, lon, _ = ecef2geodetic(x + d * u, y + d * v, z + d * w, deg=deg) return lat, lon, d diff --git a/pymap3d/lox.py b/pymap3d/lox.py index 05671af..acd469e 100644 --- a/pymap3d/lox.py +++ b/pymap3d/lox.py @@ -1,149 +1,192 @@ -""" isometric latitude, meridian distance """ import numpy as np from .ellipsoid import Ellipsoid +from typing import Tuple +from .latitude import geodetic2isometric, geodetic2rectifying, rectifying2geodetic, geodetic2authalic, authalic2geodetic +from .rsphere import rsphere_rectifying +from .rcurve import rcurve_parallel +from .util import wrapToPi, sph2cart, cart2sph +from numpy import arctan, sqrt, tan, sign, sin, cos, arctan2, arcsin, nan, pi, radians, degrees +__all__ = ['loxodrome_inverse', 'loxodrome_direct', 'medianarc', 'departure', 'meanm'] -def isometric(lat: float, ell: Ellipsoid = None, deg: bool = True): +def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float, + ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float]: """ - computes isometric latitude of a point on an ellipsoid + computes the arc length and azimuth of the loxodrome + between two points on the surface of the reference ellipsoid - Parameters - ---------- + Parameters + ---------- - lat : float or numpy.ndarray of float - geodetic latitude - ell : Ellipsoid, optional + lat1 : float or numpy.ndarray of float + geodetic latitude of first point + lon1 : float or numpy.ndarray of float + geodetic longitude of first point + lat2 : float or numpy.ndarray of float + geodetic latitude of second point + lon2 : float or numpy.ndarray of float + geodetic longitude of second point + ell : Ellipsoid, optional reference ellipsoid (default WGS84) - deg : bool, optional + deg : bool, optional degrees input/output (False: radians in/out) - Returns - ------- - - isolat : float or numpy.ndarray of float - isometric latiude + Results + ------- - Notes - ----- + dist : float or numpy.ndarray of float + distance along loxodrome + a12 : float or numpy.ndarray of float + azimuth of loxodrome (degrees/radians) - Isometric latitude is an auxiliary latitude proportional to the spacing - of parallels of latitude on an ellipsoidal mercator projection. + Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, + School of Mathematical and Geospatial Sciences, RMIT University, January 2010 - Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, - School of Mathematical and Geospatial Sciences, RMIT University, - January 2010 + [1] Bowring, B.R., 1985, 'The geometry of the loxodrome on the + ellipsoid', The Canadian Surveyor, Vol. 39, No. 3, Autumn 1985, + pp.223-230. + [2] Snyder, J.P., 1987, Map Projections-A Working Manual. U.S. + Geological Survey Professional Paper 1395. Washington, DC: U.S. + Government Printing Office, pp.15-16 and pp. 44-45. + [3] Thomas, P.D., 1952, Conformal Projections in Geodesy and + Cartography, Special Publication No. 251, Coast and Geodetic + Survey, U.S. Department of Commerce, Washington, DC: U.S. + Government Printing Office, p. 66. """ + # set ellipsoid parameters if ell is None: ell = Ellipsoid() - f = ell.flattening # flattening of ellipsoid - if deg is True: - lat = np.deg2rad(lat) + lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2]) - e2 = f * (2 - f) # eccentricity-squared - e = np.sqrt(e2) # eccentricity of ellipsoid + # compute isometric latitude of P1 and P2 + isolat1 = geodetic2isometric(lat1, deg=False, ell=ell) + isolat2 = geodetic2isometric(lat2, deg=False, ell=ell) - x = e * np.sin(lat) - y = (1 - x) / (1 + x) - z = np.pi / 4 + lat / 2 + # compute changes in isometric latitude and longitude between points + disolat = isolat2 - isolat1 -# calculate the isometric latitude - isolat = np.log(np.tan(z) * (y**(e / 2))) + dlon = abs(wrapToPi(lon2-lon1)) + + # compute azimuth + a12 = np.arctan2(dlon, disolat) + cosaz = abs(cos(a12)) + + # compute distance along loxodromic curve + dist = meridianarc(lat2, lat1, deg=False, ell=ell) / cosaz + + # consider degenerate case (directly east/west) + epsilon = 1e-10 + dist = np.where(cosaz < epsilon,departure(lon1, lon2, (lat1+lat2)/2, ell, deg=False), dist) if deg is True: - isolat = np.degrees(isolat) + a12 = np.degrees(a12) % 360. - return isolat + return dist, a12 -def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True): +def loxodrome_direct(lat1: float, lon1: float, rng: float, a12: float, + ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]: """ - computes the meridian distance on an ellipsoid *from the equator* to a latitude. Parameters ---------- - lat : float or numpy.ndarray of float - geodetic latitude + lat1 : float + inital geodetic latitude (degrees) + lon1 : float + initial geodetic longitude (degrees) + rng : float + ground distance (meters) + a12 : float + azimuth (degrees) clockwide from north. ell : Ellipsoid, optional - reference ellipsoid (default WGS84) - deg : bool, optional - degrees input/output (False: radians in/out) + reference ellipsoid Results ------- - mdist : float or numpy.ndarray of float - meridian distance (degrees/radians) + lat2 : float + final geodetic latitude (degrees) + lon2 : float + final geodetic longitude (degrees) + a21 : float + reverse azimuth (degrees), at final point facing back toward the intial point +""" + if ell is None: + ell = Ellipsoid() - Notes - ----- + lat1 = np.atleast_1d(lat1) + lon1 = np.atleast_1d(lon1) + rng = np.atleast_1d(rng) + a12 = np.atleast_1d(a12) - Formula given Baeschlin, C.F., 1948, - "Lehrbuch Der Geodasie", Orell Fussli Verlag, Zurich, pp.47-50. + if rng.ndim != 1 or a12.ndim != 1: + raise ValueError('Range and azimuth must be scalar or vector') - Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, - School of Mathematical and Geospatial Sciences, RMIT University, January 2010 - """ + if lat1.size > 1 and rng.size > 1: + raise ValueError('LOXODROME_DIRECT: Variable ranges are only allowed for a single point.') + + if rng.size != a12.size and rng.size == 1: + rng = np.broadcast_to(rng, a12.size) if deg is True: - lat = np.radians(lat) + lat1, lon1, a12 = radians([lat1, lon1, a12]) + + if abs(lat1) > (pi/2): + raise ValueError('LOXODROME_DIRECT: Input lat. must be between -90 and 90 deg., inclusive.') - # set ellipsoid parameters - if ell is None: - ell = Ellipsoid() + # if range is negative, use inverse bearing + a12 = np.where(sign(rng) == -1, (a12 + pi) % (2*pi), a12) + rng = abs(rng) - a = ell.semimajor_axis - f = ell.flattening # flattening of ellipsoid + # compute rectifying sphere latitude and radius + reclat = geodetic2rectifying(lat1, ell, deg=False) - e2 = f * (2 - f) # eccentricity-squared + lat2 = np.zeros_like(lat1) + lon2 = np.zeros_like(lon1) + epsilon = 1e-10 # Set tolerance (should set this for specific machine) - # powers of eccentricity - e4 = e2 * e2 - e6 = e4 * e2 - e8 = e6 * e2 - e10 = e8 * e2 + # correct azimuths at either pole. + a12 = np.where(lat1 >= (pi/2) - epsilon, pi, a12) + a12 = np.where(lat1 <= (-pi/2) + epsilon, 0, a12) - # coefficients of series expansion for meridian distance - A = 1 + (3 / 4) * e2 + (45 / 64) * e4 + (175 / 256) * e6 + (11025 / 16384) * e8 + (43659 / 65536) * e10 - B = (3 / 4) * e2 + (15 / 16) * e4 + (525 / 512) * e6 + (2205 / 2048) * e8 + (72765 / 65536) * e10 - C = (15 / 64) * e4 + (105 / 256) * e6 + (2205 / 4096) * e8 + (10395 / 16384) * e10 - D = (35 / 512) * e6 + (315 / 2048) * e8 + (31185 / 131072) * e10 - E = (315 / 16384) * e8 + (3465 / 65536) * e10 - F = (693 / 131072) * e10 + # compute the new points + cosaz = cos(a12) + lat2 = reclat + (rng/rsphere_rectifying(ell))*cosaz # compute rectifying latitude + lat2 = rectifying2geodetic(lat2,ell,deg=False) # transform to geodetic latitude - term1 = A * lat - term2 = (B / 2) * np.sin(2 * lat) - term3 = (C / 4) * np.sin(4 * lat) - term4 = (D / 6) * np.sin(6 * lat) - term5 = (E / 8) * np.sin(8 * lat) - term6 = (F / 10) * np.sin(10 * lat) + # nudge latitudes near either pole + lat2 = np.where(lat2 >= (pi/2) - epsilon, (pi/2) - epsilon, lat2) + lat2 = np.where(lat2 <= (-pi/2) + epsilon, (-pi/2) + epsilon , lat2) + lat2 = np.where(abs(cosaz)<=epsilon, lat1, lat2) - mdist = a * (1 - e2) * (term1 - term2 + term3 - term4 + term5 - term6) + newiso = geodetic2isometric(lat2,ell,deg=False) + iso = geodetic2isometric(lat1,ell,deg=False) + + dlon = np.where(abs(cosaz)<=epsilon, sign(sin(a12)) * rng/rcurve_parallel(lat1, ell, deg=False), tan(a12) * (newiso - iso)) + lon2 = lon1 + dlon - return mdist + lon2 = wrapToPi(lon2) + a21 = (a12 + pi)%(2*pi) + if deg is True: + lat2, lon2, a21 = degrees([lat2, lon2, a21]) -def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float, - ell: Ellipsoid = None, deg: bool = True): - """ - computes the arc length and azimuth of the loxodrome - between two points on the surface of the reference ellipsoid + return lat2, lon2, a21 + +def meridianarc(lat1: float, lat2: float, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + Computes the meridian distance on an ellipsoid between two latitudes. Parameters ---------- - lat1 : float or numpy.ndarray of float - geodetic latitude of first point - lon1 : float or numpy.ndarray of float - geodetic longitude of first point - lat2 : float or numpy.ndarray of float - geodetic latitude of second point - lon2 : float or numpy.ndarray of float - geodetic longitude of second point + lat1, lat2 : float or numpy.ndarray of float + geodetic latitudes ell : Ellipsoid, optional reference ellipsoid (default WGS84) deg : bool, optional @@ -152,51 +195,49 @@ def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float, Results ------- - lox_s : float or numpy.ndarray of float - distance along loxodrome - az12 : float or numpy.ndarray of float - azimuth of loxodrome (degrees/radians) - - Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, - School of Mathematical and Geospatial Sciences, RMIT University, January 2010 + dist : float or numpy.ndarray of float + distance (units same as ellipsoid) + ''' - [1] Bowring, B.R., 1985, 'The geometry of the loxodrome on the - ellipsoid', The Canadian Surveyor, Vol. 39, No. 3, Autumn 1985, - pp.223-230. - [2] Snyder, J.P., 1987, Map Projections-A Working Manual. U.S. - Geological Survey Professional Paper 1395. Washington, DC: U.S. - Government Printing Office, pp.15-16 and pp. 44-45. - [3] Thomas, P.D., 1952, Conformal Projections in Geodesy and - Cartography, Special Publication No. 251, Coast and Geodetic - Survey, U.S. Department of Commerce, Washington, DC: U.S. - Government Printing Office, p. 66. - """ + if deg is True: + lat1, lat2 = np.radians([lat1, lat2]) # set ellipsoid parameters if ell is None: ell = Ellipsoid() + + rlat1 = geodetic2rectifying(lat1, ell, deg=False) + rlat2 = geodetic2rectifying(lat2, ell, deg=False) - if deg is True: - lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2]) + return rsphere_rectifying(ell) * abs(rlat2 - rlat1) - # compute isometric latitude of P1 and P2 - isolat1 = isometric(lat1, deg=False, ell=ell) - isolat2 = isometric(lat2, deg=False, ell=ell) +def departure(lon1: float, lon2: float, lat: float, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + Computes is the distance along a specific parallel between two meridians. + ''' - # compute changes in isometric latitude and longitude between points - disolat = isolat2 - isolat1 - dlon = lon2 - lon1 + if ell is None: + ell = Ellipsoid() - # compute azimuth - az12 = np.arctan2(dlon, disolat) + if deg is True: + lon1, lon2, lat = np.radians([lon1, lon2, lat]) + + return rcurve_parallel(lat, ell, deg=False) * abs(wrapToPi(lon2-lon1)) - # compute distance along loxodromic curve - m1 = meridian_dist(lat1, deg=False, ell=ell) - m2 = meridian_dist(lat2, deg=False, ell=ell) - dm = m2 - m1 - lox_s = dm / np.cos(az12) + +def meanm(lat: float, lon: float, ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float]: + '''Computes geographic mean for geographic points on an ellipsoid''' + + if ell is None: + ell = Ellipsoid() if deg is True: - az12 = np.degrees(az12) % 360. + lon, lat = radians([lon, lat]) + + lat = geodetic2authalic(lat) + x,y,z = sph2cart(lon,lat,np.ones_like(lat)) + latbar, lonbar, _ =cart2sph(np.sum(x),np.sum(y),np.sum(z)) + latbar = authalic2geodetic(latbar, ell, deg=False) + lonbar = wrapToPi(lonbar) - return lox_s, az12 + return latbar, lonbar diff --git a/pymap3d/rcurve.py b/pymap3d/rcurve.py new file mode 100644 index 0000000..08b6809 --- /dev/null +++ b/pymap3d/rcurve.py @@ -0,0 +1,47 @@ +"""compute radii of curvature for an ellipsoid""" +import numpy as np +from .ellipsoid import Ellipsoid +from numpy import radians, sin, cos, sqrt + +__all__ = ['rcurve_parallel', 'rcurve_meridian', 'rcurve_transverse'] + +def rcurve_parallel(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + computes the radius of the small circle encompassing the globe at the specified latitude + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + f1 = 1 - ell.eccentricity**2 + f2 = 1 - (ell.eccentricity * cos(lat))**2 + return cos(lat) * ell.semimajor_axis * sqrt(f1 / f2) + +def rcurve_meridian(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + '''computes the meridional radius of curvature for the ellipsoid''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + f1 = ell.semimajor_axis * (1 - ell.eccentricity ** 2) + f2 = 1 - (ell.eccentricity * sin(lat)) ** 2 + return f1 / sqrt(f2 ** 3) + +def rcurve_transverse(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + '''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''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + return ell.semimajor_axis / (1 - (ell.eccentricity * sin(lat)) ** 2) \ No newline at end of file diff --git a/pymap3d/rsphere.py b/pymap3d/rsphere.py new file mode 100644 index 0000000..c6d040a --- /dev/null +++ b/pymap3d/rsphere.py @@ -0,0 +1,109 @@ +""" compute radii of auxiliary spheres""" +from typing import Tuple +import numpy as np +from .ellipsoid import Ellipsoid +from numpy import radians, sin, cos, tan, arctan, log, hypot, degrees, arctan2, sqrt, pi +from .rcurve import rcurve_meridian, rcurve_transverse + +__all__ = ['rsphere_eqavol', 'rsphere_authalic', 'rsphere_rectifying', 'rsphere_euler', 'rsphere_curve', 'rsphere_triaxial', 'rsphere_biaxial'] + +def rsphere_eqavol(ell: Ellipsoid = None) -> float: + '''computes the radius of the sphere with equal volume as the ellipsoid''' + if ell is None: + ell = Ellipsoid() + + return ell.semimajor_axis * (1 - ell.flattening/3 - ell.flattening**2/9) + +def rsphere_authalic(ell: Ellipsoid = None) -> float: + '''computes the radius of the sphere with equal surface area as the ellipsoid''' + if ell is None: + ell = Ellipsoid() + + if ell.eccentricity > 0: + f1 = ell.semimajor_axis**2 / 2 + f2 = (1 - ell.eccentricity**2) / 2*ell.eccentricity + f3 = log((1+ell.eccentricity) / (1-ell.eccentricity)) + return sqrt(f1 * (1 + f2 * f3)) + else: + return ell.semimajor_axis + +def rsphere_rectifying(ell: Ellipsoid = None) -> float: + '''computes the radius of the sphere with equal meridional distances as the ellipsoid''' + if ell is None: + ell = Ellipsoid() + + return ((ell.semimajor_axis**(3/2) + ell.semiminor_axis**(3/2))/2)**(2/3) + + +def rsphere_euler(lat1, lon1, lat2, lon2, ell: Ellipsoid = None, deg: bool = True) -> float: + '''computes the Euler radii of curvature at the midpoint of the + great circle arc defined by the endpoints (lat1,lon1) and (lat2,lon2)''' + + if ell is None: + ell = Ellipsoid() + + if deg is False: + lat1, lon1, lat2, lon2 = degrees(lat1, lon1, lat2, lon2) + + latmid = lat1 + (lat2 - lat1)/2 # compute the midpoint + + # compute azimuth + _,_,az = vdist(lat1,lon1,lat2,lon2,ell) + + # compute meridional and transverse radii of curvature + + rho = rcurve_meridian(latmid, ell, deg=True) + nu = rcurve_transverse(latmid, ell, deg=True) + + az = radians(az) + den = rho * sin(az)**2 + nu * cos(az)**2; + + # compute radius of the arc from point 1 to point 2 + return rho * nu / den + +def rsphere_curve(lat, ell: Ellipsoid = None, deg: bool = True, method='mean') -> float: + '''computes the arithmetic average of the transverse and meridional + radii of curvature at a specified latitude point''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + rho = rcurve_meridian(latmid, ell, deg=False) + nu = rcurve_transverse(latmid, ell, deg=False) + + if method=='mean': + return (rho + nu) / 2 + elif method=='norm': + return sqrt(rho * nu) + else: + raise Exception("pymap3d.rsphere.curve: method must be mean or norm") + +def rsphere_triaxial(ell: Ellipsoid = None, method='mean') -> float: + '''computes triaxial average of the semimajor and semiminor axes of the ellipsoid''' + + if ell is None: + ell = Ellipsoid() + + if method=='mean': + return (2*ell.semimajor_axis + ell.semiminor_axis) / 3 + elif method=='norm': + return (ell.semimajor_axis**2 * ell.semiminor_axis) ** (1/3) + else: + raise Exception("pymap3d.rsphere.rsphere_triaxial: method must be mean or norm") + + +def rsphere_biaxial(ell: Ellipsoid = None, method='mean') -> float: + '''computes biaxial average of the semimajor and semiminor axes of the ellipsoid''' + + if ell is None: + ell = Ellipsoid() + + if method=='mean': + return (ell.semimajor_axis + ell.semiminor_axis) / 2 + elif method=='norm': + return sqrt(ell.semimajor_axis * ell.semiminor_axis) + else: + raise Exception("pymap3d.rsphere.rsphere_biaxial: method must be mean or norm") diff --git a/pymap3d/sidereal.py b/pymap3d/sidereal.py index 95387c5..8df0bcc 100644 --- a/pymap3d/sidereal.py +++ b/pymap3d/sidereal.py @@ -44,9 +44,9 @@ def datetime2sidereal(time: datetime, usevallado = usevallado or Time is None if usevallado: jd = juliandate(str2dt(time)) -# %% Greenwich Sidereal time RADIANS +# Greenwich Sidereal time RADIANS gst = julian2sidereal(jd) -# %% Algorithm 15 p. 188 rotate GST to LOCAL SIDEREAL TIME +# Algorithm 15 p. 188 rotate GST to LOCAL SIDEREAL TIME tsr = gst + lon_radians else: tsr = Time(time).sidereal_time(kind='apparent', @@ -122,7 +122,7 @@ def julian2sidereal(Jdate: float) -> float: tsr = np.empty(jdate.size) for i, jd in enumerate(jdate): - # %% Vallado Eq. 3-42 p. 184, Seidelmann 3.311-1 + # Vallado Eq. 3-42 p. 184, Seidelmann 3.311-1 tUT1 = (jd - 2451545.0) / 36525. # Eqn. 3-47 p. 188 diff --git a/pymap3d/util.py b/pymap3d/util.py new file mode 100644 index 0000000..1ed9e88 --- /dev/null +++ b/pymap3d/util.py @@ -0,0 +1,42 @@ +from numpy import arctan2, hypot, cos, sin, pi +import numpy as np +from typing import Tuple + +__all__ = ['cart2pol', 'pol2cart', 'cart2sph', 'sph2cart', 'wrapToPi', 'wrapTo2Pi'] + +def cart2pol(x: float, y: float) -> Tuple[float, float]: + '''Transform Cartesian to polar coordinates''' + return np.arctan2(y, x), np.hypot(x, y) + +def pol2cart(theta: float, rho: float) -> Tuple[float, float]: + '''Transform polar to Cartesian coordinates''' + return rho * np.cos(theta), rho * np.sin(theta) + +def cart2sph(x: float, y: float, z: float) -> Tuple[float, float, float]: + '''Transform Cartesian to spherical coordinates''' + hxy = np.hypot(x, y) + r = np.hypot(hxy, z) + el = np.arctan2(z, hxy) + az = np.arctan2(y, x) + return az, el, r + +def sph2cart(az: float, el: float, r: float) -> Tuple[float, float, float]: + '''Transform spherical to Cartesian coordinates''' + rcos_theta = r * np.cos(el) + x = rcos_theta * np.cos(az) + y = rcos_theta * np.sin(az) + z = r * np.sin(el) + return x, y, z + +def wrapToPi(x: float) -> float: + '''wraps angles to [-pi, pi)''' + return ((x + pi) % (2 * pi)) - pi + +def wrapTo2Pi(x: float) -> float: + '''wraps angles to [0, 2*pi]''' + return x % (2 * pi) + + + + + diff --git a/pymap3d/vallado.py b/pymap3d/vallado.py index 04e9673..2b92676 100644 --- a/pymap3d/vallado.py +++ b/pymap3d/vallado.py @@ -68,7 +68,7 @@ def azel2radec(az_deg: float, el_deg: float, el = radians(el) lat = radians(lat) lon = radians(lon) -# %% Vallado "algorithm 28" p 268 +# Vallado "algorithm 28" p 268 dec = arcsin(sin(el) * sin(lat) + cos(el) * cos(lat) * cos(az)) lha = arctan2(-(sin(az) * cos(el)) / cos(dec), @@ -135,11 +135,11 @@ def radec2azel(ra_deg: float, dec_deg: float, lon = radians(lon) lst = datetime2sidereal(time, lon) # RADIANS -# %% Eq. 4-11 p. 267 LOCAL HOUR ANGLE +# Eq. 4-11 p. 267 LOCAL HOUR ANGLE lha = lst - ra -# %% #Eq. 4-12 p. 267 +# #Eq. 4-12 p. 267 el = arcsin(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(lha)) -# %% combine Eq. 4-13 and 4-14 p. 268 +# combine Eq. 4-13 and 4-14 p. 268 az = arctan2(-sin(lha) * cos(dec) / cos(el), (sin(dec) - sin(el) * sin(lat)) / (cos(el) * cos(lat))) diff --git a/pymap3d/vincenty.py b/pymap3d/vincenty.py index 5ba72c6..ef247c7 100644 --- a/pymap3d/vincenty.py +++ b/pymap3d/vincenty.py @@ -86,7 +86,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N """ if ell is None: ell = Ellipsoid() -# %% prepare inputs +# prepare inputs lat1 = np.atleast_1d(Lat1) lat2 = np.atleast_1d(Lat2) lon1 = np.atleast_1d(Lon1) @@ -102,21 +102,21 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N if lat2.size == 1: lat2 = np.broadcast_to(lat2, lat1.shape) lon2 = np.broadcast_to(lon2, lon1.shape) -# %% Input check: +# Input check: if ((abs(lat1) > 90) | (abs(lat2) > 90)).any(): raise ValueError('Input latitudes must be in [-90, 90] degrees.') -# %% Supply WGS84 earth ellipsoid axis lengths in meters: +# Supply WGS84 earth ellipsoid axis lengths in meters: a = ell.semimajor_axis b = ell.semiminor_axis -# %% preserve true input latitudes: +# preserve true input latitudes: lat1tr = lat1.copy() lat2tr = lat2.copy() -# %% convert inputs in degrees to np.radians: +# convert inputs in degrees to np.radians: lat1 = np.radians(lat1) lon1 = np.radians(lon1) lat2 = np.radians(lat2) lon2 = np.radians(lon2) -# %% correct for errors at exact poles by adjusting 0.6 millimeters: +# correct for errors at exact poles by adjusting 0.6 millimeters: kidx = abs(pi / 2 - abs(lat1)) < 1e-10 if kidx.any(): lat1[kidx] = sign(lat1[kidx]) * (pi / 2 - (1e-10)) @@ -203,7 +203,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N dist_m = (b * A * (sigma - deltasigma)) -# %% From point #1 to point #2 +# From point #1 to point #2 # correct sign of lambda for azimuth calcs: lamb = abs(lamb) kidx = sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0 @@ -213,12 +213,12 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N a12 = arctan2(numer, denom) kidx = a12 < 0 a12[kidx] = a12[kidx] + 2 * pi - # %% from poles + # from poles a12[lat1tr <= -90] = 0 a12[lat1tr >= 90] = pi az = np.degrees(a12) -# %% From point #2 to point #1 +# From point #2 to point #1 # correct sign of lambda for azimuth calcs: lamb = abs(lamb) kidx = sign(sin(lon1 - lon2)) * sign(sin(lamb)) < 0 @@ -228,7 +228,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N a21 = arctan2(numer, denom) kidx = a21 < 0 a21[kidx] = a21[kidx] + 2 * pi - # %% backwards from poles: + # backwards from poles: a21[lat2tr >= 90] = pi a21[lat2tr <= -90] = 0. backaz = np.degrees(a21) From ea1c3e64507e6af535b64f487a048aa055105970 Mon Sep 17 00:00:00 2001 From: rpavlick Date: Sun, 11 Aug 2019 13:59:34 -0700 Subject: [PATCH 2/5] Adding rcurve.py containing several functions that compute the radii of curvature for an ellipsoid. Adding rsphere.py containing several functions that compute the radii of auxiliary spheres. Adding latitude.py containing several functions that compute auxiliary latitudes and their inverse functions. Adding 'loxodrome_direct' to compute the end point of a loxodromic curve from a start lat/lon, azimuth, and distance. Fixing 'loxodrome_inverse', it previously failed on the degenerate case of east/west azimuth. Adding 'meanm' which computes the geographic mean of points on a spheroid. Adding 'departure' which computes the distance between two longitudes at a specified latitude. Adding 'meridianarc' which computes the distance between two latitudes along a specified longitude. Added thirdflattening and eccentricty attributes to ellipsoid definitions. Removed all instances '%%' in comments. --- pymap3d/__init__.py | 6 +- pymap3d/ecef.py | 2 +- pymap3d/eci.py | 8 +- pymap3d/ellipsoid.py | 22 +- pymap3d/latitude.py | 573 +++++++++++++++++++++++++++++++++++++++++++ pymap3d/los.py | 4 +- pymap3d/lox.py | 289 ++++++++++++---------- pymap3d/rcurve.py | 46 ++++ pymap3d/rsphere.py | 107 ++++++++ pymap3d/sidereal.py | 6 +- pymap3d/util.py | 41 ++++ pymap3d/vallado.py | 8 +- pymap3d/vincenty.py | 20 +- 13 files changed, 982 insertions(+), 150 deletions(-) create mode 100644 pymap3d/latitude.py create mode 100644 pymap3d/rcurve.py create mode 100644 pymap3d/rsphere.py create mode 100644 pymap3d/util.py diff --git a/pymap3d/__init__.py b/pymap3d/__init__.py index 2b733a1..3767edd 100644 --- a/pymap3d/__init__.py +++ b/pymap3d/__init__.py @@ -52,6 +52,10 @@ from .enu import enu2geodetic, geodetic2enu, aer2enu, enu2aer # noqa: F401 from .aer import ecef2aer, aer2ecef, geodetic2aer, aer2geodetic, eci2aer, aer2eci # noqa: F401 from .los import lookAtSpheroid # noqa: F401 - from .lox import isometric, meridian_dist, loxodrome_inverse # noqa: F401 + from .lox import loxodrome_direct, loxodrome_inverse, meridianarc, departure, meanm + from .latitude import geodetic2isometric, isometric2geodetic, geodetic2rectifying, rectifying2geodetic, geodetic2conformal, conformal2geodetic, geodetic2parametric, parametric2geodetic, geodetic2geocentric, geocentric2geodetic, geodetic2authalic, authalic2geodetic + from .rcurve import rcurve_meridian, rcurve_parallel, rcurve_transverse + from .rsphere import rsphere_authalic, rsphere_biaxial, rsphere_curve, rsphere_eqavol, rsphere_euler, rsphere_rectifying, rsphere_triaxial + from .util import wrapToPi, wrapTo2Pi, sph2cart, cart2sph, pol2cart, cart2pol else: # pure Python only from .math import * # type: ignore # noqa: F401, F403 diff --git a/pymap3d/ecef.py b/pymap3d/ecef.py index 7814d7f..f7961b7 100644 --- a/pymap3d/ecef.py +++ b/pymap3d/ecef.py @@ -126,7 +126,7 @@ def ecef2geodetic(x: float, y: float, z: float, eps = ((ell.semiminor_axis * u - ell.semimajor_axis * huE + E**2) * sin(Beta)) / (ell.semimajor_axis * huE * 1 / cos(Beta) - E**2 * cos(Beta)) Beta += eps -# %% final output +# final output lat = arctan(ell.semimajor_axis / ell.semiminor_axis * tan(Beta)) lon = arctan2(y, x) diff --git a/pymap3d/eci.py b/pymap3d/eci.py index 38026fc..f8ae26b 100644 --- a/pymap3d/eci.py +++ b/pymap3d/eci.py @@ -41,14 +41,14 @@ def eci2ecef(x: float, y: float, z: float = None, z : float target z ECEF coordinate """ -# %% +# # FIXME: temporary old API, which was a single N x 3 numpy.ndarray if z is None and isinstance(y, (str, datetime)): time = y z = x[:, 2] y = x[:, 1] x = x[:, 0] -# %% +# x = np.atleast_1d(x) y = np.atleast_1d(y) z = np.atleast_1d(z) @@ -112,14 +112,14 @@ def ecef2eci(x: float, y: float, z: float = None, z : float target z ECI coordinate """ -# %% +# # FIXME: temporary old API, which was a single N x 3 numpy.ndarray if z is None and isinstance(y, (str, datetime)): time = y z = x[:, 2] y = x[:, 1] x = x[:, 0] -# %% +# x = np.atleast_1d(x) y = np.atleast_1d(y) z = np.atleast_1d(z) diff --git a/pymap3d/ellipsoid.py b/pymap3d/ellipsoid.py index 1ed4939..8fe713f 100644 --- a/pymap3d/ellipsoid.py +++ b/pymap3d/ellipsoid.py @@ -1,4 +1,4 @@ - +from math import sqrt class Ellipsoid: """ @@ -26,19 +26,27 @@ def __init__(self, model='wgs84'): self.semimajor_axis = 6378137. self.flattening = 1 / 298.2572235630 self.semiminor_axis = self.semimajor_axis * (1 - self.flattening) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'wgs72': self.semimajor_axis = 6378135. self.flattening = 1 / 298.26 self.semiminor_axis = self.semimajor_axis * (1 - self.flattening) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'grs80': """https://en.wikipedia.org/wiki/GRS_80""" self.semimajor_axis = 6378137. self.flattening = 1 / 298.257222100882711243 self.semiminor_axis = self.semimajor_axis * (1 - self.flattening) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'clarke1866': self.semimajor_axis = 6378206.4 self.semiminor_axis = 6356583.8 self.flattening = -(self.semiminor_axis / self.semimajor_axis - 1) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'mars': """ https://tharsis.gsfc.nasa.gov/geodesy.html @@ -46,18 +54,26 @@ def __init__(self, model='wgs84'): self.semimajor_axis = 3396900. self.semiminor_axis = 3376097.80585952 self.flattening = 1 / 163.295274386012 + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'moon': self.semimajor_axis = 1738000. self.semiminor_axis = self.semimajor_axis self.flattening = 0. + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'venus': self.semimajor_axis = 6051000. self.semiminor_axis = self.semimajor_axis + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) self.flattening = 0. elif model == 'jupiter': self.semimajor_axis = 71492000. self.flattening = 1 / 15.415446277169725 self.flattening = -(self.semiminor_axis / self.semimajor_axis - 1) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'io': """ https://doi.org/10.1006/icar.1998.5987 @@ -65,9 +81,13 @@ def __init__(self, model='wgs84'): self.semimajor_axis = 1829.7 self.flattening = 1 / 131.633093525179 self.semiminor_axis = self.semimajor_axis * (1 - self.flattening) + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) elif model == 'pluto': self.semimajor_axis = 1187000. self.semiminor_axis = self.semimajor_axis self.flattening = 0. + self.thirdflattening = (self.semimajor_axis - self.semiminor_axis) / (self.semimajor_axis + self.semiminor_axis) + self.eccentricity = sqrt(2*self.flattening * self.flattening**2) else: raise NotImplementedError('{} model not implemented, let us know and we will add it (or make a pull request)'.format(model)) diff --git a/pymap3d/latitude.py b/pymap3d/latitude.py new file mode 100644 index 0000000..38f106c --- /dev/null +++ b/pymap3d/latitude.py @@ -0,0 +1,573 @@ +""" compute auxiliary latitudes and their inverse""" +import numpy as np +from .ellipsoid import Ellipsoid +from numpy import radians, sin, tan, exp, arctan, log, degrees, sqrt, pi + +__all__ = ['geodetic2isometric', 'isometric2geodetic', 'geodetic2rectifying', 'rectifying2geodetic', 'geodetic2conformal', 'conformal2geodetic', 'geodetic2parametric', 'parametric2geodetic', 'geodetic2geocentric', 'geocentric2geodetic', 'geodetic2authalic', 'authalic2geodetic'] + +def geodetic2isometric(lat: float, ell: Ellipsoid = None, deg: bool = True) -> float: + """ + computes isometric latitude of a point on an ellipsoid + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + isolat : float or numpy.ndarray of float + isometric latiude + + Notes + ----- + + Isometric latitude is an auxiliary latitude proportional to the spacing + of parallels of latitude on an ellipsoidal mercator projection. + + Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, + School of Mathematical and Geospatial Sciences, RMIT University, + January 2010 + """ + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = np.deg2rad(lat) + + x = ell.eccentricity * np.sin(lat) + y = (1 - x)/(1 + x) + z = np.pi/4 + lat/2 + +# calculate the isometric latitude + isolat = np.log(np.tan(z) * (y**(ell.eccentricity /2))) + + if deg is True: + isolat = np.degrees(isolat) + + return isolat + +def isometric2geodetic(isolat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from isometric latitude to geodetic latitude + + Parameters + ---------- + + isolat : float or numpy.ndarray of float + isometric latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + n = ell.thirdflattening + f1 = 3*n/2 - 27*n**3/32 + f2 = 21*n**2/16 - 55*n**4/32 + f3 = 151*n**3/96 + f4 = 1097*n**4/512 + + if deg is True: + isolat = radians(isolat) + + cnflat = 2 * arctan(exp(isolat)) - (pi/2) + lat = conformal2geodetic(cnflat, ell, deg=False) + + if deg is True: + lat = degrees(lat) + + return lat + +def geodetic2rectifying(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from geodetic latitude to rectifying latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + reclat : float or numpy.ndarray of float + rectifying latiude + + Notes + ----- + + Rectifying latitude is an auxiliary latitude used to map an ellipsoid to a + sphere such that correct distances along meridians are preserved. + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + n = ell.thirdflattening + f1 = 3*n/2 - 9*n**3/16 + f2 = 15*n**2/16 - 15*n**4/32 + f3 = 35*n**3/48 + f4 = 315*n**4/512 + + if deg is True: + lat = radians(lat) + + reclat = lat - f1*sin(2*lat) + f2*sin(4*lat) - f3*sin(6*lat) + f4*sin(8*lat) + + if deg is True: + reclat = degrees(reclat) + + return reclat + + +def rectifying2geodetic(reclat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from rectifying latitude to geodetic latitude + + Parameters + ---------- + + reclat : float or numpy.ndarray of float + rectifying latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + n = ell.thirdflattening + f1 = 3*n/2 - 27*n**3/32 + f2 = 21*n**2/16 - 55*n**4/32 + f3 = 151*n**3/96 + f4 = 1097*n**4/512 + + if deg is True: + reclat = radians(reclat) + + lat = reclat + f1*sin(2*reclat) + f2*sin(4*reclat) + f3*sin(6*reclat) + f4*sin(8*reclat) + + if deg is True: + lat = degrees(lat) + + return lat + + +def geodetic2conformal(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from geodetic latitude to conformal latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + cnflat : float or numpy.ndarray of float + conformal latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + e = ell.eccentricity + f1 = 1 - e*sin(lat) + f2 = 1 + e*sin(lat) + f3 = 1 - sin(lat) + f4 = 1 + sin(lat) + + # compute conformal latitudes with correction for points at +90 + cnflat = np.where(f3 == 0, pi/2, 2 * arctan(sqrt((f4/f3) * ((f1/f2)**e))) - (pi/2)) + + if deg is True: + cnflat = degrees(cnflat) + + return cnflat + +def conformal2geodetic(cnflat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from conformal latitude to geodetic latitude + + Parameters + ---------- + + cnflat : float or numpy.ndarray of float + conformal latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + cnflat = radians(cnflat) + + e = ell.eccentricity + f1 = e**2/2 + 5*e**4/24 + e**6/12 + 13*e**8/360 + f2 = 7*e**4/48 + 29*e**6/240 + 811*e**8/11520 + f3 = 7*e**6/120 + 81*e**8/1120 + f4 = 4279*e**8/161280 + + lat = cnflat + f1*sin(2*cnflat) + f2*sin(4*cnflat) + f3*sin(6*cnflat) + f4*sin(8*cnflat) + + if deg is True: + lat = degrees(lat) + + return lat + +def geodetic2parametric(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from geodetic latitude to parametric latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + parlat : float or numpy.ndarray of float + parametric latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + epsilon = 1E-10 + + lat = np.where(lat == pi/2, pi/2 - epsilon, lat) + lat = np.where(lat == -pi/2, -pi/2 + epsilon, lat) + + parlat = arctan(sqrt(1-(ell.eccentricity)**2) * tan(lat)) + + if deg is True: + parlat = degrees(parlat) + + return parlat + +def parametric2geodetic(parlat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from parametric latitude to geodetic latitude + + Parameters + ---------- + + parlat : float or numpy.ndarray of float + parametric latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + parlat = radians(parlat) + + epsilon = 1E-10 + + parlat = np.where(parlat == pi/2, pi/2 - epsilon, parlat) + parlat = np.where(parlat == -pi/2, -pi/2 + epsilon, parlat) + + lat = arctan(tan(parlat)/ sqrt(1-(ell.eccentricity)**2)) + + if deg is True: + lat = degrees(lat) + + return lat + +def geodetic2geocentric(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from geodetic latitude to geocentric latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + cenlat : float or numpy.ndarray of float + geocentric latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + epsilon = 1E-10 + + lat = np.where(lat == pi/2, pi/2 - epsilon, lat) + lat = np.where(lat == -pi/2, -pi/2 + epsilon, lat) + + cenlat = arctan((1-(ell.eccentricity)**2) * tan(lat)) + + if deg is True: + cenlat = degrees(cenlat) + + return cenlat + +def geocentric2geodetic(cenlat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from geocentric latitude to geodetic latitude + + Parameters + ---------- + + cenlat : float or numpy.ndarray of float + geocentric latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + cenlat = radians(cenlat) + + epsilon = 1E-10 + + cenlat = np.where(cenlat == pi/2, pi/2 - epsilon, cenlat) + cenlat = np.where(cenlat == -pi/2, -pi/2 + epsilon, cenlat) + + lat = arctan(tan(cenlat)/ (1-(ell.eccentricity)**2)) + + if deg is True: + lat = degrees(lat) + + return lat + +def geodetic2authalic(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from geodetic latitude to authalic latitude + + Parameters + ---------- + + lat : float or numpy.ndarray of float + geodetic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + autlat : float or numpy.ndarray of float + authalic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + e = ell.eccentricity + f1 = e**2/3 + 31*e**4/180 + 59*e**6/560 + f2 = 17*e**4/360 + 61*e**6/1260 + f3 = 383*e**6/45360 + + autlat = lat - f1*sin(2*lat) + f2*sin(4*lat) - f3*sin(6*lat) + + if deg is True: + autlat = degrees(autlat) + + return autlat + +def authalic2geodetic(autlat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + converts from authalic latitude to geodetic latitude + + Parameters + ---------- + + autlat : float or numpy.ndarray of float + authalic latitude + ell : Ellipsoid, optional + reference ellipsoid (default WGS84) + deg : bool, optional + degrees input/output (False: radians in/out) + + Returns + ------- + + lat : float or numpy.ndarray of float + geodetic latiude + + Notes + ----- + + Equations from J. P. Snyder, "Map Projections - A Working Manual", + US Geological Survey Professional Paper 1395, US Government Printing + Office, Washington, DC, 1987, pp. 13-18. +''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + autlat = radians(autlat) + + e = ell.eccentricity + f1 = e**2/3 + 31*e**4/180 + 517*e**6/5040 + f2 = 23*e**4/360 + 251*e**6/3780 + f3 = 761*e**6/45360 + + lat = autlat + f1*sin(2*autlat) + f2*sin(4*autlat) + f3*sin(6*autlat) + + if deg is True: + lat = degrees(lat) + + return lat \ No newline at end of file diff --git a/pymap3d/los.py b/pymap3d/los.py index ef175fb..bfe47a0 100644 --- a/pymap3d/los.py +++ b/pymap3d/los.py @@ -72,13 +72,13 @@ def lookAtSpheroid(lat0: float, lon0: float, h0: float, az: float, tilt: float, 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 +# Return nan if radical < 0 or d < 0 because LOS vector does not point towards Earth with np.errstate(invalid='ignore'): d = np.where(radical > 0, (value - a * b * c * np.sqrt(radical)) / magnitude, np.nan) d[d < 0] = np.nan -# %% cartesian to ellipsodal +# cartesian to ellipsodal lat, lon, _ = ecef2geodetic(x + d * u, y + d * v, z + d * w, deg=deg) return lat, lon, d diff --git a/pymap3d/lox.py b/pymap3d/lox.py index 05671af..08a7f04 100644 --- a/pymap3d/lox.py +++ b/pymap3d/lox.py @@ -1,149 +1,192 @@ -""" isometric latitude, meridian distance """ import numpy as np from .ellipsoid import Ellipsoid +from typing import Tuple +from .latitude import geodetic2isometric, geodetic2rectifying, rectifying2geodetic, geodetic2authalic, authalic2geodetic +from .rsphere import rsphere_rectifying +from .rcurve import rcurve_parallel +from .util import wrapToPi, sph2cart, cart2sph +from numpy import tan, sign, sin, cos, pi, radians, degrees +__all__ = ['loxodrome_inverse', 'loxodrome_direct', 'meridianarc', 'departure', 'meanm'] -def isometric(lat: float, ell: Ellipsoid = None, deg: bool = True): +def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float, + ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float]: """ - computes isometric latitude of a point on an ellipsoid + computes the arc length and azimuth of the loxodrome + between two points on the surface of the reference ellipsoid - Parameters - ---------- + Parameters + ---------- - lat : float or numpy.ndarray of float - geodetic latitude - ell : Ellipsoid, optional + lat1 : float or numpy.ndarray of float + geodetic latitude of first point + lon1 : float or numpy.ndarray of float + geodetic longitude of first point + lat2 : float or numpy.ndarray of float + geodetic latitude of second point + lon2 : float or numpy.ndarray of float + geodetic longitude of second point + ell : Ellipsoid, optional reference ellipsoid (default WGS84) - deg : bool, optional + deg : bool, optional degrees input/output (False: radians in/out) - Returns - ------- - - isolat : float or numpy.ndarray of float - isometric latiude + Results + ------- - Notes - ----- + dist : float or numpy.ndarray of float + distance along loxodrome + a12 : float or numpy.ndarray of float + azimuth of loxodrome (degrees/radians) - Isometric latitude is an auxiliary latitude proportional to the spacing - of parallels of latitude on an ellipsoidal mercator projection. + Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, + School of Mathematical and Geospatial Sciences, RMIT University, January 2010 - Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, - School of Mathematical and Geospatial Sciences, RMIT University, - January 2010 + [1] Bowring, B.R., 1985, 'The geometry of the loxodrome on the + ellipsoid', The Canadian Surveyor, Vol. 39, No. 3, Autumn 1985, + pp.223-230. + [2] Snyder, J.P., 1987, Map Projections-A Working Manual. U.S. + Geological Survey Professional Paper 1395. Washington, DC: U.S. + Government Printing Office, pp.15-16 and pp. 44-45. + [3] Thomas, P.D., 1952, Conformal Projections in Geodesy and + Cartography, Special Publication No. 251, Coast and Geodetic + Survey, U.S. Department of Commerce, Washington, DC: U.S. + Government Printing Office, p. 66. """ + # set ellipsoid parameters if ell is None: ell = Ellipsoid() - f = ell.flattening # flattening of ellipsoid - if deg is True: - lat = np.deg2rad(lat) + lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2]) - e2 = f * (2 - f) # eccentricity-squared - e = np.sqrt(e2) # eccentricity of ellipsoid + # compute isometric latitude of P1 and P2 + isolat1 = geodetic2isometric(lat1, deg=False, ell=ell) + isolat2 = geodetic2isometric(lat2, deg=False, ell=ell) - x = e * np.sin(lat) - y = (1 - x) / (1 + x) - z = np.pi / 4 + lat / 2 + # compute changes in isometric latitude and longitude between points + disolat = isolat2 - isolat1 -# calculate the isometric latitude - isolat = np.log(np.tan(z) * (y**(e / 2))) + dlon = abs(wrapToPi(lon2-lon1)) + + # compute azimuth + a12 = np.arctan2(dlon, disolat) + cosaz = abs(cos(a12)) + + # compute distance along loxodromic curve + dist = meridianarc(lat2, lat1, deg=False, ell=ell) / cosaz + + # consider degenerate case (directly east/west) + epsilon = 1e-10 + dist = np.where(cosaz < epsilon,departure(lon1, lon2, (lat1+lat2)/2, ell, deg=False), dist) if deg is True: - isolat = np.degrees(isolat) + a12 = np.degrees(a12) % 360. - return isolat + return dist, a12 -def meridian_dist(lat: float, ell: Ellipsoid = None, deg: bool = True): +def loxodrome_direct(lat1: float, lon1: float, rng: float, a12: float, + ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float, float]: """ - computes the meridian distance on an ellipsoid *from the equator* to a latitude. Parameters ---------- - lat : float or numpy.ndarray of float - geodetic latitude + lat1 : float + inital geodetic latitude (degrees) + lon1 : float + initial geodetic longitude (degrees) + rng : float + ground distance (meters) + a12 : float + azimuth (degrees) clockwide from north. ell : Ellipsoid, optional - reference ellipsoid (default WGS84) - deg : bool, optional - degrees input/output (False: radians in/out) + reference ellipsoid Results ------- - mdist : float or numpy.ndarray of float - meridian distance (degrees/radians) + lat2 : float + final geodetic latitude (degrees) + lon2 : float + final geodetic longitude (degrees) + a21 : float + reverse azimuth (degrees), at final point facing back toward the intial point +""" + if ell is None: + ell = Ellipsoid() - Notes - ----- + lat1 = np.atleast_1d(lat1) + lon1 = np.atleast_1d(lon1) + rng = np.atleast_1d(rng) + a12 = np.atleast_1d(a12) - Formula given Baeschlin, C.F., 1948, - "Lehrbuch Der Geodasie", Orell Fussli Verlag, Zurich, pp.47-50. + if rng.ndim != 1 or a12.ndim != 1: + raise ValueError('Range and azimuth must be scalar or vector') - Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, - School of Mathematical and Geospatial Sciences, RMIT University, January 2010 - """ + if lat1.size > 1 and rng.size > 1: + raise ValueError('LOXODROME_DIRECT: Variable ranges are only allowed for a single point.') + + if rng.size != a12.size and rng.size == 1: + rng = np.broadcast_to(rng, a12.size) if deg is True: - lat = np.radians(lat) + lat1, lon1, a12 = radians([lat1, lon1, a12]) + + if abs(lat1) > (pi/2): + raise ValueError('LOXODROME_DIRECT: Input lat. must be between -90 and 90 deg., inclusive.') - # set ellipsoid parameters - if ell is None: - ell = Ellipsoid() + # if range is negative, use inverse bearing + a12 = np.where(sign(rng) == -1, (a12 + pi) % (2*pi), a12) + rng = abs(rng) - a = ell.semimajor_axis - f = ell.flattening # flattening of ellipsoid + # compute rectifying sphere latitude and radius + reclat = geodetic2rectifying(lat1, ell, deg=False) - e2 = f * (2 - f) # eccentricity-squared + lat2 = np.zeros_like(lat1) + lon2 = np.zeros_like(lon1) + epsilon = 1e-10 # Set tolerance (should set this for specific machine) - # powers of eccentricity - e4 = e2 * e2 - e6 = e4 * e2 - e8 = e6 * e2 - e10 = e8 * e2 + # correct azimuths at either pole. + a12 = np.where(lat1 >= (pi/2) - epsilon, pi, a12) + a12 = np.where(lat1 <= (-pi/2) + epsilon, 0, a12) - # coefficients of series expansion for meridian distance - A = 1 + (3 / 4) * e2 + (45 / 64) * e4 + (175 / 256) * e6 + (11025 / 16384) * e8 + (43659 / 65536) * e10 - B = (3 / 4) * e2 + (15 / 16) * e4 + (525 / 512) * e6 + (2205 / 2048) * e8 + (72765 / 65536) * e10 - C = (15 / 64) * e4 + (105 / 256) * e6 + (2205 / 4096) * e8 + (10395 / 16384) * e10 - D = (35 / 512) * e6 + (315 / 2048) * e8 + (31185 / 131072) * e10 - E = (315 / 16384) * e8 + (3465 / 65536) * e10 - F = (693 / 131072) * e10 + # compute the new points + cosaz = cos(a12) + lat2 = reclat + (rng/rsphere_rectifying(ell))*cosaz # compute rectifying latitude + lat2 = rectifying2geodetic(lat2,ell,deg=False) # transform to geodetic latitude - term1 = A * lat - term2 = (B / 2) * np.sin(2 * lat) - term3 = (C / 4) * np.sin(4 * lat) - term4 = (D / 6) * np.sin(6 * lat) - term5 = (E / 8) * np.sin(8 * lat) - term6 = (F / 10) * np.sin(10 * lat) + # nudge latitudes near either pole + lat2 = np.where(lat2 >= (pi/2) - epsilon, (pi/2) - epsilon, lat2) + lat2 = np.where(lat2 <= (-pi/2) + epsilon, (-pi/2) + epsilon , lat2) + lat2 = np.where(abs(cosaz)<=epsilon, lat1, lat2) - mdist = a * (1 - e2) * (term1 - term2 + term3 - term4 + term5 - term6) + newiso = geodetic2isometric(lat2,ell,deg=False) + iso = geodetic2isometric(lat1,ell,deg=False) + + dlon = np.where(abs(cosaz)<=epsilon, sign(sin(a12)) * rng/rcurve_parallel(lat1, ell, deg=False), tan(a12) * (newiso - iso)) + lon2 = lon1 + dlon - return mdist + lon2 = wrapToPi(lon2) + a21 = (a12 + pi)%(2*pi) + if deg is True: + lat2, lon2, a21 = degrees([lat2, lon2, a21]) -def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float, - ell: Ellipsoid = None, deg: bool = True): - """ - computes the arc length and azimuth of the loxodrome - between two points on the surface of the reference ellipsoid + return lat2, lon2, a21 + +def meridianarc(lat1: float, lat2: float, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + Computes the meridian distance on an ellipsoid between two latitudes. Parameters ---------- - lat1 : float or numpy.ndarray of float - geodetic latitude of first point - lon1 : float or numpy.ndarray of float - geodetic longitude of first point - lat2 : float or numpy.ndarray of float - geodetic latitude of second point - lon2 : float or numpy.ndarray of float - geodetic longitude of second point + lat1, lat2 : float or numpy.ndarray of float + geodetic latitudes ell : Ellipsoid, optional reference ellipsoid (default WGS84) deg : bool, optional @@ -152,51 +195,49 @@ def loxodrome_inverse(lat1: float, lon1: float, lat2: float, lon2: float, Results ------- - lox_s : float or numpy.ndarray of float - distance along loxodrome - az12 : float or numpy.ndarray of float - azimuth of loxodrome (degrees/radians) - - Based on Deakin, R.E., 2010, 'The Loxodrome on an Ellipsoid', Lecture Notes, - School of Mathematical and Geospatial Sciences, RMIT University, January 2010 + dist : float or numpy.ndarray of float + distance (units same as ellipsoid) + ''' - [1] Bowring, B.R., 1985, 'The geometry of the loxodrome on the - ellipsoid', The Canadian Surveyor, Vol. 39, No. 3, Autumn 1985, - pp.223-230. - [2] Snyder, J.P., 1987, Map Projections-A Working Manual. U.S. - Geological Survey Professional Paper 1395. Washington, DC: U.S. - Government Printing Office, pp.15-16 and pp. 44-45. - [3] Thomas, P.D., 1952, Conformal Projections in Geodesy and - Cartography, Special Publication No. 251, Coast and Geodetic - Survey, U.S. Department of Commerce, Washington, DC: U.S. - Government Printing Office, p. 66. - """ + if deg is True: + lat1, lat2 = np.radians([lat1, lat2]) # set ellipsoid parameters if ell is None: ell = Ellipsoid() + + rlat1 = geodetic2rectifying(lat1, ell, deg=False) + rlat2 = geodetic2rectifying(lat2, ell, deg=False) - if deg is True: - lat1, lon1, lat2, lon2 = np.radians([lat1, lon1, lat2, lon2]) + return rsphere_rectifying(ell) * abs(rlat2 - rlat1) - # compute isometric latitude of P1 and P2 - isolat1 = isometric(lat1, deg=False, ell=ell) - isolat2 = isometric(lat2, deg=False, ell=ell) +def departure(lon1: float, lon2: float, lat: float, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + Computes is the distance along a specific parallel between two meridians. + ''' - # compute changes in isometric latitude and longitude between points - disolat = isolat2 - isolat1 - dlon = lon2 - lon1 + if ell is None: + ell = Ellipsoid() - # compute azimuth - az12 = np.arctan2(dlon, disolat) + if deg is True: + lon1, lon2, lat = np.radians([lon1, lon2, lat]) + + return rcurve_parallel(lat, ell, deg=False) * abs(wrapToPi(lon2-lon1)) - # compute distance along loxodromic curve - m1 = meridian_dist(lat1, deg=False, ell=ell) - m2 = meridian_dist(lat2, deg=False, ell=ell) - dm = m2 - m1 - lox_s = dm / np.cos(az12) + +def meanm(lat: float, lon: float, ell: Ellipsoid = None, deg: bool = True) -> Tuple[float, float]: + '''Computes geographic mean for geographic points on an ellipsoid''' + + if ell is None: + ell = Ellipsoid() if deg is True: - az12 = np.degrees(az12) % 360. + lon, lat = radians([lon, lat]) + + lat = geodetic2authalic(lat) + x,y,z = sph2cart(lon,lat,np.ones_like(lat)) + latbar, lonbar, _ =cart2sph(np.sum(x),np.sum(y),np.sum(z)) + latbar = authalic2geodetic(latbar, ell, deg=False) + lonbar = wrapToPi(lonbar) - return lox_s, az12 + return latbar, lonbar diff --git a/pymap3d/rcurve.py b/pymap3d/rcurve.py new file mode 100644 index 0000000..90de10e --- /dev/null +++ b/pymap3d/rcurve.py @@ -0,0 +1,46 @@ +"""compute radii of curvature for an ellipsoid""" +from .ellipsoid import Ellipsoid +from numpy import radians, sin, cos, sqrt + +__all__ = ['rcurve_parallel', 'rcurve_meridian', 'rcurve_transverse'] + +def rcurve_parallel(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + ''' + computes the radius of the small circle encompassing the globe at the specified latitude + ''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + f1 = 1 - ell.eccentricity**2 + f2 = 1 - (ell.eccentricity * cos(lat))**2 + return cos(lat) * ell.semimajor_axis * sqrt(f1 / f2) + +def rcurve_meridian(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + '''computes the meridional radius of curvature for the ellipsoid''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + f1 = ell.semimajor_axis * (1 - ell.eccentricity ** 2) + f2 = 1 - (ell.eccentricity * sin(lat)) ** 2 + return f1 / sqrt(f2 ** 3) + +def rcurve_transverse(lat, ell: Ellipsoid = None, deg: bool = True) -> float: + '''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''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + return ell.semimajor_axis / (1 - (ell.eccentricity * sin(lat)) ** 2) \ No newline at end of file diff --git a/pymap3d/rsphere.py b/pymap3d/rsphere.py new file mode 100644 index 0000000..1321f68 --- /dev/null +++ b/pymap3d/rsphere.py @@ -0,0 +1,107 @@ +""" compute radii of auxiliary spheres""" +from .ellipsoid import Ellipsoid +from numpy import radians, sin, cos, tan, log, degrees, sqrt +from .rcurve import rcurve_meridian, rcurve_transverse + +__all__ = ['rsphere_eqavol', 'rsphere_authalic', 'rsphere_rectifying', 'rsphere_euler', 'rsphere_curve', 'rsphere_triaxial', 'rsphere_biaxial'] + +def rsphere_eqavol(ell: Ellipsoid = None) -> float: + '''computes the radius of the sphere with equal volume as the ellipsoid''' + if ell is None: + ell = Ellipsoid() + + return ell.semimajor_axis * (1 - ell.flattening/3 - ell.flattening**2/9) + +def rsphere_authalic(ell: Ellipsoid = None) -> float: + '''computes the radius of the sphere with equal surface area as the ellipsoid''' + if ell is None: + ell = Ellipsoid() + + if ell.eccentricity > 0: + f1 = ell.semimajor_axis**2 / 2 + f2 = (1 - ell.eccentricity**2) / 2*ell.eccentricity + f3 = log((1+ell.eccentricity) / (1-ell.eccentricity)) + return sqrt(f1 * (1 + f2 * f3)) + else: + return ell.semimajor_axis + +def rsphere_rectifying(ell: Ellipsoid = None) -> float: + '''computes the radius of the sphere with equal meridional distances as the ellipsoid''' + if ell is None: + ell = Ellipsoid() + + return ((ell.semimajor_axis**(3/2) + ell.semiminor_axis**(3/2))/2)**(2/3) + + +def rsphere_euler(lat1, lon1, lat2, lon2, ell: Ellipsoid = None, deg: bool = True) -> float: + '''computes the Euler radii of curvature at the midpoint of the + great circle arc defined by the endpoints (lat1,lon1) and (lat2,lon2)''' + + if ell is None: + ell = Ellipsoid() + + if deg is False: + lat1, lon1, lat2, lon2 = degrees(lat1, lon1, lat2, lon2) + + latmid = lat1 + (lat2 - lat1)/2 # compute the midpoint + + # compute azimuth + _,_,az = vdist(lat1,lon1,lat2,lon2,ell) + + # compute meridional and transverse radii of curvature + + rho = rcurve_meridian(latmid, ell, deg=True) + nu = rcurve_transverse(latmid, ell, deg=True) + + az = radians(az) + den = rho * sin(az)**2 + nu * cos(az)**2; + + # compute radius of the arc from point 1 to point 2 + return rho * nu / den + +def rsphere_curve(lat, ell: Ellipsoid = None, deg: bool = True, method='mean') -> float: + '''computes the arithmetic average of the transverse and meridional + radii of curvature at a specified latitude point''' + + if ell is None: + ell = Ellipsoid() + + if deg is True: + lat = radians(lat) + + rho = rcurve_meridian(lat, ell, deg=False) + nu = rcurve_transverse(lat, ell, deg=False) + + if method=='mean': + return (rho + nu) / 2 + elif method=='norm': + return sqrt(rho * nu) + else: + raise Exception("pymap3d.rsphere.curve: method must be mean or norm") + +def rsphere_triaxial(ell: Ellipsoid = None, method='mean') -> float: + '''computes triaxial average of the semimajor and semiminor axes of the ellipsoid''' + + if ell is None: + ell = Ellipsoid() + + if method=='mean': + return (2*ell.semimajor_axis + ell.semiminor_axis) / 3 + elif method=='norm': + return (ell.semimajor_axis**2 * ell.semiminor_axis) ** (1/3) + else: + raise Exception("pymap3d.rsphere.rsphere_triaxial: method must be mean or norm") + + +def rsphere_biaxial(ell: Ellipsoid = None, method='mean') -> float: + '''computes biaxial average of the semimajor and semiminor axes of the ellipsoid''' + + if ell is None: + ell = Ellipsoid() + + if method=='mean': + return (ell.semimajor_axis + ell.semiminor_axis) / 2 + elif method=='norm': + return sqrt(ell.semimajor_axis * ell.semiminor_axis) + else: + raise Exception("pymap3d.rsphere.rsphere_biaxial: method must be mean or norm") diff --git a/pymap3d/sidereal.py b/pymap3d/sidereal.py index 95387c5..8df0bcc 100644 --- a/pymap3d/sidereal.py +++ b/pymap3d/sidereal.py @@ -44,9 +44,9 @@ def datetime2sidereal(time: datetime, usevallado = usevallado or Time is None if usevallado: jd = juliandate(str2dt(time)) -# %% Greenwich Sidereal time RADIANS +# Greenwich Sidereal time RADIANS gst = julian2sidereal(jd) -# %% Algorithm 15 p. 188 rotate GST to LOCAL SIDEREAL TIME +# Algorithm 15 p. 188 rotate GST to LOCAL SIDEREAL TIME tsr = gst + lon_radians else: tsr = Time(time).sidereal_time(kind='apparent', @@ -122,7 +122,7 @@ def julian2sidereal(Jdate: float) -> float: tsr = np.empty(jdate.size) for i, jd in enumerate(jdate): - # %% Vallado Eq. 3-42 p. 184, Seidelmann 3.311-1 + # Vallado Eq. 3-42 p. 184, Seidelmann 3.311-1 tUT1 = (jd - 2451545.0) / 36525. # Eqn. 3-47 p. 188 diff --git a/pymap3d/util.py b/pymap3d/util.py new file mode 100644 index 0000000..48a856c --- /dev/null +++ b/pymap3d/util.py @@ -0,0 +1,41 @@ +from numpy import arctan2, hypot, cos, sin, pi +from typing import Tuple + +__all__ = ['cart2pol', 'pol2cart', 'cart2sph', 'sph2cart', 'wrapToPi', 'wrapTo2Pi'] + +def cart2pol(x: float, y: float) -> Tuple[float, float]: + '''Transform Cartesian to polar coordinates''' + return arctan2(y, x), hypot(x, y) + +def pol2cart(theta: float, rho: float) -> Tuple[float, float]: + '''Transform polar to Cartesian coordinates''' + return rho * cos(theta), rho * sin(theta) + +def cart2sph(x: float, y: float, z: float) -> Tuple[float, float, float]: + '''Transform Cartesian to spherical coordinates''' + hxy = hypot(x, y) + r = hypot(hxy, z) + el = arctan2(z, hxy) + az = arctan2(y, x) + return az, el, r + +def sph2cart(az: float, el: float, r: float) -> Tuple[float, float, float]: + '''Transform spherical to Cartesian coordinates''' + rcos_theta = r * cos(el) + x = rcos_theta * cos(az) + y = rcos_theta * sin(az) + z = r * sin(el) + return x, y, z + +def wrapToPi(x: float) -> float: + '''wraps angles to [-pi, pi)''' + return ((x + pi) % (2 * pi)) - pi + +def wrapTo2Pi(x: float) -> float: + '''wraps angles to [0, 2*pi]''' + return x % (2 * pi) + + + + + diff --git a/pymap3d/vallado.py b/pymap3d/vallado.py index 04e9673..2b92676 100644 --- a/pymap3d/vallado.py +++ b/pymap3d/vallado.py @@ -68,7 +68,7 @@ def azel2radec(az_deg: float, el_deg: float, el = radians(el) lat = radians(lat) lon = radians(lon) -# %% Vallado "algorithm 28" p 268 +# Vallado "algorithm 28" p 268 dec = arcsin(sin(el) * sin(lat) + cos(el) * cos(lat) * cos(az)) lha = arctan2(-(sin(az) * cos(el)) / cos(dec), @@ -135,11 +135,11 @@ def radec2azel(ra_deg: float, dec_deg: float, lon = radians(lon) lst = datetime2sidereal(time, lon) # RADIANS -# %% Eq. 4-11 p. 267 LOCAL HOUR ANGLE +# Eq. 4-11 p. 267 LOCAL HOUR ANGLE lha = lst - ra -# %% #Eq. 4-12 p. 267 +# #Eq. 4-12 p. 267 el = arcsin(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(lha)) -# %% combine Eq. 4-13 and 4-14 p. 268 +# combine Eq. 4-13 and 4-14 p. 268 az = arctan2(-sin(lha) * cos(dec) / cos(el), (sin(dec) - sin(el) * sin(lat)) / (cos(el) * cos(lat))) diff --git a/pymap3d/vincenty.py b/pymap3d/vincenty.py index 5ba72c6..ef247c7 100644 --- a/pymap3d/vincenty.py +++ b/pymap3d/vincenty.py @@ -86,7 +86,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N """ if ell is None: ell = Ellipsoid() -# %% prepare inputs +# prepare inputs lat1 = np.atleast_1d(Lat1) lat2 = np.atleast_1d(Lat2) lon1 = np.atleast_1d(Lon1) @@ -102,21 +102,21 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N if lat2.size == 1: lat2 = np.broadcast_to(lat2, lat1.shape) lon2 = np.broadcast_to(lon2, lon1.shape) -# %% Input check: +# Input check: if ((abs(lat1) > 90) | (abs(lat2) > 90)).any(): raise ValueError('Input latitudes must be in [-90, 90] degrees.') -# %% Supply WGS84 earth ellipsoid axis lengths in meters: +# Supply WGS84 earth ellipsoid axis lengths in meters: a = ell.semimajor_axis b = ell.semiminor_axis -# %% preserve true input latitudes: +# preserve true input latitudes: lat1tr = lat1.copy() lat2tr = lat2.copy() -# %% convert inputs in degrees to np.radians: +# convert inputs in degrees to np.radians: lat1 = np.radians(lat1) lon1 = np.radians(lon1) lat2 = np.radians(lat2) lon2 = np.radians(lon2) -# %% correct for errors at exact poles by adjusting 0.6 millimeters: +# correct for errors at exact poles by adjusting 0.6 millimeters: kidx = abs(pi / 2 - abs(lat1)) < 1e-10 if kidx.any(): lat1[kidx] = sign(lat1[kidx]) * (pi / 2 - (1e-10)) @@ -203,7 +203,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N dist_m = (b * A * (sigma - deltasigma)) -# %% From point #1 to point #2 +# From point #1 to point #2 # correct sign of lambda for azimuth calcs: lamb = abs(lamb) kidx = sign(sin(lon2 - lon1)) * sign(sin(lamb)) < 0 @@ -213,12 +213,12 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N a12 = arctan2(numer, denom) kidx = a12 < 0 a12[kidx] = a12[kidx] + 2 * pi - # %% from poles + # from poles a12[lat1tr <= -90] = 0 a12[lat1tr >= 90] = pi az = np.degrees(a12) -# %% From point #2 to point #1 +# From point #2 to point #1 # correct sign of lambda for azimuth calcs: lamb = abs(lamb) kidx = sign(sin(lon1 - lon2)) * sign(sin(lamb)) < 0 @@ -228,7 +228,7 @@ def vdist(Lat1: float, Lon1: float, Lat2: float, Lon2: float, ell: Ellipsoid = N a21 = arctan2(numer, denom) kidx = a21 < 0 a21[kidx] = a21[kidx] + 2 * pi - # %% backwards from poles: + # backwards from poles: a21[lat2tr >= 90] = pi a21[lat2tr <= -90] = 0. backaz = np.degrees(a21) From 65ef88f46884227046cca851d14fa97f1bd8506d Mon Sep 17 00:00:00 2001 From: Ryan Pavlick Date: Sun, 11 Aug 2019 19:04:22 -0700 Subject: [PATCH 3/5] Update util.py --- pymap3d/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymap3d/util.py b/pymap3d/util.py index 9cb91cb..48a856c 100644 --- a/pymap3d/util.py +++ b/pymap3d/util.py @@ -25,7 +25,7 @@ def sph2cart(az: float, el: float, r: float) -> Tuple[float, float, float]: x = rcos_theta * cos(az) y = rcos_theta * sin(az) z = r * sin(el) - return x, y, z + return x, y, z def wrapToPi(x: float) -> float: '''wraps angles to [-pi, pi)''' From 100bb7fcdd69aa6ce6e652b06a2bc9260a9fb066 Mon Sep 17 00:00:00 2001 From: Ryan Pavlick Date: Sun, 11 Aug 2019 19:37:45 -0700 Subject: [PATCH 4/5] Update latitude.py --- pymap3d/latitude.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/pymap3d/latitude.py b/pymap3d/latitude.py index 38f106c..a609367 100644 --- a/pymap3d/latitude.py +++ b/pymap3d/latitude.py @@ -85,12 +85,6 @@ def isometric2geodetic(isolat, ell: Ellipsoid = None, deg: bool = True) -> float if ell is None: ell = Ellipsoid() - n = ell.thirdflattening - f1 = 3*n/2 - 27*n**3/32 - f2 = 21*n**2/16 - 55*n**4/32 - f3 = 151*n**3/96 - f4 = 1097*n**4/512 - if deg is True: isolat = radians(isolat) @@ -570,4 +564,4 @@ def authalic2geodetic(autlat, ell: Ellipsoid = None, deg: bool = True) -> float: if deg is True: lat = degrees(lat) - return lat \ No newline at end of file + return lat From 982b946e501443c3e18a3a745ef226cfff0c1cf1 Mon Sep 17 00:00:00 2001 From: Ryan Pavlick Date: Sun, 11 Aug 2019 19:38:46 -0700 Subject: [PATCH 5/5] Update rsphere.py --- pymap3d/rsphere.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymap3d/rsphere.py b/pymap3d/rsphere.py index 1321f68..ac59f7d 100644 --- a/pymap3d/rsphere.py +++ b/pymap3d/rsphere.py @@ -1,6 +1,6 @@ """ compute radii of auxiliary spheres""" from .ellipsoid import Ellipsoid -from numpy import radians, sin, cos, tan, log, degrees, sqrt +from numpy import radians, sin, cos, log, degrees, sqrt from .rcurve import rcurve_meridian, rcurve_transverse __all__ = ['rsphere_eqavol', 'rsphere_authalic', 'rsphere_rectifying', 'rsphere_euler', 'rsphere_curve', 'rsphere_triaxial', 'rsphere_biaxial']