Coverage for pygeodesy/ecef.py: 95%
622 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-10 16:16 -0500
2# -*- coding: utf-8 -*-
4u'''I{Geocentric} Earth-Centered, Earth-Fixed (ECEF) coordinates.
6Geocentric conversions transcoded from I{Charles Karney}'s C++ class U{Geocentric
7<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Geocentric.html>} into pure
8Python class L{EcefKarney}, class L{EcefFukushima} from I{Toshio Fukushima}'s U{Fortran
9<https://www.ResearchGate.net/publication/277721539>} version, class L{EcefOlson} from I{D.K. Olson}'s
10U{paper<https://IEEEXplore.IEEE.org/document/481290>}, class L{EcefOsen} with I{Karl Osen}'s U{improvements
11<https://HAL.Science/hal-01704943v2>} of I{Jijie Zhu}'s U{algorithm <https://IEEEXplore.IEEE.org/document/303772>},
12class L{EcefSudano} based on I{John Sudano}'s U{paper<https://www.ResearchGate.net/publication/3709199>},
13class L{EcefUPC} using the I{Universitat Politècnica de Catalunya}'s U{method, page 186
14<https://GSSC.ESA.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf>}, class L{EcefVeness} transcoded
15from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal, Cartesian
16<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}, class L{EcefYou}
17implementing I{Rey-Jer You}'s U{transformations<https://www.ResearchGate.net/publication/240359424>}
18and classes L{EcefFarrell22} and L{EcefFarrell22} from I{Jay A. Farrell}'s U{Table 2.1 and 2.2
19<https://Books.Google.com/books?id=fW4foWASY6wC>}, page 29-30.
21Following is a copy of I{Karney}'s U{Detailed Description
22<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Geocentric.html>}.
24Convert between geodetic coordinates C{lat}-, C{lon}gitude and height C{h} (measured vertically
25from the surface of the ellipsoid) to geocentric C{x}, C{y} and C{z} coordinates, also known as
26I{Earth-Centered, Earth-Fixed} (U{ECEF<https://WikiPedia.org/wiki/ECEF>}).
28The origin of geocentric coordinates is at the center of the earth. The C{z} axis goes thru
29the North pole, C{lat} = 90°. The C{x} axis goes thru C{lat} = 0°, C{lon} = 0°.
31The I{local (cartesian) origin} is at (C{lat0}, C{lon0}, C{height0}). The I{local} C{x} axis points
32East, the I{local} C{y} axis points North and the I{local} C{z} axis is normal to the ellipsoid. The
33plane C{z = -height0} is tangent to the ellipsoid, hence the alternate name I{local tangent plane}.
35Forward conversion from geodetic to geocentric (ECEF) coordinates is straightforward.
37For the reverse transformation we use Hugues Vermeille's U{Direct transformation from geocentric
38coordinates to geodetic coordinates<https://DOI.org/10.1007/s00190-002-0273-6>}, J. Geodesy
39(2002) 76, page 451-454.
41Several changes have been made to ensure that the method returns accurate results for all finite
42inputs (even if h is infinite). The changes are described in Appendix B of C. F. F. Karney
43U{Geodesics on an ellipsoid of revolution<https://ArXiv.org/abs/1102.1215v1>}, Feb. 2011, 85,
44105-117 (U{preprint<https://ArXiv.org/abs/1102.1215v1>}). Vermeille similarly updated his method
45in U{An analytical method to transform geocentric into geodetic coordinates
46<https://DOI.org/10.1007/s00190-010-0419-x>}, J. Geodesy (2011) 85, page 105-117. See U{Geocentric
47coordinates<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>} for more information.
49The errors in these routines are close to round-off. Specifically, for points within 5,000 Km of
50the surface of the ellipsoid (either inside or outside the ellipsoid), the error is bounded by 7
51nm (7 nanometers) for the WGS84 ellipsoid. See U{Geocentric coordinates
52<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>} for further information on the errors.
54@note: The C{reverse} methods of all C{Ecef...} classes return by default C{INT0} as the (geodetic)
55longitude for I{polar} ECEF location C{x == y == 0}. Use keyword argument C{lon00} or property
56C{lon00} to configure that value.
58@see: Module L{ltp} and class L{LocalCartesian}, a transcription of I{Charles Karney}'s C++ class
59U{LocalCartesian<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1LocalCartesian.html>},
60for conversion between geodetic and I{local cartesian} coordinates in a I{local tangent
61plane} as opposed to I{geocentric} (ECEF) ones.
62'''
64from pygeodesy.basics import copysign0, _isin, isscalar, issubclassof, neg, map1, \
65 _xinstanceof, _xsubclassof, typename # _args_kwds_names
66from pygeodesy.constants import EPS, EPS0, EPS02, EPS1, INT0, PI, PI_2, _0_0, _0_25, \
67 _0_5, _1_0, _1_0_1T, _1_5, _2_0, _3_0, _4_0, _6_0, \
68 _copysign_1_0, _isNAN0, _over, isnon0 # PYCHOK used!
69from pygeodesy.datums import _ellipsoidal_datum, _WGS84, a_f2Tuple, _EWGS84
70from pygeodesy.ecefLocals import _EcefLocal
71# from pygeodesy.ellipsoids import a_f2Tuple, _EWGS84 # from .datums
72from pygeodesy.errors import _IndexError, LenError, _ValueError, _TypesError, \
73 _xattr, _xdatum, _xkwds
74from pygeodesy.fmath import cbrt, _fdotf, hypot, hypot1, hypot2_
75from pygeodesy.fsums import Fsum, fsumf_, Fmt, unstr
76# from pygeodesy.internals import typename # from .basics
77from pygeodesy.interns import NN, _a_, _C_, _datum_, _ellipsoid_, _f_, _height_, \
78 _lat_, _lon_, _M_, _name_, _singular_, _SPACE_, \
79 _x_, _xyz_, _y_, _z_
80from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
81from pygeodesy.named import _name__, _name1__, _NamedBase, _NamedTuple, _Pass, _xnamed
82from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \
83 PhiLam2Tuple, Vector3Tuple, Vector4Tuple
84from pygeodesy.props import deprecated_method, deprecated_property, Property_RO, \
85 property_RO, property_ROver
86# from pygeodesy.streprs import Fmt, unstr # from .fsums
87from pygeodesy.units import _isRadius, Degrees_, Height, Int, Lam, Lat, Lon, \
88 Meter, Phi, Scalar, Scalar_
89from pygeodesy.utily import atan1, atan1d, atan2, atan2d, degrees90, degrees180, \
90 sincos2, sincos2_, sincos2d_
91# from pygeodesy.vector3d import Vector3d # _MODS
93from math import acos, asin, cos, degrees, fabs, radians, sqrt
95__all__ = _ALL_LAZY.ecef
96__version__ = '26.09.28'
98_Ecef_ = 'Ecef'
99_prolate_ = 'prolate'
100_TOL = 1.e-12 # degrees > 1.e-14
101_TRIPS = 33 # 8..9 sufficient
102_xyz_y_z = _xyz_, _y_, _z_ # _args_kwds_names(_xyzn4)[:3]
105def _Degrees2Radians(tol): # for EcefUPC
106 return Degrees_(tol=tol, low=EPS, Error=EcefError).toRadians()
109class _EcefBase(_NamedBase):
110 '''(INTERNAL) Base class for C{Ecef*} convertor classes.
111 '''
112 _datum = _WGS84
113 _E = _EWGS84
114 _e_e2 = False # or True or 2-tuple for EcefYou only
115 _lon00 = INT0 # arbitrary, "polar" lon for LocalCartesian, Ltp
117 def __init__(self, a_ellipsoid=_EWGS84, f=None, lon00=INT0, **name):
118 '''New C{Ecef*} converter.
120 @arg a_ellipsoid: An ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple})
121 or a datum (L{Datum}) or the ellipsoid's equatorial
122 radius (C{meter}).
123 @kwarg f: C{None} or the ellipsoid flattening (C{scalar}), required if
124 C{B{a_ellipsoid} is scalar}.
125 @kwarg lon00: An arbitrary, I{"polar"} longitude (C{degrees}), see the
126 C{reverse} method.
127 @kwarg name: Optional C{B{name}=NN} (C{str}).
129 @raise EcefError: If B{C{a_ellipsoid}} is not an L{Ellipsoid}, L{Ellipsoid2},
130 L{a_f2Tuple} or L{Datum} instance or a positive C{scalar}
131 or if B{C{f}} is not C{scalar} and less than C{1.0}.
132 '''
133 try:
134 E = a_ellipsoid
135 if f is None:
136 pass
137 elif _isRadius(E) and isscalar(f):
138 E = a_f2Tuple(E, f)
139 else:
140 raise ValueError() # _invalid_
142 if not _isin(E, _EWGS84, _WGS84):
143 d = _ellipsoidal_datum(E, **name)
144 E = d.ellipsoid
145 if E.a < EPS or E.f > EPS1:
146 raise ValueError() # _invalid_
147 self._datum = d
148 self._E = E
150 if self._e_e2: # EcefYou
151 E = self.ellipsoid
152 e2 = E.a2 - E.b2
153 if e2 < 0 or E.f < 0:
154 raise ValueError(_prolate_)
155 self._e_e2 = sqrt(e2), e2
157 except (TypeError, ValueError) as x:
158 t = unstr(self.classname, a=a_ellipsoid, f=f)
159 raise EcefError(_SPACE_(t, _ellipsoid_), cause=x)
161 if name:
162 self.name = name
163 if lon00 is not INT0:
164 self.lon00 = lon00
166 def __eq__(self, other):
167 '''Compare this and an other Ecef.
169 @arg other: The other ecef (C{Ecef*}).
171 @return: C{True} if equal, C{False} otherwise.
172 '''
173 return other is self or (isinstance(other, type(self)) and
174 other.ellipsoid == self.ellipsoid)
176 @Property_RO
177 def datum(self):
178 '''Get the datum (L{Datum}).
179 '''
180 return self._datum
182 @Property_RO
183 def ellipsoid(self):
184 '''Get the ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
185 '''
186 return self._E
188 @Property_RO
189 def equatoradius(self):
190 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}).
191 '''
192 return self.ellipsoid.a
194 a = equatorialRadius = equatoradius # Karney property
196 @Property_RO
197 def flattening(self): # Karney property
198 '''Get the C{ellipsoid}'s flattening (C{scalar}), positive for
199 I{oblate}, negative for I{prolate} or C{0} for I{near-spherical}.
200 '''
201 return self.ellipsoid.f
203 f = flattening
205 def _forward(self, lat, lon, h, name, M=False, _philam=False): # in .ltp.LocalCartesian.forward and -.reset
206 '''(INTERNAL) Common for all C{Ecef*}.
208 @note: From C{Karney}'s, let C{v} be a unit vector located at C{(lat,
209 lon, h)}. We can express C{v} as column vectors in one of two
210 ways, C{v1} in East, North, Up (ENU) coordinates (where the
211 components are relative to a local coordinate system at C{C(lat0,
212 lon0, h0)}) or as C{v0} in geocentric C{x, y, z} coordinates.
213 Then, M{v0 = M ⋅ v1} where C{M} is the rotation matrix.
214 '''
215 if _philam: # lat, lon in radians
216 sa, ca, sb, cb = sincos2_(lat, lon)
217 lat = Lat(degrees90( lat), Error=EcefError)
218 lon = Lon(degrees180(lon), Error=EcefError)
219 else:
220 sa, ca, sb, cb = sincos2d_(lat, lon)
222 E = self.ellipsoid
223 n = E.roc1_(sa, ca) if self._e_e2 else E.roc1_(sa)
224 H = _isNAN0(h)
225 c = (H + n) * ca
226 x = cb * c
227 y = sb * c
228 z = (H + n * E.e21) * sa
230 m = self._Matrix(sa, ca, sb, cb) if M else None
231 n = self._name__(name)
232 return Ecef9Tuple(x, y, z, lat, lon, h, 0, # C=0, forward
233 m, self.datum, name=n)
235 def forward(self, latlonh, lon=None, height=0, M=False, **name):
236 '''Convert from geodetic C{(lat, lon, height)} to geocentric C{(x, y, z)}.
238 @arg latlonh: Either a C{LatLon}, an L{Ecef9Tuple} or C{scalar}
239 latitude (C{degrees}).
240 @kwarg lon: Optional C{scalar} longitude for C{scalar} B{C{latlonh}}
241 (C{degrees}).
242 @kwarg height: Optional height (C{meter}), vertically above (or below)
243 the surface of the ellipsoid.
244 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
245 @kwarg name: Optional C{B{name}=NN} (C{str}).
247 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
248 geocentric C{(x, y, z)} coordinates for the given geodetic ones
249 C{(lat, lon, height)}, case C{C} (0, forward), rotation matrix
250 C{M} (L{EcefMatrix} or C{None}) and C{datum}.
252 @raise EcefError: If B{C{latlonh}} not C{LatLon}, L{Ecef9Tuple} or
253 C{scalar} or B{C{lon}} not C{scalar} for C{scalar}
254 B{C{latlonh}} or C{abs(lat)} exceeds 90°.
256 @note: Use method C{.forward_} to specify C{lat} and C{lon} in C{radians}
257 and avoid double angle conversions.
258 '''
259 llhn = _llhn4(latlonh, lon, height, **name)
260 return self._forward(*llhn, M=M)
262 def forward_(self, phi, lam, height=0, M=False, **name):
263 '''Like method C{.forward} except with geodetic lat- and longitude given
264 in I{radians}.
266 @arg phi: Latitude in I{radians} (C{scalar}).
267 @arg lam: Longitude in I{radians} (C{scalar}).
268 @kwarg height: Optional height (C{meter}), vertically above (or below)
269 the surface of the ellipsoid.
270 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
271 @kwarg name: Optional C{B{name}=NN} (C{str}).
273 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
274 with C{lat} set to C{degrees90(B{phi})} and C{lon} to
275 C{degrees180(B{lam})}.
277 @raise EcefError: If B{C{phi}} or B{C{lam}} invalid or not C{scalar}.
278 '''
279 try: # like function C{_llhn4} below
280 plhn = Phi(phi), Lam(lam), Height(height), _name__(name)
281 except (TypeError, ValueError) as x:
282 raise EcefError(phi=phi, lam=lam, height=height, cause=x)
283 return self._forward(*plhn, M=M, _philam=True)
285 @property_ROver
286 def _Geocentrics(self):
287 '''(INTERNAL) Get the valid geocentric classes. I{once}.
288 '''
289 return (Ecef9Tuple, # overwrite property_ROver
290 _MODS.vector3d.Vector3d) # _MODS.cartesianBase.CartesianBase
292 @property
293 def lon00(self):
294 '''Get the I{"polar"} longitude (C{degrees}), see method C{reverse}.
295 '''
296 return self._lon00
298 @lon00.setter # PYCHOK setter!
299 def lon00(self, lon00):
300 '''Set the I{"polar"} longitude (C{degrees}), see method C{reverse}.
301 '''
302 self._lon00 = _Lon00(lon00)
304 def _Matrix(self, *sa_ca_sb_cb):
305 '''(INTERNAL) Create a rotation L{EcefMatrix}.
306 '''
307 return self._xnamed(EcefMatrix(*sa_ca_sb_cb))
309 def _polon(self, y, x, p, lon00=None, name=NN): # PYCHOK name unused
310 '''(INTERNAL) Handle I{"polar"} longitude.
311 '''
312 return atan2d(y, x) if p else (self.lon00 if lon00 is None else _Lon00(lon00))
314 def reverse(self, xyz, y=None, z=None, M=False, **lon00_name):
315 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}.
317 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF X
318 coordinate (C{meter}).
319 @kwarg y: ECEF Y coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
320 @kwarg z: ECEF Z coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
321 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
322 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument
323 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude
324 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.
326 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
327 geodetic coordinates C{(lat, lon, height)} for the given geocentric
328 ones C{(x, y, z)}, case indicator C{C} (C{int} 1..6), rotation matrix
329 C{M} (L{EcefMatrix} or C{None}) and C{datum}.
331 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
332 not C{scalar} for C{scalar} B{C{xyz}}.
333 '''
334 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, **lon00_name)
336 E = self.ellipsoid
337 i = None
339 sa, ca, sb, cb, h, p, C = _norm7(y, x, z, E)
340 if C: # PYCHOK no cover
341 pass # too high, too far
343 elif p < EPS: # near polar
344 p = 0 # force lon00
345 sa = _copysign_1_0(z)
346 ca = _0_0
347 h = fabs(z) - E.b
348 C = 2 # polar
350 elif E.e4: # E.isEllipsoidal
351 s, t = _equatorial2(E, p, z)
352 r = fsumf_(s, t, -E.e4)
353 if t or r > 0:
354 try:
355 sa, ca, h, i, C = self._reverse5(1, h, p, z, y, x, r, s, t) # PYLINT E1111
356 except (TypeError, ValueError) as X:
357 t = unstr(self.reverse, x=x, y=y, z=z)
358 raise EcefError(t, cause=X)
360 else: # near equatorial plane: e = E.e4 * q == 0 and r <= 0
361 # This leads to k = 0 (oblate, equatorial plane) and k + E.e^2 = 0
362 # (prolate, rotation axis) and the generation of 0/0 in the general
363 # formulas for phi and h, using the general formula and division
364 # by 0 in formula for h. Handle this case by taking the limits:
365 # f > 0: z -> 0, k -> E.e2 * sqrt(q) / sqrt(E.e4 - s)
366 # f < 0: r -> 0, k + E.e2 -> -E.e2 * sqrt(q) / sqrt(E.e4 - s)
367 sa, ca, h, C = _equatorial4(E, s, z)
369 else: # E.isSpherical: E.e4 == 0
370 # Dealing with underflow in the general case with E.e2 = 0 is
371 # difficult. Origin maps to North pole, same as with ellipsoid.
372 sa, ca, _ = _norm3((z if h else _1_0), p)
373 h -= E.a
374 C = 4 # spherical
376 lat = atan1d(sa, ca)
377 # lon00 <https://GitHub.com/mrJean1/PyGeodesy/issues/77>
378 lon = self._polon(sb, cb, p, **lon00_name)
379 m = self._Matrix(sa, ca, sb, cb) if M else None
380 return Ecef9Tuple(x, y, z, lat, lon, h, C, m, self.datum,
381 iteration=i, name=self._name__(name)) # PYCHOK return
383 def _reverse5(self, *C_h_p_z_y_x_r_s_t): # PYCHOK no cover
384 '''I{Must be overloaded}.'''
385 self._notOverloaded(*C_h_p_z_y_x_r_s_t)
387 def toStr(self, prec=9, **unused): # PYCHOK signature
388 '''Return this C{Ecef*} as a string.
390 @kwarg prec: Precision, number of decimal digits (0..9).
392 @return: This C{Ecef*} (C{str}).
393 '''
394 return self.attrs(_a_, _f_, _datum_, _name_, prec=prec) # _ellipsoid_
397class EcefError(_ValueError):
398 '''An ECEF or C{Ecef*} related issue.
399 '''
400 pass # PYLINT W0107
403class EcefFarrell21(_EcefBase):
404 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)
405 coordinates based on I{Jay A. Farrell}'s U{Table 2.1<https://Books.Google.com/
406 books?id=fW4foWASY6wC>}, page 29, aka the I{Heikkinen application} of U{Ferrari's
407 solution<https://WikiPedia.org/wiki/Geographic_coordinate_conversion>}.
409 @see: Classes L{EcefFarrell22} and L{EcefVeness}.
410 '''
411 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature
412 E = self.ellipsoid
413 a = E.a
414 a2 = E.a2
415 b2 = E.b2
416 e2 = E.e2
417 e2_ = E.e2abs * E.a2_b2 # (E.e * E.a_b)**2 = 0.0820944... WGS84
418 e4 = E.e4
420 z2 = z**2 # names as page 29
421 ez = z2 * (_1_0 - e2) # E.e2s2(z)
423 p2 = p**2
424 G = p2 + ez - e2 * (a2 - b2) # p2 + ez - e4 * a2
425 F = b2 * z2 * 54
426 c = e4 * p2 * F / G**3
427 s = sqrt(c * (c + _2_0))
428 c = cbrt(s + c + _1_0)
429 G *= fsumf_(c, _1_0, _1_0 / c) # k
430 P = F / (G**2 * _3_0)
431 Q = sqrt(_2_0 * e4 * P + _1_0)
432 Q1 = Q + _1_0
433 s = fsumf_(a2 * (Q1 / Q) * _0_5,
434 -P * ez / (Q1 * Q),
435 -P * p2 * _0_5)
436 r = p * P * e2 / Q1 - sqrt(s)
437 r = p + r * e2
438 v = b2 / (sqrt(r**2 + ez) * a) # z0 / z
440 h = hypot(r, z) * (_1_0 - v)
441 z += e2_ * v * z # lat = atan1d(z, p)
442 return z, p, h, None, C
443 # note, phi and lam are swapped on page 29
446class EcefFarrell22(_EcefBase):
447 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)
448 coordinates based on I{Jay A. Farrell}'s U{Table 2.2<https://Books.Google.com/
449 books?id=fW4foWASY6wC>}, page 30, aka I{B.R. Bowring's closed-form}.
451 @see: Classes L{EcefFarrell21} and L{EcefVeness}.
452 '''
453 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature
454 E = self.ellipsoid
455 a, b = E.a, E.b
456 s, c, _ = _norm3(z * a, p * b) # Bowring theta
457 sa, ca = _latB2(E, s, c, z, p) # Bowring lat
458 if ca:
459 h = p / fabs(ca)
460 if sa:
461 h -= E.roc1_(sa)
462 else:
463 h -= a
464# C = 3 # XXX 1?
465 else:
466 h = fabs(z) - b
467 C = 2
468 # lat = atan1d(sa, ca)
469 return sa, ca, h, None, C
470 # note, phi and lam are swapped on page 30
473class EcefFukushima(_EcefBase):
474 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
475 transcoded from I{Toshio Fukushima}'s U{Fortran<https://www.ResearchGate.net/publication/277721539>}
476 implementation.
478 @see: Fukushima, T. U{Transformation from Cartesian to Geodetic Coordinates Accelerated by
479 Halley’s Method<https://www.researchgate.net/publication/227215135>} and Eleiche, M.
480 U{A comparison between Fukushima-Halley algorithm and Trilateration algorithm for
481 geodetic conversion<https://link.Springer.com/article/10.1007/s12145-022-00779-7>}.
482 '''
483 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature
484 E = self.ellipsoid
485 a = E.a
486 e2 = E.e2
487 e4 = _1_5 * E.e4 # e4T
488 ec = _1_0 - E.f # sqrt(_1_0 - e2)
489# assert (a * ec) == E.b
491 za = fabs(z)
492 s0 = za / a
493 zc = s0 * ec
494 pn = p / a
495 # Newton Correction Factors
496 c0 = pn * ec
497 c2 = c0**2
498# c3 = c2 * c0
499 s2 = s0**2
500# s3 = s2 * s0
501# a2 = s2 + c2
502 a0 = hypot(s0, c0) # sqrt(a2)
503 a3 = a0**3 # a0 * a2
504 d0 = a3 * zc + e2 * s2 * s0 # s3
505 f0 = a3 * pn - e2 * c2 * c0 # c3
506 # Halley Correction Factor
507 b0 = e4 * s2 * c2 * pn * (a0 - ec)
508 sa = d0 * f0 - b0 * s0
509 ca = (f0**2 - b0 * c0) * ec
511 # lat = atan1d(sa, ca)
512 h = hypot(ec * sa, ca)
513 h = fsumf_(p * ca, za * sa, -h * a)
514 h = _over(h, hypot(sa, ca))
515 return sa, ca, h, None, C
518class EcefKarney(_EcefBase):
519 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
520 transcoded from I{Karney}'s C++ U{Geocentric<https://GeographicLib.SourceForge.io/C++/doc/
521 classGeographicLib_1_1Geocentric.html>} methods.
523 @note: In general, there are multiple solutions and the result which minimizes C{height} is
524 returned, i.e., the C{(lat, lon)} corresponding to the closest point on the ellipsoid.
525 If there are still multiple solutions with different latitudes (applies only if C{z}
526 = 0), then the solution with C{lat} > 0 is returned. If there are still multiple
527 solutions with different longitudes (applies only if C{x} = C{y} = 0), then C{lon00}
528 is returned. The returned C{lon} is in the range [−180°, 180°] and C{height} is not
529 below M{−E.a * (1 − E.e2) / sqrt(1 − E.e2 * sin(lat)**2)}. Like C{forward} above,
530 M{v1 = Transpose(M) ⋅ v0}.
531 '''
532 def _reverse5(self, C, h, p, z, y, x, r, s, q): # PYCHOK unused y, x
533 E = self.ellipsoid
534 e = E.e4 * q # p renamed to s
535 # Avoid possible division by zero when r = 0 by multiplying
536 # equations for s and t by r^3 and r, respectively.
537 d = s = e * s / _4_0 # s = r^3 * s
538 u = r = r / _6_0
539 r2 = r**2
540 r3 = r2 * r
541 t3 = r3 + s
542 d *= t3 + r3
543 if d < 0:
544 # t is complex, but the way u is defined, the result is real.
545 # There are three possible cube roots. We choose the root
546 # which avoids cancellation. Note, d < 0 implies r < 0.
547 u += cos(atan2(sqrt(-d), -t3) / _3_0) * r * _2_0
548 else:
549 # Pick the sign on the sqrt to maximize abs(t3). This
550 # minimizes loss of precision due to cancellation. The
551 # result is unchanged because of the way the t is used
552 # in definition of u.
553 if d > 0:
554 t3 += copysign0(sqrt(d), t3) # t3 = (r * t)^3
555 # N.B. cbrt always returns the real root, cbrt(-8) = -2.
556 t = cbrt(t3) # t = r * t
557 if t: # t can be zero; but then r2 / t -> 0.
558 u = fsumf_(u, t, r2 / t)
559 v = sqrt(u**2 + e) # guaranteed positive
560 # Avoid loss of accuracy when u < 0. Underflow doesn't occur in
561 # E.e4 * q / (v - u) because u ~ e^4 when q is small and u < 0.
562 u = (e / (v - u)) if u < 0 else (u + v) # u+v, guaranteed positive
563 # Need to guard against w going negative due to roundoff in u - q.
564 w = E.e2abs * (u - q) / (_2_0 * v)
565 # Rearrange expression for k to avoid loss of accuracy due to
566 # subtraction. Division by 0 not possible because u > 0, w >= 0.
567 k1 = k2 = (u / (sqrt(w**2 + u) + w)) if w > 0 else sqrt(u)
568 if E.f < 0:
569 k1 -= E.e2
570 else:
571 k2 += E.e2
572 sa, ca, h = _norm3(z / k1, p / k2)
573 h *= k1 - E.e21
574 return sa, ca, h, None, C
577class EcefOlson(_EcefBase):
578 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
579 based on I{D.K. Olson}'s U{method<https://IEEEXplore.IEEE.org/document/481290>}.
581 @see: U{Possibly wrong<https://PossiblyWrong.Wordpress.com/2014/02/14/>} and U{geodesic-fortran
582 <https://JacobWilliams.GitHub.io/geodesic-fortran/proc/olson.html>}.
583 '''
584 @Property_RO
585 def _E_a1_a4(self):
586 E = self.ellipsoid
587 a1 = E.a * E.e2
588 a2 = a1**2
589 a3 = a1 * E.e2 * _0_5
590 a4 = a2 * 2.5
591 return E, a1, a2, a3, a4
593 def _reverse5(self, C, r, p, z, *unused): # PYCHOK signature
594 E, a1, a2, a3, a4 = self._E_a1_a4
595 _1_r = _1_0 / r
597 u = a2 * _1_r
598 v = a3 - a4 * _1_r
599 za = fabs(z)
600 try:
601 s = za * _1_r
602 c = p * _1_r
603 if c > 0.5477225575: # == sqrt(0.3)
604 c *= (a1 + u + s**2 * v) * c * _1_r
605 s *= _1_0 + c
606 ra = asin(s)
607 s2 = s**2
608 c = sqrt(_1_0 - s2)
609 else:
610 u -= a3
611 s *= (a1 - u - c**2 * v) * s * _1_r
612 c *= _1_0 - s
613 ra = acos(c)
614 s2 = _1_0 - c**2
615 s = sqrt(s2)
616 g = _1_0 - E.e2 * s2
617 rg = E.a / sqrt(g)
618 rf = rg * E.e21
619 u = p - rg * c
620 v = za - rf * s
622 h = c * u + s * v
623 m = c * v - s * u
624 r = m * g / (h * g + rf)
625 h += m * r * _0_5
626 ra += r
627 if z < 0:
628 ra = -ra
629 sa, ca = sincos2(ra)
631 except (TypeError, ValueError): # domain error
632 sa, ca, h, C = _subterrestrial4(E, z, p)
633 return sa, ca, h, None, C
636class EcefOsen(_EcefBase):
637 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
638 based on I{Karl Osen}'s U{boost<https://HAL.Science/hal-01704943v2>} of I{Jijie Zhu}'s U{algorithm
639 <https://IEEEXplore.IEEE.org/document/303772>}.
640 '''
641 @Property_RO
642 def _E_Hmin_L_L2_E21a2(self):
643 E = self.ellipsoid
644 Hmin = pow(E.e, 12) / _6_0 # Hmin 64-bit
645 L = E.e2 * _0_5
646 L2 = E.e4 * _0_25 # L**2
647 E21a2 = E.e21 / E.a2
648 return E, Hmin, L, L2, E21a2
650 def _reverse5(self, C, _, w, z, *unused): # PYCHOK signature
651 E, Hmin, L, L2, E21a2 = self._E_Hmin_L_L2_E21a2
652 w2 = w**2
653 m = w2 / E.a2 # invaa
654 n = z**2 * E21a2 # p1meedaa
655 mn = m + n # mpn
656 P = (mn - E.e4) / _6_0 # L2 * _4_0
657 G = m * n * L2
658 H = P**3 * _2_0 + G
659 if H < Hmin: # below 43 Km from center
660 sa, ca, h, C = _subterrestrial4(E, z, w)
661 else:
662 i = -(mn * _0_5 + L2)
663 B = cbrt((H + G) * _0_5 + sqrt(H * G)) # pow(H + G + 2 * sqrt(H * G), inv3) * invcbrt2
664 B = i / _3_0 - B - P**2 / B
665 k = (L2 - mn) * L2
667 # compute left part of t
668 t2 = sqrt(B**2 - k) # t1
669 t3 = t2 - (B + i) * _0_5
670 # compute right part of t
671 t5 = (B - i) * _0_5
672 # t5 may accidentally drop just below zero due to
673 # numeric turbulence at lats near +/- 45.3 degrees
674 t6 = sqrt(fabs(t5))
675 # add left and right parts
676 t = sqrt(t3) + (t6 if m < n else -t6)
678 # Newton-Raphson to compute t correction
679 g = (m - n) * E.e2 # L * 2
680 F = ((t**2 + i * _2_0) * t + g) * t + k
681 D = (t**3 + i * t) * _4_0 + g # dFdt
682 dt = -F / D
684 # compute latitude (range -PI/2..PI/2)
685 t += dt
686 u = t + L
687 v = t - L
688 sa = z * u # zu
689 ca = w * v # wv
691 uv = u * v # invuv
692 h = hypot(w - ca / uv,
693 z - sa / uv * E.e21) # p1mee
694 if u < 1:
695 h = -h
696 return sa, ca, h, None, C
699class EcefSudano(_EcefBase):
700 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
701 based on I{John J. Sudano}'s U{paper<https://www.ResearchGate.net/publication/3709199>}.
702 '''
703 _TOL = \
704 _tol = EPS
706 def reverse(self, xyz, y=None, z=None, M=False, tol=EPS, **lon00_name): # PYCHOK tol
707 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using
708 I{Sudano}'s U{iterative method<https://www.ResearchGate.net/publication/3709199>}.
710 @kwarg tol: Convergence tolerance for C{sin(latitude)} (C{scalar}).
712 @see: L{Parent method<_EcefBase.reverse>} for all other information.
714 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}} not
715 C{scalar} for C{scalar} B{C{xyz}} or no convergence for C{B{tol}}.
716 '''
717 if tol != self._TOL:
718 self._tol = Scalar_(tol=tol, low=EPS, Error=EcefError)
719 return _EcefBase.reverse(self, xyz, y=y, z=z, M=M, **lon00_name)
721 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature
722 E = self.ellipsoid
723 e = E.e2 * E.a
724 d = e - p
726 sa, ca, _ = _norm3(fabs(z), p * E.e21)
727 # Sudano's Eq (A-6) and (A-7) refactored/reduced,
728 # replacing Rn from Eq (A-4) with n = E.a / ca:
729 # N = ca**2 * ((z + E.e2 * n * sa) * ca - p * sa)
730 # = ca**2 * (z * ca + E.e2 * E.a * sa - p * sa)
731 # = ca**2 * (z * ca + (E.e2 * E.a - p) * sa)
732 # D = ca**3 * (E.e2 * n / E.e2s2(sa)) - p
733 # = ca**2 * (E.e2 * E.a / E.e2s2(sa) - p / ca**2)
734 # N / D = (z * ca + (E.e2 * E.a - p) * sa) /
735 # (E.e2 * E.a / E.e2s2(sa) - p / ca**2)
736 tol = self._tol
737 _S2 = Fsum(sa).fsum2f_
738 for i in range(1, _TRIPS): # 6+ max
739 ca2 = _1_0 - sa**2
740 if ca2 < EPS02:
741 break
742 D = p / ca2 - e / E.e2s2(sa)
743 if fabs(D) < EPS0:
744 break
745 ca = sqrt(ca2)
746 sa, D = _S2(z * ca / D, d * sa / D)
747 if fabs(D) < tol:
748 break
749 else: # PYCHOK no cover
750 raise ValueError(Fmt.no_convergence(fabs(D), tol))
752 sa = copysign0(sa, z)
753 # lat = atan1d(sa, ca)
754 # h = (fabs(z) + p - E.a * cos(a + E.e21) * sa / ca) / (ca + sa)
755 # Sudano's Eq (7) doesn't produce the correct height, ...
756 h = _heightB(E, sa, ca, z, p) # ... use Veness' (Bowring eqn 7)
757 return sa, ca, h, i, C
759 @deprecated_property
760 def tolerance(self):
761 '''DEPRECATED on 2025.08.22, use keyword argument C{tol}.'''
762 return self._tol
764 @tolerance.setter # PYCHOK setter!
765 def tolerance(self, tol):
766 self._tol = Scalar_(tolerance=tol, low=EPS, Error=EcefError)
769class EcefUPC(_EcefBase):
770 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates based on
771 I{UPC}'s U{method<https://GSSC.ESA.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion>}.
772 '''
773 _TOL = _TOL
774 _tol = _Degrees2Radians(_TOL)
776 def reverse(self, xyz, y=None, z=None, M=False, tol=_TOL, **lon00_name): # PYCHOK tol
777 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using I{UPC}'s
778 U{iterative method<https://GSSC.ESA.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf>}, page 186.
780 @kwarg tol: Convergence tolerance for the C{latitude} (C{degrees}).
782 @see: L{Parent method<_EcefBase.reverse>} for all other information.
784 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}} not
785 C{scalar} for C{scalar} B{C{xyz}} or no convergence for C{B{tol}}.
786 '''
787 if tol != _TOL:
788 self._tol = _Degrees2Radians(tol)
789 return _EcefBase.reverse(self, xyz, y=y, z=z, M=M, **lon00_name)
791 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature
792 E = self.ellipsoid
793 a = E.a
794 e2 = E.e2 # signed
796 za = fabs(z)
797 ph_ = atan1(za, E.e21 * p)
798 tol = self._tol
799 for i in range(1, _TRIPS): # 5..6 max
800 s, c = sincos2(ph_)
801 N = a / sqrt(_1_0 - s**2 * e2) # N + h == N + p / c - N == p / c
802 ca = p - N * c * e2 # == p * (1 - N * e2 / (N + h)) == p * (1 - N * e2 * c / p)
803 ph = atan1(za, ca) # atan1(z / p, 1 - N * e2 / (N + h)) == atan1(z, ca)
804 r = fabs(ph - ph_)
805 if r < tol:
806 # lat = copysign0(degrees(ph), z)
807 # == atan1d(z, ca)
808 h = p / c - N
809 break
810 ph_ = ph
811 else: # PYCHOK no cover
812 r, tol = map1(degrees, r, tol)
813 raise ValueError(Fmt.no_convergence(r, tol))
814 return z, ca, h, i, C
817class EcefVeness(_EcefBase):
818 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
819 transcoded from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal, Cartesian<https://
820 www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}.
822 @note: Uses B. R. Bowring’s formulation for μm precision in concise form U{The accuracy of
823 geodetic latitude and height equations<https://www.ResearchGate.net/publication/233668213>},
824 Survey Review, Vol 28, 218, Oct 1985.
826 @see: U{A Guide to Coordinate Systems in Great Britain<https://www.OrdnanceSurvey.co.UK/documents/
827 resources/guide-coordinate-systems-great-britain.pdf>}, section I{B) Converting between 3D
828 Cartesian and ellipsoidal latitude, longitude and height coordinates}.
830 @see: Toms, Ralph M. U{An Efficient Algorithm for Geocentric to Geodetic Coordinate Conversion
831 <https://www.OSTI.gov/scitech/biblio/110235>}, Sept 1995 and U{An Improved Algorithm for
832 Geocentric to Geodetic Coordinate Conversion<https://www.OSTI.gov/scitech/servlets/purl/231228>},
833 Apr 1996, both from Lawrence Livermore National Laboratory (LLNL).
834 '''
835 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature
836 # assert h >= p > 0 # h = hypot(z, p)
837 E = self.ellipsoid
838 B = E.b * E.e22
839 # parametric latitude, Bowring eqn 17, replaced
840 t = (E.b * z) / (E.a * p) * (B / h + _1_0) # theta
841 c = _1_0 / hypot1(t) # t == atan2(z * E.a, p * E.b)
842 s = c * t # s, c == sincos2(t)
843 sa, ca = _latB2(E, s, c, z, p) # lat, Bowring eqn 18
844 # lat = atan1d(sa, ca)
845 h = _heightB(E, sa, ca, z, p) # height, Bowring eqn 7
846 return sa, ca, h, None, C
849class EcefYou(_EcefBase):
850 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
851 using I{Rey-Jer You}'s U{transformation<https://www.ResearchGate.net/publication/240359424>}
852 for I{non-prolate} ellipsoids.
854 @see: Featherstone, W.E., Claessens, S.J. U{Closed-form transformation between geodetic and
855 ellipsoidal coordinates<https://Espace.Curtin.edu.AU/bitstream/handle/20.500.11937/11589/
856 115114_9021_geod2ellip_final.pdf>} Studia Geophysica et Geodaetica, 2008, 52, pages 1-18
857 and U{PyMap3D<https://PyPI.org/project/pymap3d>}.
858 '''
859 _e_e2 = True # or 2-tuple C{(e, e2)}
861 def _reverse5(self, C, h, p, z, y, x, *unused): # PYCHOK signature
862 E = self.ellipsoid
863 a, b = E.a, E.b
864 e, e2 = self._e_e2
866 u = hypot2_(x, y, z) - e2
867 u += hypot(u, e * z * _2_0)
868 u *= _0_5
869 if u > EPS02:
870 u = sqrt(u)
871 q = hypot(u, e)
872 B = atan1(q * z, u * p) # beta0 = atan(q / u * z / p)
873 sB, cB = sincos2(B)
874 if cB and sB:
875 q *= a
876 d = (q / cB - e2 * cB) / sB
877 if isnon0(d):
878 B += fsumf_(u * b, -q, e2) / d
879 sB, cB = sincos2(B)
880 elif u < (-EPS02):
881 raise EcefError(u=u, txt=_singular_)
882 else: # near polar # PYCHOK no cover
883 sB, cB, C = _copysign_1_0(z), _0_0, 2
885 h = hypot(p - a * cB, z - b * sB)
886 if hypot2_(x, y, z * E.a_b) < E.a2: # or lat < 0 or z < 0
887 h = neg(h) # inside ellipsoid
888 # lat = atand(E.a_b * tan(B)) == atan1d(a * sB, b * cB)
889 return (a * sB), (b * cB), h, None, C
892class EcefMatrix(_NamedTuple):
893 '''A rotation matrix known as I{East-North-Up (ENU) to ECEF}.
895 @see: U{From ENU to ECEF<https://WikiPedia.org/wiki/
896 Geographic_coordinate_conversion#From_ECEF_to_ENU>} and
897 U{Issue #74<https://Github.com/mrJean1/PyGeodesy/issues/74>}.
898 '''
899 _Names_ = ('_0_0_', '_0_1_', '_0_2_', # row-major order
900 '_1_0_', '_1_1_', '_1_2_',
901 '_2_0_', '_2_1_', '_2_2_')
902 _Units_ = (Scalar,) * len(_Names_)
904 def _validate(self, **unused): # PYCHOK unused
905 '''(INTERNAL) Allow C{_Names_} with leading underscore.
906 '''
907 _NamedTuple._validate(self, underOK=True)
909 def __new__(cls, sa, ca, sb, cb, *_more, **name):
910 '''New L{EcefMatrix} matrix.
912 @arg sa: C{sin(phi)} (C{float}).
913 @arg ca: C{cos(phi)} (C{float}).
914 @arg sb: C{sin(lambda)} (C{float}).
915 @arg cb: C{cos(lambda)} (C{float}).
916 @arg _more: (INTERNAL) from C{.multiply}.
918 @raise EcefError: If B{C{sa}}, B{C{ca}}, B{C{sb}} or
919 B{C{cb}} outside M{[-1.0, +1.0]}.
920 '''
921 t = sa, ca, sb, cb
922 if _more: # all 9 matrix elements ...
923 t += _more # ... from .multiply
925 elif max(map(fabs, t)) > _1_0:
926 raise EcefError(unstr(EcefMatrix, *t))
928 else: # build matrix from the following quaternion operations
929 # qrot(lam, [0,0,1]) * qrot(phi, [0,-1,0]) * [1,1,1,1]/2
930 # or
931 # qrot(pi/2 + lam, [0,0,1]) * qrot(-pi/2 + phi, [-1,0,0])
932 # where
933 # qrot(t,v) = [cos(t/2), sin(t/2)*v[1], sin(t/2)*v[2], sin(t/2)*v[3]]
935 # Local X axis (East) in geocentric coords
936 # M[0] = -slam; M[3] = clam; M[6] = 0;
937 # Local Y axis (North) in geocentric coords
938 # M[1] = -clam * sphi; M[4] = -slam * sphi; M[7] = cphi;
939 # Local Z axis (Up) in geocentric coords
940 # M[2] = clam * cphi; M[5] = slam * cphi; M[8] = sphi;
941 t = (-sb, -cb * sa, cb * ca,
942 cb, -sb * sa, sb * ca,
943 _0_0, ca, sa)
945 return _NamedTuple.__new__(cls, *t, **name)
947 def column(self, column):
948 '''Get this matrix' B{C{column}} 0, 1 or 2 as C{3-tuple}.
949 '''
950 if 0 <= column < 3:
951 return self[column::3]
952 raise _IndexError(column=column)
954 @property_RO
955 def _columns(self):
956 for c in range(3):
957 yield self[c::3]
959 def copy(self, **unused): # PYCHOK signature
960 '''Make a shallow or deep copy of this instance.
962 @return: The copy (C{This class} or subclass thereof).
963 '''
964 return self.classof(*self)
966 __copy__ = __deepcopy__ = copy
968 @Property_RO
969 def matrix3(self):
970 '''Get this matrix' rows (C{3-tuple} of 3 C{3-tuple}s).
971 '''
972 return tuple(self._rows)
974 @Property_RO
975 def matrixTransposed3(self):
976 '''Get this matrix' I{Transposed} rows (C{3-tuple} of 3 C{3-tuple}s).
977 '''
978 return tuple(self._columns)
980 def multiply(self, other):
981 '''Matrix multiply M{M0' ⋅ M} this matrix I{Transposed} with an other matrix.
983 @arg other: The other matrix (L{EcefMatrix}).
985 @return: The matrix product (L{EcefMatrix}).
987 @raise TypeError: If B{C{other}} is not an L{EcefMatrix}.
988 '''
989 _xinstanceof(EcefMatrix, other=other)
990 # like LocalCartesian.MatrixMultiply, C{self.matrixTransposed3 X other.matrix3}
991 # <https://GeographicLib.SourceForge.io/C++/doc/LocalCartesian_8cpp_source.html>
992 X = (_fdotf(t, *c) for t in self._columns for c in other._columns)
993 return _xnamed(EcefMatrix(*X), typename(EcefMatrix.multiply))
995 def rotate(self, xyz, *xyz0):
996 '''Forward rotation M{M0' ⋅ ([x, y, z] - [x0, y0, z0])'}.
998 @arg xyz: Local C{(x, y, z)} coordinates (C{3-tuple}).
999 @arg xyz0: Optional, local C{(x0, y0, z0)} origin (C{3-tuple}).
1001 @return: Rotated C{(x, y, z)} location (C{3-tuple}).
1003 @raise LenError: Unequal C{len(B{xyz})} and C{len(B{xyz0})}.
1004 '''
1005 if xyz0:
1006 if len(xyz0) != len(xyz):
1007 raise LenError(self.rotate, xyz0=len(xyz0), xyz=len(xyz))
1008 xyz = tuple(s - s0 for s, s0 in zip(xyz, xyz0))
1010 # x' = M[0] * x + M[3] * y + M[6] * z
1011 # y' = M[1] * x + M[4] * y + M[7] * z
1012 # z' = M[2] * x + M[5] * y + M[8] * z
1013 return tuple(_fdotf(xyz, *c) for c in self._columns)
1015 def row(self, row):
1016 '''Get this matrix' B{C{row}} 0, 1 or 2 as C{3-tuple}.
1017 '''
1018 if 0 <= row < 3:
1019 r = row * 3
1020 return self[r:r+3]
1021 raise _IndexError(row=row)
1023 @property_RO
1024 def _rows(self):
1025 for r in (0, 3, 6):
1026 yield self[r:r+3]
1028 def unrotate(self, xyz, *xyz0):
1029 '''Inverse rotation M{[x0, y0, z0] + M0 ⋅ [x,y,z]'}.
1031 @arg xyz: Local C{(x, y, z)} coordinates (C{3-tuple}).
1032 @arg xyz0: Optional, local C{(x0, y0, z0)} origin (C{3-tuple}).
1034 @return: Unrotated C{(x, y, z)} location (C{3-tuple}).
1036 @raise LenError: Unequal C{len(B{xyz})} and C{len(B{xyz0})}.
1037 '''
1038 if xyz0:
1039 if len(xyz0) != len(xyz):
1040 raise LenError(self.unrotate, xyz0=len(xyz0), xyz=len(xyz))
1041 _xyz = _1_0_1T + xyz
1042 # x' = x0 + M[0] * x + M[1] * y + M[2] * z
1043 # y' = y0 + M[3] * x + M[4] * y + M[5] * z
1044 # z' = z0 + M[6] * x + M[7] * y + M[8] * z
1045 xyz_ = (_fdotf(_xyz, s, *r) for s, r in zip(xyz0, self._rows))
1046 else:
1047 # x' = M[0] * x + M[1] * y + M[2] * z
1048 # y' = M[3] * x + M[4] * y + M[5] * z
1049 # z' = M[6] * x + M[7] * y + M[8] * z
1050 xyz_ = (_fdotf(xyz, *r) for r in self._rows)
1051 return tuple(xyz_)
1054class Ecef9Tuple(_NamedTuple, _EcefLocal):
1055 '''9-Tuple C{(x, y, z, lat, lon, height, C, M, datum)} with I{geocentric} C{x},
1056 C{y} and C{z} plus I{geodetic} C{lat}, C{lon} and C{height}, case C{C} and
1057 optionally, rotation matrix C{M} (L{EcefMatrix} or C{None}) and C{datum},
1058 with C{lat} and C{lon} in C{degrees} and C{x}, C{y}, C{z} and C{height} in
1059 C{meter}, conventionally. Case C{C=0} means C{x, y,z} from foward, C{C=1}
1060 C{lat, lon, height} from reverse, C{C=2} near-polar C{lat, lon}, C{C=3}
1061 near-equatorial C{lat, lon}, C{C=4} spherical C{lat, lon}, C{C=5} if the
1062 C{height} exceeds C{datum}'s C{ellipsoid.heightMax} and C{C=6} for points
1063 too subterrestrial.
1064 '''
1065 _Names_ = (_x_, _y_, _z_, _lat_, _lon_, _height_, _C_, _M_, _datum_)
1066 _Units_ = ( Meter, Meter, Meter, Lat, Lon, Height, Int, _Pass, _Pass)
1068 @property_ROver
1069 def _CartesianBase(self):
1070 '''(INTERNAL) Get class C{CartesianBase}, I{once}.
1071 '''
1072 return _MODS.cartesianBase.CartesianBase # overwrite property_ROver
1074 @deprecated_method
1075 def convertDatum(self, datum2): # for backward compatibility
1076 '''DEPRECATED, use method L{toDatum}.'''
1077 return self.toDatum(datum2)
1079 @property_RO
1080 def _ecef9(self): # in ._EcefLocal._Ltp_ecef2local
1081 return self
1083 @property_RO
1084 def ellipsoid(self):
1085 '''Get the ellipsoid (L{Ellipsoid}).
1086 '''
1087 return (self.datum or _WGS84).ellipsoid
1089 @Property_RO
1090 def lam(self):
1091 '''Get the longitude in C{radians} (C{float}).
1092 '''
1093 return self.philam.lam
1095 @Property_RO
1096 def lamVermeille(self):
1097 '''Get the longitude in C{radians} M{[-PI*3/2..+PI*3/2]} after U{Vermeille
1098 <https://Search.ProQuest.com/docview/639493848>} (2004), page 95.
1100 @see: U{Karney<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>},
1101 U{Vermeille<https://Search.ProQuest.com/docview/847292978>} 2011, pp 112-113, 116
1102 and U{Featherstone, et.al.<https://Search.ProQuest.com/docview/872827242>}, page 7.
1103 '''
1104 x, y = self.x, self.y
1105 a = fabs(y)
1106 if a > EPS0:
1107 r = PI_2 - atan2(x, hypot(x, a) + a) * _2_0
1108 if y < 0:
1109 r = -r
1110 else: # y == 0
1111 r = PI if x < 0 else _0_0
1112 return Lam(Vermeille=r)
1114 @Property_RO
1115 def latlon(self):
1116 '''Get the lat-, longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).
1117 '''
1118 return LatLon2Tuple(self.lat, self.lon, name=self.name)
1120 @Property_RO
1121 def latlonheight(self):
1122 '''Get the lat-, longitude in C{degrees} and height (L{LatLon3Tuple}C{(lat, lon, height)}).
1123 '''
1124 return self.latlon.to3Tuple(self.height)
1126 @Property_RO
1127 def latlonheightdatum(self):
1128 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}).
1129 '''
1130 return self.latlonheight.to4Tuple(self.datum)
1132 @Property_RO
1133 def latlonVermeille(self):
1134 '''Get the latitude and I{Vermeille} longitude in C{degrees [-225..+225]} (L{LatLon2Tuple}C{(lat, lon)}).
1136 @see: Property C{lonVermeille}.
1137 '''
1138 v = self.lonVermeille
1139 return LatLon2Tuple(self.lat, v, name=self.name + v.name)
1141 @Property_RO
1142 def lonVermeille(self):
1143 '''Get the longitude in C{degrees [-225..+225]} after U{Vermeille
1144 <https://Search.ProQuest.com/docview/639493848>} 2004, page 95.
1146 @see: Property C{lamVermeille}.
1147 '''
1148 return Lon(Vermeille=degrees(self.lamVermeille))
1150 @Property_RO
1151 def Mx(self):
1152 '''Compute rotation matrix (L{EcefMatrix}), seperate from C{M}.
1153 '''
1154 sa, ca, sb, cb, _, _, _ = _norm7(self.y, self.x, self.z, self.ellipsoid)
1155 return EcefMatrix(sa, ca, sb, cb, name=self.name)
1157 @Property_RO
1158 def phi(self):
1159 '''Get the latitude in C{radians} (C{float}).
1160 '''
1161 return self.philam.phi
1163 @Property_RO
1164 def philam(self):
1165 '''Get the lat-, longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).
1166 '''
1167 return PhiLam2Tuple(radians(self.lat), radians(self.lon), name=self.name)
1169 @Property_RO
1170 def philamheight(self):
1171 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
1172 '''
1173 return self.philam.to3Tuple(self.height)
1175 @Property_RO
1176 def philamheightdatum(self):
1177 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLam4Tuple}C{(phi, lam, height, datum)}).
1178 '''
1179 return self.philamheight.to4Tuple(self.datum)
1181 @Property_RO
1182 def philamVermeille(self):
1183 '''Get the latitude and I{Vermeille} longitude in C{radians [-PI*3/2..+PI*3/2]} (L{PhiLam2Tuple}C{(phi, lam)}).
1185 @see: Property C{lamVermeille}.
1186 '''
1187 v = self.lamVermeille
1188 return PhiLam2Tuple(radians(self.lat), v, name=self.name + v.name)
1190 phiVermeille = phi
1192 def toCartesian(self, Cartesian=None, **Cartesian_kwds):
1193 '''Return the geocentric C{(x, y, z)} coordinates as an ellipsoidal or spherical
1194 C{Cartesian}.
1196 @kwarg Cartesian: Optional class to return C{(x, y, z)} (L{ellipsoidalKarney.Cartesian},
1197 L{ellipsoidalNvector.Cartesian}, L{ellipsoidalVincenty.Cartesian},
1198 L{sphericalNvector.Cartesian} or L{sphericalTrigonometry.Cartesian})
1199 or C{None}.
1200 @kwarg Cartesian_kwds: Optionally, additional B{C{Cartesian}} keyword arguments, ignored
1201 if C{B{Cartesian} is None}.
1203 @return: A B{C{Cartesian}} instance or a L{Vector4Tuple}C{(x, y, z, h)} if C{B{Cartesian}
1204 is None}.
1206 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}} item.
1207 '''
1208 if _isin(Cartesian, None, Vector4Tuple):
1209 r = self.xyzh
1210 elif Cartesian is Vector3Tuple:
1211 r = self.xyz
1212 else:
1213 _xsubclassof(self._CartesianBase, Cartesian=Cartesian)
1214 r = Cartesian(self, **_name1__(Cartesian_kwds, _or_nameof=self))
1215 return r
1217 def toDatum(self, datum2, **name):
1218 '''Convert this C{Ecef9Tuple} to an other datum.
1220 @arg datum2: Datum to convert I{to} (L{Datum}).
1221 @kwarg name: Optional C{B{name}=NN} (C{str}).
1223 @return: The converted 9-Tuple (C{Ecef9Tuple}).
1225 @raise TypeError: The B{C{datum2}} is not a L{Datum}.
1226 '''
1227 n = _name__(name, _or_nameof=self)
1228 if _isin(self.datum, None, datum2): # PYCHOK _Names_
1229 r = self.copy(name=n)
1230 else:
1231 c = self._CartesianBase(self, datum=self.datum, name=n) # PYCHOK _Names_
1232 # c.toLatLon converts datum, x, y, z, lat, lon, etc.
1233 # and returns another Ecef9Tuple iff LatLon is None
1234 r = c.toLatLon(datum=datum2, LatLon=None)
1235 return r
1237 def toLatLon(self, LatLon=None, **LatLon_kwds):
1238 '''Return the geodetic C{(lat, lon, height[, datum])} coordinates.
1240 @kwarg LatLon: Optional class to return C{(lat, lon, height[, datum])} or C{None}.
1241 @kwarg LatLon_kwds: Optional B{C{height}}, B{C{datum}} and other B{C{LatLon}}
1242 keyword arguments.
1244 @return: A B{C{LatLon}} instance or if C{B{LatLon} is None}, a L{LatLon4Tuple}C{(lat,
1245 lon, height, datum)} or L{LatLon3Tuple}C{(lat, lon, height)} if C{datum} is
1246 specified or not.
1248 @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}} item.
1249 '''
1250 lat, lon, D = self.lat, self.lon, self.datum # PYCHOK Ecef9Tuple
1251 kwds = _name1__(LatLon_kwds, _or_nameof=self)
1252 kwds = _xkwds(kwds, height=self.height, datum=D) # PYCHOK Ecef9Tuple
1253 d = kwds.get(_datum_, LatLon)
1254 if LatLon is None:
1255 r = LatLon3Tuple(lat, lon, kwds[_height_], name=kwds[_name_])
1256 if d is not None:
1257 # assert d is not LatLon
1258 r = r.to4Tuple(d) # checks type(d)
1259 else:
1260 if d is None:
1261 _ = kwds.pop(_datum_) # remove None datum
1262 r = LatLon(lat, lon, **kwds)
1263 _xdatum(_xattr(r, datum=D), D)
1264 return r
1266 def toVector(self, Vector=None, **Vector_kwds):
1267 '''Return these geocentric C{(x, y, z)} coordinates as vector.
1269 @kwarg Vector: Optional vector class to return C{(x, y, z)} or C{None}.
1270 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword arguments,
1271 ignored if C{B{Vector} is None}.
1273 @return: A B{C{Vector}} instance or a L{Vector3Tuple}C{(x, y, z)} if
1274 C{B{Vector} is None}.
1276 @raise TypeError: Invalid B{C{Vector}} or B{C{Vector_kwds}} item.
1278 @see: Propertes C{xyz} and C{xyzh}
1279 '''
1280 return self.xyz if Vector is None else Vector(
1281 *self.xyz, **_name1__(Vector_kwds, _or_nameof=self)) # PYCHOK Ecef9Tuple
1283# def _T_x_M(self, T):
1284# '''(INTERNAL) Update M{self.M = T.multiply(self.M)}.
1285# '''
1286# return self.dup(M=T.multiply(self.M))
1288 @Property_RO
1289 def xyz(self):
1290 '''Get the geocentric C{(x, y, z)} coordinates (L{Vector3Tuple}C{(x, y, z)}).
1291 '''
1292 return Vector3Tuple(self.x, self.y, self.z, name=self.name)
1294 @Property_RO
1295 def xyzh(self):
1296 '''Get the geocentric C{(x, y, z)} coordinates and C{height} (L{Vector4Tuple}C{(x, y, z, h)})
1297 '''
1298 return self.xyz.to4Tuple(self.height)
1301def _4Ecef(this, Ecef): # in .datums.Datum.ecef, .ellipsoids.Ellipsoid.ecef
1302 '''Return an ECEF converter for C{this} L{Datum} or L{Ellipsoid}.
1303 '''
1304 if Ecef is None:
1305 Ecef = EcefKarney
1306 else:
1307 _xinstanceof(*_Ecefs, Ecef=Ecef)
1308 return Ecef(this, name=this.name)
1311def _equatorial2(E, p, z):
1312 '''(INTERNAL) Equatorial plane from C{EcefKarney}.
1313 '''
1314 # Treat prolate spheroids by swapping p and z here and by
1315 # switching the arguments to phi = atan2(...) at the end
1316 # of method C{EcefKarney._reverse5}
1317 p = (p / E.a)**2
1318 z = (z / E.a)**2 * E.e21
1319 return (z, p) if E.f < 0 else (p, z)
1322def _equatorial4(E, s, z):
1323 '''(INTERNAL) C{C=3}, near equatorial plane from C{EcefKarney}.
1324 '''
1325 t = E.e4 - s
1326 if E.f < 0:
1327 s, t = t, s
1328 e = E.a
1329 else:
1330 e = E.b2_a
1331 sa, ca, h = _norm3(*map1(sqrt, E._1_e21 * t, s))
1332 if z < 0: # for tiny negative z, not for prolate
1333 sa = neg(sa)
1334 h *= neg(e / E.e2abs)
1335 return sa, ca, h, 3 # equatorial
1338def _heightB(E, sa, ca, z, p): # in ecef.EcefSudano, ecec.EcefVeness
1339 '''(INTERNAL) Height above ellipsoid (Bowring eqn 7) at C{lat}.
1340 '''
1341 # sa, ca = sincos2d(lat)
1342 # p = hypot(x, y) # distance to polar axis
1344 # r = a / E.e2s(sa) # length of normal terminated by polar axis
1345 # h = p * ca + z * sa - (a * a / r)
1346 return (p * ca + fabs(z * sa) - E.a * E.e2s(sa)) if sa else (p - E.a)
1349def _latB2(E, s, c, z, p):
1350 '''(INTERNAL) Geodetic latitude, Bowring eqn 18.
1351 '''
1352 sa, ca, _ = _norm3(z + s**3 * E.b * E.e22,
1353 p - c**3 * E.a * E.e2)
1354 return sa, ca
1357def _llhn4(latlonh, lon, height, suffix=NN, Error=EcefError, **name): # in .ltp
1358 '''(INTERNAL) Get a C{(lat, lon, h, name)} 4-tuple.
1359 '''
1360 try:
1361 lat, lon = latlonh.lat, latlonh.lon
1362 h = _xattr(latlonh, height=_xattr(latlonh, h=height))
1363 n = _name__(name, _or_nameof=latlonh) # == latlonh._name__(name)
1364 except AttributeError:
1365 lat, h, n = latlonh, height, _name__(**name)
1366 try:
1367 return Lat(lat), Lon(lon), Height(h), n
1368 except (TypeError, ValueError) as x:
1369 t = _lat_, _lon_, _height_
1370 if suffix:
1371 t = (_ + suffix for _ in t)
1372 d = dict(zip(t, (lat, lon, h)))
1373 raise Error(cause=x, **d)
1376def _Lon00(lon00):
1377 '''(INTERNAL) Return C{lon00}.
1378 '''
1379 return Lon(lon00=lon00, clip=0)
1382def _norm3(y, x, eps=0):
1383 '''(INTERNAL) Return C{y, x, h} normalized.
1384 '''
1385 h = hypot(y, x) # EPS0, EPS_2
1386 return (y / h, x / h, h) if h > eps else (_0_0, _1_0, h) # copysign_1_0(x)
1389def _norm7(y, x, z=0, E=_EWGS84):
1390 '''(INTERNAL) Return C{phi, lam, h, p, C}.
1391 '''
1392 sb, cb, p = _norm3(y, x) # lam, distance to polar axis
1393 sa, ca, h = _norm3(z, p) # phi, distance to earth center
1394 if h > E.heightMax:
1395 # We are really far away (> 12M light years). Treat the earth
1396 # as a point and h above as an acceptable approximation to the
1397 # height. This avoids overflow, e.g., in the computation of d
1398 # below. It's possible that h has overflowed to INF, that's OK.
1399 # Treat finite x, y, but p overflows to +INF by scaling by 2.
1400 sb, cb, p = _norm3(y * _0_5, x * _0_5)
1401 sa, ca, _ = _norm3(z * _0_5, p)
1402 C = 5 # near polar
1403 else:
1404 C = None # tbd
1405 return sa, ca, sb, cb, h, p, C
1408def _subterrestrial4(E, z, p):
1409 '''(INTERNAL) Case C{C=6}, too subterrestrial.
1410 '''
1411 sa, ca, h = _norm3(z, p) # like _norm7
1412 h -= E.b if fabs(z) > p else E.a
1413 return sa, ca, h, 6 # near center
1416def _xEcef(Ecef): # PYCHOK .latlonBase
1417 '''(INTERNAL) Validate B{C{Ecef}} I{class}.
1418 '''
1419 if issubclassof(Ecef, _EcefBase):
1420 return Ecef
1421 raise _TypesError(_Ecef_, Ecef, *_Ecefs)
1424# kwd lon00 unused but will throw a TypeError if misspelled, etc.
1425def _xyzn4(xyz, y, z, Types, Error=EcefError, lon00=0, # PYCHOK unused
1426 _xyz_y_z_names=_xyz_y_z, **name): # in .ltp, .pychlv
1427 '''(INTERNAL) Get an C{(x, y, z, name)} 4-tuple.
1428 '''
1429 try:
1430 n = _name__(name, _or_nameof=xyz) # == xyz._name__(name)
1431 try:
1432 t = xyz.x, xyz.y, xyz.z, n
1433 if not isinstance(xyz, Types):
1434 raise _TypesError(_xyz_y_z_names[0], xyz, *Types)
1435 except AttributeError:
1436 t = map1(float, xyz, y, z) + (n,)
1437 except (TypeError, ValueError) as x:
1438 d = dict(zip(_xyz_y_z_names, (xyz, y, z)))
1439 raise Error(cause=x, **d)
1440 return t
1441# assert _xyz_y_z == _args_kwds_names(_xyzn4)[:3]
1444_Ecefs = tuple(_ for _ in locals().values()
1445 if issubclassof(_, _EcefBase) and
1446 _ is not _EcefBase)
1447__all__ += _ALL_DOCS(_EcefBase)
1449# **) MIT License
1450#
1451# Copyright (C) 2016-2027 -- mrJean1 at Gmail -- All Rights Reserved.
1452#
1453# Permission is hereby granted, free of charge, to any person obtaining a
1454# copy of this software and associated documentation files (the "Software"),
1455# to deal in the Software without restriction, including without limitation
1456# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1457# and/or sell copies of the Software, and to permit persons to whom the
1458# Software is furnished to do so, subject to the following conditions:
1459#
1460# The above copyright notice and this permission notice shall be included
1461# in all copies or substantial portions of the Software.
1462#
1463# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1464# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1465# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1466# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1467# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1468# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1469# OTHER DEALINGS IN THE SOFTWARE.