Coverage for pygeodesy/etm.py : 81%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
I{Charles Karney's} C++ class U{TransverseMercatorExact <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1TransverseMercatorExact.html>}, abbreviated as C{TMExact} below.
Python class L{ExactTransverseMercator} implements the C{Exact Transverse Mercator} (ETM) projection. Instances of class L{Etm} represent ETM C{easting, nothing} locations.
Following is a copy of Karney's U{TransverseMercatorExact.hpp <https://GeographicLib.SourceForge.io/html/TransverseMercatorExact_8hpp_source.html>} file C{Header}.
Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2017) and licensed under the MIT/X11 License. For more information, see the U{GeographicLib<https://GeographicLib.SourceForge.io/>} documentation.
The method entails using the C{Thompson Transverse Mercator} as an intermediate projection. The projections from the intermediate coordinates to C{phi, lam} and C{x, y} are given by elliptic functions. The inverse of these projections are found by Newton's method with a suitable starting guess.
The relevant section of L.P. Lee's paper U{Conformal Projections Based On Jacobian Elliptic Functions<https://DOI.org/10.3138/X687-1574-4325-WM62>} is part V, pp 67--101. The C++ implementation and notation closely follow Lee, with the following exceptions::
Lee here Description
x/a xi Northing (unit Earth)
y/a eta Easting (unit Earth)
s/a sigma xi + i * eta
y x Easting
x y Northing
k e Eccentricity
k^2 mu Elliptic function parameter
k'^2 mv Elliptic function complementary parameter
m k Scale
zeta zeta Complex longitude = Mercator = chi in paper
s sigma Complex GK = zeta in paper
Minor alterations have been made in some of Lee's expressions in an attempt to control round-off. For example, C{atanh(sin(phi))} is replaced by C{asinh(tan(phi))} which maintains accuracy near C{phi = pi/2}. Such changes are noted in the code. ''' # make sure int/int division yields float quotient raise ImportError('%s 1/2 == %d' % ('division', division))
_NamedBase, _xnamed _TypeError _toXtm8, _to7zBlldfn
fmod, radians, sinh, sqrt, tan
except ImportError: # no geographiclib
from pygeodesy.fmath import NAN as _NAN
def _diff182(deg0, deg): # mimick Math.AngDiff '''Compute C{deg - deg0}, reduced to C{[-180,180]} accurately. ''' d, t = _sum2(_wrap180(-deg0), _wrap180(deg)) d = _wrap180(d) if d == 180 and t > 0: d = -180 return _sum2(d, t)
def _fix90(deg): # mimick Math.LatFix '''Replace angles outside [-90,90] by NaN. ''' return _NAN if abs(deg) > 90 else deg
def _sum2(u, v): # mimick Math::sum, actually sum2 '''Error free transformation of a C{sum}.
@return: 2-Tuple C{(sum_u_plus_v, residual)}.
@note: The C{residual} can be the same as one of the first two arguments. ''' s = u + v r = s - v t = s - r r -= u t -= v t = -(r + t) # u + v = s + t # = round(u + v) + t return s, t
def _wrap180(deg): # mimick Math.AngNormalize '''Reduce angle to (-180,180] ''' # with Python 2.7.16 and 3.7.3 on macOS 10.13.6 # fmod( 0, 360) == 0.0 # fmod( 360, 360) == 0.0 # fmod(-0, 360) == 0.0 # fmod(-0.0, 360) == -0.0 # fmod(-360, 360) == -0.0 # however, using the % operator ... # 0 % 360 == 0 # 360 % 360 == 0 # 360.0 % 360 == 0.0 # -0 % 360 == 0 # -360 % 360 == 0 # -0.0 % 360 == 0.0 # -360.0 % 360 == 0.0
# On Windows 32-bit with Python 2.7, math.fmod(-0.0, 360) # == +0.0. This fixes this bug. See also Math::AngNormalize # in the C++ library. Math::sincosd has a similar fix. d = fmod(deg, 360) if deg else deg return (d + 360) if d <= -180 else ( d if d <= 180 else (d - 360))
'''Exact Transverse Mercator (ETM) parse, projection or other L{Etm} issue. '''
'''Exact Transverse Mercator (ETM) coordinate, a sub-class of L{Utm}, a Universal Transverse Mercator (UTM) coordinate using the L{ExactTransverseMercator} projection for highest accuracy.
@note: Conversion of L{Etm} coordinates to and from (geodetic) lat- and longitude is 3-4 times slower than L{Utm}.
@see: Karney's U{Detailed Description<https://GeographicLib.SourceForge.io/ html/classGeographicLib_1_1TransverseMercatorExact.html#details>}. '''
datum=Datums.WGS84, falsed=True, convergence=None, scale=None, name=''): '''New L{Etm} coordinate.
@param zone: Longitudinal UTM zone (C{int}, 1..60) or zone with/-out (latitudinal) Band letter (C{str}, '01C'..'60X'). @param hemisphere: Northern or southern hemisphere (C{str}, C{'N[orth]'} or C{'S[outh]'}). @param easting: Easting, see B{C{falsed}} (C{meter}). @param northing: Northing, see B{C{falsed}} (C{meter}). @keyword band: Optional, (latitudinal) band (C{str}, 'C'..'X'). @keyword datum: Optional, this coordinate's datum (L{Datum}). @keyword falsed: Both B{C{easting}} and B{C{northing}} are falsed (C{bool}). @keyword convergence: Optional meridian convergence, bearing off grid North, clockwise from true North (C{degrees}) or C{None}. @keyword scale: Optional grid scale factor (C{scalar}) or C{None}. @keyword name: Optional name (C{str}).
@raise EllipticError: No convergence.
@raise ETMError: Invalid B{C{zone}}, B{C{hemishere}} or B{C{band}}.
@example:
>>> import pygeodesy >>> u = pygeodesy.Etm(31, 'N', 448251, 5411932) ''' band=band, datum=datum, falsed=falsed, convergence=convergence, scale=scale, name=name)
def exactTM(self): '''Get the ETM projection (L{ExactTransverseMercator}). '''
def exactTM(self, exactTM): '''Set the ETM projection (L{ExactTransverseMercator}). '''
or exactTM.flattening != E.f: raise ETMError('%r vs %r' % (exactTM, E))
'''Parse a string to a ETM coordinate.
@return: The coordinate (L{Etm}).
@see: Function L{parseETM5} in this module L{etm}. ''' return parseETM5(strETM, datum=self.datum, Etm=self.classof)
'''Convert this ETM coordinate to an (ellipsoidal) geodetic point.
@keyword LatLon: Optional, ellipsoidal (sub-)class to return the point (C{LatLon}) or C{None}. @keyword unfalse: Unfalse B{C{easting}} and B{C{northing}} if falsed (C{bool}).
@return: This ETM coordinate as (B{C{LatLon}}) or a L{LatLonDatum5Tuple}C{(lat, lon, datum, convergence, scale)} if B{C{LatLon}} is C{None}.
@raise EllipticError: No convergence.
@raise TypeError: If B{C{LatLon}} is not ellipsoidal.
@example:
>>> from pygeodesy import ellipsoidalVincenty as eV, Etm >>> u = Etm(31, 'N', 448251.795, 5411932.678) >>> ll = u.toLatLon(eV.LatLon) # 48°51′29.52″N, 002°17′40.20″E ''' # double check that this and exactTM's ellipsoids stil match raise ETMError('%r vs %r' % (xTM._E, d.ellipsoid))
return self._latlon5(LatLon)
# f = unfalse == self.falsed # == unfalse and self.falsed or (not unfalse and not self.falsed) # == unfalse if self.falsed else not unfalse # == unfalse if self.falsed else f f = unfalse
'''(INTERNAL) See C{.toLatLon}, C{toEtm8}, C{_toXtm8}. '''
'''Coopy this ETM to a UTM coordinate.
@return: The UTM coordinate (L{Utm}). ''' return self._xnamed(self._xcopy2(Utm))
'''A Python version of Karney's U{TransverseMercatorExact <https://GeographicLib.SourceForge.io/html/TransverseMercatorExact_8cpp_source.html>} C++ class, a numerically exact transverse mercator projection, referred to as C{TMExact} here.
@see: C{TMExact(real a, real f, real k0, bool extendp)}. '''
'''New L{ExactTransverseMercator} projection.
@keyword datum: The datum and ellipsoid to use (C{Datum}). @keyword lon0: The central meridian (C{degrees180}). @keyword k0: The central scale factor (C{float}). @keyword extendp: Use the extended domain (C{bool}). @keyword name: Optional name for the projection (C{str}).
@raise EllipticError: No convergence.
@raise ETMError: Invalid B{C{k0}}.
@raise TypeError: Invalid B{C{datum}}.
@note: The maximum error for all 255.5K U{TMcoords.dat <https://Zenodo.org/record/32470>} tests (with C{0 <= lat <= 84} and C{0 <= lon}) is C{5.2e-08 .forward} or 52 nano-meter easting and northing and C{3.8e-13 .reverse} or 0.38 pico-degrees lat- and longitude (with Python 3.7.3, 2.7.16, PyPy6 3.5.3 and PyPy6 2.7.13, all in 64-bit on macOS 10.13.6 High Sierra). ''' self.name = name
def datum(self): '''Get the datum (L{Datum}) or C{None}. '''
def datum(self, datum): '''Set the datum and ellipsoid (L{Datum}).
@raise EllipticError: No convergence.
@raise TypeError: Invalid B{C{datum}}. '''
def extendp(self): '''Get using the extended domain (C{bool}). '''
def flattening(self): '''Get the flattening (C{float}). '''
'''Forward projection, from geographic to transverse Mercator.
@param lat: Latitude of point (C{degrees}). @param lon: Longitude of point (C{degrees}). @keyword lon0: Central meridian of the projection (C{degrees}).
@return: L{EasNorExact4Tuple}C{(easting, northing, convergence, scale)} in C{meter}, C{meter}, C{degrees} and C{scalar}.
@see: C{void TMExact::Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k)}.
@raise EllipticError: No convergence. ''' # Explicitly enforce the parity if lat < 0: _lat, lat = True, -lat if lon < 0: _lon, lon = True, -lon if lon > 90: backside = True if lat == 0: _lat = True lon = 180 - lon
# u,v = coordinates for the Thompson TM, Lee 54 u, v = self._Eu_K, 0 u, v = 0, self._Ev_K else: # tau = tan(phi), taup = sinh(psi)
xi = 2 * self._Eu_E - xi
g, k = lon, self._k0 else: # Recompute (T, L) from (u, v) to improve accuracy of Scale
g = 180 - g y, g = -y, -g x, g = -x, -g
def k0(self): '''Get the central scale factor (C{float}), aka I{C{scale0}}. '''
def k0(self, k0): '''Set the central scale factor (C{float}), aka I{C{scale0}}.
@raise EllipticError: Invalid B{C{k0}}. ''' raise ETMError('%s invalid: %r' % ('k0', k0))
def lon0(self): '''Get the central meridian (C{degrees180}). ''' return self._lon0
def lon0(self, lon0): '''Set the central meridian (C{degrees180}). '''
def majoradius(self): '''Get the major (equatorial) radius, semi-axis (C{float}). '''
'''(INTERNAL) Get elliptic functions and pre-compute frequently used values.
@raise EllipticError: No convergence. ''' # assert e2 == e**2
'''Reverse projection, from transverse Mercator to geographic.
@param x: Easting of point (C{meters}). @param y: Northing of point (C{meters}). @keyword lon0: Central meridian of the projection (C{degrees}).
@return: L{LatLonExact4Tuple}C{(lat, lon, convergence, scale)} in C{degrees}, C{degrees180}, C{degrees} and C{scalar}.
@see: C{void TMExact::Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k)}
@raise EllipticError: No convergence. ''' # undoes the steps in .forward. if y < 0: _lat, xi = True, -xi if x < 0: _lon, eta = True, -eta if xi > self._Eu_E: backside = True xi = 2 * self._Eu_E - xi
# u,v = coordinates for the Thompson TM, Lee 54 u, v = 0, self._Ev_K else:
else: lat, lon = 90, 0 g, k = 0, self._k0
lon, g = 180 - lon, 180 - g lat, g = -lat, -g lon, g = -lon, -g
''' @note: Argument B{C{d2}} is C{_mu * cnu**2 + _mv * cnv**2} from C{._sigma3} or C{._zeta3}.
@return: 2-Tuple C{(convergence, scale)}.
@see: C{void TMExact::Scale(real tau, real /*lam*/, real snu, real cnu, real dnu, real snv, real cnv, real dnv, real &gamma, real &k)}. ''' # Lee 55.12 -- negated for our sign convention. g gives # the bearing (clockwise from true north) of grid north # Lee 55.13 with nu given by Lee 9.1 -- in sqrt change # the numerator from # # (1 - snu^2 * dnv^2) to (_mv * snv^2 + cnu^2 * dnv^2) # # to maintain accuracy near phi = 90 and change the # denomintor from # (dnu^2 + dnv^2 - 1) to (_mu * cnu^2 + _mv * cnv^2) # # to maintain accuracy near phi = 0, lam = 90 * (1 - e). # Similarly rewrite sqrt term in 9.1 as # # _mv + _mu * c^2 instead of 1 - _mu * sin(phi)^2
''' @return: 3-Tuple C{(xi, eta, d2)}.
@see: C{void TMExact::sigma(real /*u*/, real snu, real cnu, real dnu, real v, real snv, real cnv, real dnv, real &xi, real &eta)}.
@raise EllipticError: No convergence. ''' # Lee 55.4 writing # dnu^2 + dnv^2 - 1 = _mu * cnu^2 + _mv * cnv^2
''' @return: 2-Tuple C{(du, dv)}.
@see: C{void TMExact::dwdsigma(real /*u*/, real snu, real cnu, real dnu, real /*v*/, real snv, real cnv, real dnv, real &du, real &dv)}. ''' # Reciprocal of 55.9: dw / ds = dn(w)^2/_mv, # expanding complex dn(w) using A+S 16.21.4
'''Invert C{sigma} using Newton's method.
@return: 2-Tuple C{(u, v)}.
@see: C{void TMExact::sigmainv(real xi, real eta, real &u, real &v)}.
@raise EllipticError: No convergence. ''' # min iterations = 2, max = 7, mean = 3.9 else: raise EllipticError('no %s convergence' % ('sigmaInv',))
'''Starting point for C{sigmaInv}.
@return: 3-Tuple C{(u, v, trip)}.
@see: C{bool TMExact::sigmainv0(real xi, real eta, real &u, real &v)}. ''' eta - self._Ev_KE): # sigma as a simple pole at # w = w0 = Eu.K() + i * Ev.K() # and sigma is approximated by # sigma = (Eu.E() + i * Ev.KE()) + 1/(w - w0)
xi < self._Eu_E_1_4): # At w = w0 = i * Ev.K(), we have # sigma = sigma0 = i * Ev.KE() # sigma' = sigma'' = 0 # # including the next term in the Taylor series gives: # sigma = sigma0 - _mv / 3 * (w - w0)^3 # # When inverting this, we map arg(w - w0) = [-pi/2, -pi/6] # to arg(sigma - sigma0) = [-pi/2, pi/2] # mapping arg = [-pi/2, -pi/6] to [-pi/2, pi/2]
# Error using this guess is about 0.068 * rad^(5/3) # Map the range [-90, 180] in sigma space to [-90, 0] in # w space. See discussion in zetainv0 on the cut for ang.
else: # use w = sigma * Eu.K/Eu.E (correct in the limit _e -> 0)
'''Return a C{str} representation.
@param kwds: Optional, keyword arguments. ''' d = dict(datum=self.datum.name, lon0=self.lon0, k0=self.k0, extendp=self.extendp) if self.name: d['name'] = self.name if kwds: d.update(kwds) return ', '.join('%s=%s' % t for t in sorted(d.items()))
''' @return: 3-Tuple C{(taup, lambda, d2)}.
@see: C{void TMExact::zeta(real /*u*/, real snu, real cnu, real dnu, real /*v*/, real snv, real cnv, real dnv, real &taup, real &lam)} ''' # Lee 54.17 but write # atanh(snu * dnv) = asinh(snu * dnv / sqrt(cnu^2 + _mv * snu^2 * snv^2)) # atanh(_e * snu / dnv) = asinh(_e * snu / sqrt(_mu * cnu^2 + _mv * cnv^2)) # Overflow value s.t. atan(overflow) = pi/2 # psi = asinh(t1) - asinh(t2) # taup = sinh(psi) atan2(e * cnu * snv, dnu * cnv)) if (d1 > 0 and d2 > 0) else 0
''' @return: 2-Tuple C{(du, dv)}.
@see: C{void TMExact::dwdzeta(real /*u*/, real snu, real cnu, real dnu, real /*v*/, real snv, real cnv, real dnv, real &du, real &dv)}. ''' # Lee 54.21 but write # (1 - dnu^2 * snv^2) = (cnv^2 + _mu * snu^2 * snv^2) # (see A+S 16.21.4)
'''Invert C{zeta} using Newton's method.
@return: 2-Tuple C{(u, v)}.
@see: C{void TMExact::zetainv(real taup, real lam, real &u, real &v)}.
@raise EllipticError: No convergence. ''' # min iterations = 2, max = 6, mean = 4.0 else: raise EllipticError('no %s convergence' % ('zetaInv',))
'''Starting point for C{zetaInv}.
@return: 3-Tuple C{(u, v, trip)}.
@see: C{bool TMExact::zetainv0(real psi, real lam, # radians real &u, real &v)}. ''' and psi < lam - self._1_e_PI_2): # N.B. this branch is normally not taken because psi < 0 # is converted psi > 0 by Forward. # # There's a log singularity at w = w0 = Eu.K() + i * Ev.K(), # corresponding to the south pole, where we have, approximately # # psi = _e + i * pi/2 - _e * atanh(cos(i * (w - w0)/(1 + _mu/2))) # # Inverting this gives: h = sinh(1 - psi / self._e) a = (PI_2 - lam) / self._e s, c = sincos2(a) u = self._Eu_K - asinh(s / hypot(c, h)) * self._mu_2_1 v = self._Ev_K - atan2(c, h) * self._mu_2_1
# At w = w0 = i * Ev.K(), we have # # zeta = zeta0 = i * (1 - _e) * pi/2 # zeta' = zeta'' = 0 # # including the next term in the Taylor series gives: # # zeta = zeta0 - (_mv * _e) / 3 * (w - w0)^3 # # When inverting this, we map arg(w - w0) = [-90, 0] to # arg(zeta - zeta0) = [-90, 180]
# Error using this guess is about 0.21 * (rad/e)^(5/3) # atan2(dlam-psi, psi+dlam) + 45d gives arg(zeta - zeta0) # in range [-135, 225). Subtracting 180 (since multiplier # is negative) makes range [-315, 45). Multiplying by 1/3 # (for cube root) gives range [-105, 15). In particular # the range [-90, 180] in zeta space maps to [-90, 0] in # w space as required.
else: # Use spherical TM, Lee 12.6 -- writing C{atanh(sin(lam) / # cosh(psi)) = asinh(sin(lam) / hypot(cos(lam), sinh(psi)))}. # This takes care of the log singularity at C{zeta = Eu.K()}, # corresponding to the north pole. # But scale to put 90, 0 on the right place
'''Parse a string representing a UTM coordinate, consisting of C{"zone[band] hemisphere easting northing"}.
@param strUTM: A UTM coordinate (C{str}). @keyword datum: Optional datum to use (L{Datum}). @keyword Etm: Optional (sub-)class to return the UTM coordinate (L{Etm}) or C{None}. @keyword falsed: Both easting and northing are falsed (C{bool}). @keyword name: Optional B{C{Etm}} name (C{str}).
@return: The UTM coordinate (B{C{Etm}}) or a L{UtmUps5Tuple}C{(zone, hemipole, easting, northing, band)} if B{C{Etm}} is C{None}. The C{hemipole} is the hemisphere C{'N'|'S'}.
@raise ETMError: Invalid B{C{strUTM}}.
@example:
>>> u = parseETM5('31 N 448251 5411932') >>> u.toStr2() # [Z:31, H:N, E:448251, N:5411932] >>> u = parseETM5('31 N 448251.8 5411932.7') >>> u.toStr() # 31 N 448252 5411933 ''' r = _parseUTM5(strUTM, ETMError) if Etm is not None: z, h, e, n, B = r r = Etm(z, h, e, n, band=B, datum=datum, falsed=falsed) return _xnamed(r, name)
zone=None, **cmoff): '''Convert a lat-/longitude point to an ETM coordinate.
@param latlon: Latitude (C{degrees}) or an (ellipsoidal) geodetic C{LatLon} point. @keyword lon: Optional longitude (C{degrees}) or C{None}. @keyword datum: Optional datum for this ETM coordinate, overriding B{C{latlon}}'s datum (C{Datum}). @keyword Etm: Optional (sub-)class to return the ETM coordinate (L{Etm}) or C{None}. @keyword falsed: False both easting and northing (C{bool}). @keyword name: Optional B{C{Utm}} name (C{str}). @keyword zone: Optional UTM zone to enforce (C{int} or C{str}). @keyword cmoff: DEPRECATED, use B{C{falsed}}. Offset longitude from the zone's central meridian (C{bool}).
@return: The ETM coordinate (B{C{Etm}}) or a L{UtmUps8Tuple}C{(zone, hemipole, easting, northing, band, datum, convergence, scale)} if B{C{Etm}} is C{None} or not B{C{falsed}}. The C{hemipole} is the C{'N'|'S'} hemisphere.
@raise EllipticError: No convergence.
@raise ETMError: Invalid B{C{zone}}.
@raise TypeError: If B{C{latlon}} is not ellipsoidal.
@raise RangeError: If B{C{lat}} outside the valid UTM bands or if B{C{lat}} or B{C{lon}} outside the valid range and L{rangerrors} set to C{True}.
@raise ValueError: If B{C{lon}} value is missing or if B{C{latlon}} is invalid. ''' falsed, name, zone, ETMError, **cmoff)
# **) MIT License # # Copyright (C) 2016-2020 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE. |