Coverage for pygeodesy/ecef.py: 95%
507 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-05-02 13:53 -0400
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>}
8into pure Python class L{EcefKarney}, class L{EcefSudano} based on I{John Sudano}'s
9U{paper<https://www.ResearchGate.net/publication/
103709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>},
11class L{EcefUPC} using the I{Universitat Politècnica de Catalunya}'s U{method, page 186
12<https://GSSC.ESA.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf>}, class L{EcefVeness}
13transcoded from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal, Cartesian
14<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}, class L{EcefYou}
15implementing I{Rey-Jer You}'s U{transformations<https://www.ResearchGate.net/publication/240359424>}
16and classes L{EcefFarrell22} and L{EcefFarrell22} from I{Jay A. Farrell}'s U{Table 2.1 and 2.2
17<https://Books.Google.com/books?id=fW4foWASY6wC>}, page 29-30.
19Following is a copy of I{Karney}'s U{Detailed Description
20<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Geocentric.html>}.
22Convert between geodetic coordinates C{lat}-, C{lon}gitude and height C{h} (measured vertically
23from the surface of the ellipsoid) to geocentric C{x}, C{y} and C{z} coordinates, also known as
24I{Earth-Centered, Earth-Fixed} (U{ECEF<https://WikiPedia.org/wiki/ECEF>}).
26The origin of geocentric coordinates is at the center of the earth. The C{z} axis goes thru
27the North pole, C{lat} = 90°. The C{x} axis goes thru C{lat} = 0°, C{lon} = 0°.
29The I{local (cartesian) origin} is at (C{lat0}, C{lon0}, C{height0}). The I{local} C{x} axis points
30East, the I{local} C{y} axis points North and the I{local} C{z} axis is normal to the ellipsoid. The
31plane C{z = -height0} is tangent to the ellipsoid, hence the alternate name I{local tangent plane}.
33Forward conversion from geodetic to geocentric (ECEF) coordinates is straightforward.
35For the reverse transformation we use Hugues Vermeille's U{I{Direct transformation from geocentric
36coordinates to geodetic coordinates}<https://DOI.org/10.1007/s00190-002-0273-6>}, J. Geodesy
37(2002) 76, page 451-454.
39Several changes have been made to ensure that the method returns accurate results for all finite
40inputs (even if h is infinite). The changes are described in Appendix B of C. F. F. Karney
41U{I{Geodesics on an ellipsoid of revolution}<https://ArXiv.org/abs/1102.1215v1>}, Feb. 2011, 85,
42105-117 (U{preprint<https://ArXiv.org/abs/1102.1215v1>}). Vermeille similarly updated his method
43in U{I{An analytical method to transform geocentric into geodetic coordinates}
44<https://DOI.org/10.1007/s00190-010-0419-x>}, J. Geodesy (2011) 85, page 105-117. See U{Geocentric
45coordinates<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>} for more information.
47The errors in these routines are close to round-off. Specifically, for points within 5,000 Km of
48the surface of the ellipsoid (either inside or outside the ellipsoid), the error is bounded by 7
49nm (7 nanometers) for the WGS84 ellipsoid. See U{Geocentric coordinates
50<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>} for further information on the errors.
52@note: The C{reverse} methods of all C{Ecef...} classes return by default C{INT0} as the (geodetic)
53longitude for I{polar} ECEF location C{x == y == 0}. Use keyword argument C{lon00} or property
54C{lon00} to configure that value.
56@see: Module L{ltp} and class L{LocalCartesian}, a transcription of I{Charles Karney}'s C++ class
57U{LocalCartesian<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1LocalCartesian.html>},
58for conversion between geodetic and I{local cartesian} coordinates in a I{local tangent
59plane} as opposed to I{geocentric} (ECEF) ones.
60'''
62from pygeodesy.basics import copysign0, _isin, isscalar, issubclassof, neg, map1, \
63 _xinstanceof, _xsubclassof, typename # _args_kwds_names
64from pygeodesy.constants import EPS, EPS0, EPS02, EPS1, INT0, PI, PI_2, _0_0, \
65 _0_5, _1_0, _1_0_1T, _2_0, _3_0, _4_0, _6_0, _90_0, \
66 _copysign_1_0, isnon0 # PYCHOK used!
67from pygeodesy.datums import _ellipsoidal_datum, _WGS84, a_f2Tuple, _EWGS84
68from pygeodesy.ecefLocals import _EcefLocal
69# from pygeodesy.ellipsoids import a_f2Tuple, _EWGS84 # from .datums
70from pygeodesy.errors import _IndexError, LenError, _ValueError, _TypesError, \
71 _xattr, _xdatum, _xkwds, _xkwds_get
72from pygeodesy.fmath import cbrt, _fdotf, hypot, hypot1, hypot2_, sqrt0
73from pygeodesy.fsums import Fsum, fsumf_, Fmt, unstr
74# from pygeodesy.internals import typename # from .basics
75from pygeodesy.interns import NN, _a_, _C_, _datum_, _ellipsoid_, _f_, _height_, \
76 _lat_, _lon_, _M_, _name_, _singular_, _SPACE_, \
77 _x_, _xyz_, _y_, _z_
78from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
79from pygeodesy.named import _name__, _name1__, _NamedBase, _NamedTuple, _Pass, _xnamed
80from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \
81 PhiLam2Tuple, Vector3Tuple, Vector4Tuple
82from pygeodesy.props import deprecated_method, deprecated_property, Property_RO, \
83 property_RO, property_ROver
84# from pygeodesy.streprs import Fmt, unstr # from .fsums
85from pygeodesy.units import _isRadius, Degrees, Degrees_, Height, Int, Lam, Lat, \
86 Lon, Meter, Phi, Scalar, Scalar_
87from pygeodesy.utily import atan1, atan1d, atan2, atan2d, degrees90, degrees180, \
88 sincos2, sincos2_, sincos2d_
89# from pygeodesy.vector3d import Vector3d # _MODS
91from math import cos, degrees, fabs, radians, sqrt
93__all__ = _ALL_LAZY.ecef
94__version__ = '26.03.28'
96_Ecef_ = 'Ecef'
97_prolate_ = 'prolate'
98_TOL = 1.e-12 # degrees > 1.e-14
99_TRIPS = 33 # 8..9 sufficient
100_xyz_y_z = _xyz_, _y_, _z_ # _args_kwds_names(_xyzn4)[:3]
103class EcefError(_ValueError):
104 '''An ECEF or C{Ecef*} related issue.
105 '''
106 pass
109class _EcefBase(_NamedBase):
110 '''(INTERNAL) Base class for C{Ecef*} conversion classes.
111 '''
112 _datum = _WGS84
113 _E = _EWGS84
114 _isYou = False
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: A (non-prolate) ellipsoid (L{Ellipsoid}, L{Ellipsoid2},
121 L{Datum} or L{a_f2Tuple}) or C{scalar} ellipsoid's
122 equatorial radius (C{meter}).
123 @kwarg f: C{None} or the ellipsoid flattening (C{scalar}), required
124 for C{scalar} B{C{a_ellipsoid}}, C{B{f}=0} represents a
125 sphere, negative B{C{f}} a prolate ellipsoid.
126 @kwarg lon00: An arbitrary, I{"polar"} longitude (C{degrees}), see the
127 C{reverse} method.
128 @kwarg name: Optional C{B{name}=NN} (C{str}).
130 @raise EcefError: If B{C{a_ellipsoid}} not L{Ellipsoid}, L{Ellipsoid2},
131 L{Datum} or L{a_f2Tuple} or C{scalar} or B{C{f}} not
132 C{scalar} or if C{scalar} B{C{a_ellipsoid}} not positive
133 or B{C{f}} not less than 1.0.
134 '''
135 try:
136 E = a_ellipsoid
137 if f is None:
138 pass
139 elif _isRadius(E) and isscalar(f):
140 E = a_f2Tuple(E, f)
141 else:
142 raise ValueError() # _invalid_
144 if not _isin(E, _EWGS84, _WGS84):
145 d = _ellipsoidal_datum(E, **name)
146 E = d.ellipsoid
147 if E.a < EPS or E.f > EPS1:
148 raise ValueError() # _invalid_
149 self._datum = d
150 self._E = E
152 except (TypeError, ValueError) as x:
153 t = unstr(self.classname, a=a_ellipsoid, f=f)
154 raise EcefError(_SPACE_(t, _ellipsoid_), cause=x)
156 if name:
157 self.name = name
158 if lon00 is not INT0:
159 self.lon00 = lon00
161 def __eq__(self, other):
162 '''Compare this and an other Ecef.
164 @arg other: The other ecef (C{Ecef*}).
166 @return: C{True} if equal, C{False} otherwise.
167 '''
168 return other is self or (isinstance(other, type(self)) and
169 other.ellipsoid == self.ellipsoid)
171 @Property_RO
172 def datum(self):
173 '''Get the datum (L{Datum}).
174 '''
175 return self._datum
177 @Property_RO
178 def ellipsoid(self):
179 '''Get the ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
180 '''
181 return self._E
183 @Property_RO
184 def equatoradius(self):
185 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}).
186 '''
187 return self.ellipsoid.a
189 a = equatorialRadius = equatoradius # Karney property
191 @Property_RO
192 def flattening(self): # Karney property
193 '''Get the C{ellipsoid}'s flattening (C{scalar}), positive for
194 I{oblate}, negative for I{prolate} or C{0} for I{near-spherical}.
195 '''
196 return self.ellipsoid.f
198 f = flattening
200 def _forward(self, lat, lon, h, name, M=False, _philam=False): # in .ltp.LocalCartesian.forward and -.reset
201 '''(INTERNAL) Common for all C{Ecef*}.
202 '''
203 if _philam: # lat, lon in radians
204 sa, ca, sb, cb = sincos2_(lat, lon)
205 lat = Lat(degrees90( lat), Error=EcefError)
206 lon = Lon(degrees180(lon), Error=EcefError)
207 else:
208 sa, ca, sb, cb = sincos2d_(lat, lon)
210 E = self.ellipsoid
211 n = E.roc1_(sa, ca) if self._isYou else E.roc1_(sa)
212 c = (h + n) * ca
213 x = cb * c
214 y = sb * c
215 z = (h + n * E.e21) * sa
217 m = self._Matrix(sa, ca, sb, cb) if M else None
218 return Ecef9Tuple(x, y, z, lat, lon, h,
219 0, m, self.datum, # C=0, always
220 name=self._name__(name))
222 def forward(self, latlonh, lon=None, height=0, M=False, **name):
223 '''Convert from geodetic C{(lat, lon, height)} to geocentric C{(x, y, z)}.
225 @arg latlonh: Either a C{LatLon}, an L{Ecef9Tuple} or C{scalar}
226 latitude (C{degrees}).
227 @kwarg lon: Optional C{scalar} longitude for C{scalar} B{C{latlonh}}
228 (C{degrees}).
229 @kwarg height: Optional height (C{meter}), vertically above (or below)
230 the surface of the ellipsoid.
231 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
232 @kwarg name: Optional C{B{name}=NN} (C{str}).
234 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
235 geocentric C{(x, y, z)} coordinates for the given geodetic ones
236 C{(lat, lon, height)}, case C{C} 0, optional C{M} (L{EcefMatrix})
237 and C{datum} if available.
239 @raise EcefError: If B{C{latlonh}} not C{LatLon}, L{Ecef9Tuple} or
240 C{scalar} or B{C{lon}} not C{scalar} for C{scalar}
241 B{C{latlonh}} or C{abs(lat)} exceeds 90°.
243 @note: Use method C{.forward_} to specify C{lat} and C{lon} in C{radians}
244 and avoid double angle conversions.
245 '''
246 llhn = _llhn4(latlonh, lon, height, **name)
247 return self._forward(*llhn, M=M)
249 def forward_(self, phi, lam, height=0, M=False, **name):
250 '''Like method C{.forward} except with geodetic lat- and longitude given
251 in I{radians}.
253 @arg phi: Latitude in I{radians} (C{scalar}).
254 @arg lam: Longitude in I{radians} (C{scalar}).
255 @kwarg height: Optional height (C{meter}), vertically above (or below)
256 the surface of the ellipsoid.
257 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
258 @kwarg name: Optional C{B{name}=NN} (C{str}).
260 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
261 with C{lat} set to C{degrees90(B{phi})} and C{lon} to
262 C{degrees180(B{lam})}.
264 @raise EcefError: If B{C{phi}} or B{C{lam}} invalid or not C{scalar}.
265 '''
266 try: # like function C{_llhn4} below
267 plhn = Phi(phi), Lam(lam), Height(height), _name__(name)
268 except (TypeError, ValueError) as x:
269 raise EcefError(phi=phi, lam=lam, height=height, cause=x)
270 return self._forward(*plhn, M=M, _philam=True)
272 @property_ROver
273 def _Geocentrics(self):
274 '''(INTERNAL) Get the valid geocentric classes. I{once}.
275 '''
276 return (Ecef9Tuple, # overwrite property_ROver
277 _MODS.vector3d.Vector3d) # _MODS.cartesianBase.CartesianBase
279 @property
280 def lon00(self):
281 '''Get the I{"polar"} longitude (C{degrees}), see method C{reverse}.
282 '''
283 return self._lon00
285 @lon00.setter # PYCHOK setter!
286 def lon00(self, lon00):
287 '''Set the I{"polar"} longitude (C{degrees}), see method C{reverse}.
288 '''
289 self._lon00 = Degrees(lon00=lon00)
291 def _Matrix(self, *sa_ca_sb_cb):
292 '''(INTERNAL) Create a rotation L{EcefMatrix}.
293 '''
294 return self._xnamed(EcefMatrix(*sa_ca_sb_cb))
296 def _polon(self, y, x, p, **lon00_name):
297 '''(INTERNAL) Handle I{"polar"} longitude.
298 '''
299 return atan2d(y, x) if p else _xkwds_get(lon00_name, lon00=self.lon00)
301 def reverse(self, xyz, y=None, z=None, M=False, **tol_lon00_name): # PYCHOK no cover
302 '''I{Must be overloaded}.'''
303 self._notOverloaded(xyz, y=y, z=z, M=M, **tol_lon00_name)
305 def _reversError(self, x, y, z, r, tol):
306 '''(INTERNAL) Convergence error.
307 '''
308 t = unstr(self.reverse, x=x, y=y, z=z)
309 return EcefError(t, txt=Fmt.no_convergence(r, tol))
311 def toStr(self, prec=9, **unused): # PYCHOK signature
312 '''Return this C{Ecef*} as a string.
314 @kwarg prec: Precision, number of decimal digits (0..9).
316 @return: This C{Ecef*} (C{str}).
317 '''
318 return self.attrs(_a_, _f_, _datum_, _name_, prec=prec) # _ellipsoid_
320 def _xyzllhCpn9(self, xyz, y, z, **lon00_name):
321 '''(INTERNAL) Get C{x, y, z} and determine case C{C}, C{lat}, C{lon}, etc.
322 '''
323 x, y, z, n = _xyzn4(xyz, y, z, self._Geocentrics, **lon00_name)
325 s, c, p = _norm3(y, x) # distance to polar axis
326 if p < EPS0: # polar
327 lat = copysign0(_90_0, z)
328 h = fabs(z) - self.ellipsoid.b
329 C = 2
330 p = 0 # force lon00
331 elif fabs(z) < EPS0: # equatorial
332 lat = _0_0
333 h = p - self.ellipsoid.a
334 C = 3
335 else: # see _norm7 below, EcefKarney
336 h = hypot(z, p) # distance to earth center
337 if h > self.ellipsoid._heightMax:
338 lat = atan1d(z / h, p / h)
339 C = 4 # too high
340# p *= _0_5
341 else: # pass h for EcefVeness
342 lat = None # -90..90
343 C = 1 # normal
344 lon = self._polon(s, c, p, **lon00_name)
345 return x, y, z, lat, lon, h, C, p, self._name__(n)
348class EcefFarrell21(_EcefBase):
349 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)
350 coordinates based on I{Jay A. Farrell}'s U{Table 2.1<https://Books.Google.com/
351 books?id=fW4foWASY6wC>}, page 29.
352 '''
354 def reverse(self, xyz, y=None, z=None, M=None, **lon00_name): # PYCHOK unused M
355 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using
356 I{Farrell}'s U{Table 2.1<https://Books.Google.com/books?id=fW4foWASY6wC>},
357 page 29, aka the I{Heikkinen application} of U{Ferrari's solution
358 <https://WikiPedia.org/wiki/Geographic_coordinate_conversion>}.
360 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
361 coordinate (C{meter}).
362 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
363 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
364 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
365 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument
366 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude
367 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.
369 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
370 geodetic coordinates C{(lat, lon, height)} for the given geocentric
371 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum}.
373 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
374 not C{scalar} for C{scalar} B{C{xyz}} or C{sqrt} domain or
375 zero division error.
377 @see: L{EcefFarrell22} and L{EcefVeness}.
378 '''
379 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)
380 if lat is None:
381 E = self.ellipsoid
382 a = E.a
383 a2 = E.a2
384 b2 = E.b2
385 e2 = E.e2
386 e2_ = E.e2abs * E.a2_b2 # (E.e * E.a_b)**2 = 0.0820944... WGS84
387 e4 = E.e4
389 try: # names as page 29
390 z2 = z**2
391 ez = z2 * (_1_0 - e2) # E.e2s2(z)
393 p2 = p**2
394 G = p2 + ez - e2 * (a2 - b2) # p2 + ez - e4 * a2
395 F = b2 * z2 * 54
396 c = e4 * p2 * F / G**3
397 s = sqrt(c * (c + _2_0))
398 c = cbrt(s + c + _1_0)
399 G *= fsumf_(c, _1_0, _1_0 / c) # k
400 P = F / (G**2 * _3_0)
401 Q = sqrt(_2_0 * e4 * P + _1_0)
402 Q1 = Q + _1_0
403 s = fsumf_(a2 * (Q1 / Q) * _0_5,
404 -P * ez / (Q * Q1),
405 -P * p2 * _0_5)
406 r = p * P * e2 / Q1 - sqrt(s)
407 r = p + r * e2
408 v = b2 / (sqrt(r**2 + ez) * a) # z0 / z
410 h = hypot(r, z) * (_1_0 - v)
411 lat = atan1d((e2_ * v + _1_0) * z, p)
412 # note, phi and lam are swapped on page 29
414 except (ValueError, ZeroDivisionError) as X:
415 raise EcefError(x=x, y=y, z=z, cause=X)
417 return Ecef9Tuple(x, y, z, lat, lon, h,
418 C, None, self.datum, name=name)
421class EcefFarrell22(_EcefBase):
422 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)
423 coordinates based on I{Jay A. Farrell}'s U{Table 2.2<https://Books.Google.com/
424 books?id=fW4foWASY6wC>}, page 30.
425 '''
427 def reverse(self, xyz, y=None, z=None, M=None, **lon00_name): # PYCHOK unused M
428 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using
429 I{Farrell}'s U{Table 2.2<https://Books.Google.com/books?id=fW4foWASY6wC>},
430 page 30.
432 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
433 coordinate (C{meter}).
434 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
435 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
436 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
437 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument
438 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude
439 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.
441 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
442 geodetic coordinates C{(lat, lon, height)} for the given geocentric
443 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum}.
445 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
446 not C{scalar} for C{scalar} B{C{xyz}} or C{sqrt} domain or
447 zero division error.
449 @see: L{EcefFarrell21} and L{EcefVeness}.
450 '''
451 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)
452 if lat is None:
453 E = self.ellipsoid
454 a, b = E.a, E.b
455 s, c, _ = _norm3(z * a, p * b) # Bowring
456 s, c, _ = _norm3(z + s**3 * b * E.e22,
457 p - c**3 * a * E.e2)
458 lat = atan1d(s, c)
459 h = (p / fabs(c) - (E.roc1_(s) if s else a)) if c else (fabs(z) - b)
460 # note, phi and lam are swapped on page 30
462 return Ecef9Tuple(x, y, z, lat, lon, h,
463 C, None, self.datum, name=name)
466class EcefKarney(_EcefBase):
467 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)
468 coordinates transcoded from I{Karney}'s C++ U{Geocentric<https://GeographicLib.SourceForge.io/
469 C++/doc/classGeographicLib_1_1Geocentric.html>} methods.
471 @note: On methods C{.forward} and C{.forwar_}, let C{v} be a unit vector located
472 at C{(lat, lon, h)}. We can express C{v} as column vectors in one of two
473 ways, C{v1} in East, North, Up (ENU) coordinates (where the components are
474 relative to a local coordinate system at C{C(lat0, lon0, h0)}) or as C{v0}
475 in geocentric C{x, y, z} coordinates. Then, M{v0 = M ⋅ v1} where C{M} is
476 the rotation matrix.
477 '''
479 def reverse(self, xyz, y=None, z=None, M=False, **lon00_name):
480 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}.
482 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
483 coordinate (C{meter}).
484 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
485 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
486 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
487 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument
488 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude
489 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.
491 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
492 geodetic coordinates C{(lat, lon, height)} for the given geocentric
493 ones C{(x, y, z)}, case C{C}, optional C{M} (L{EcefMatrix}) and
494 C{datum} if available.
496 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
497 not C{scalar} for C{scalar} B{C{xyz}}.
499 @note: In general, there are multiple solutions and the result which minimizes
500 C{height} is returned, i.e., the C{(lat, lon)} corresponding to the
501 closest point on the ellipsoid. If there are still multiple solutions
502 with different latitudes (applies only if C{z} = 0), then the solution
503 with C{lat} > 0 is returned. If there are still multiple solutions with
504 different longitudes (applies only if C{x} = C{y} = 0), then C{lon00} is
505 returned. The returned C{lon} is in the range [−180°, 180°] and C{height}
506 is not below M{−E.a * (1 − E.e2) / sqrt(1 − E.e2 * sin(lat)**2)}. Like
507 C{forward} above, M{v1 = Transpose(M) ⋅ v0}.
508 '''
509 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, **lon00_name)
511 E = self.ellipsoid
512 f = E.f
514 sa, ca, sb, cb, R, h, C = _norm7(y, x, z, E)
515 if C: # PYCHOK no cover
516 pass # too high, far
517 elif E.e4: # E.isEllipsoidal
518 # Treat prolate spheroids by swapping R and Z here and by
519 # switching the arguments to phi = atan2(...) at the end.
520 p = (R / E.a)**2
521 q = (z / E.a)**2 * E.e21
522 if f < 0:
523 p, q = q, p
524 r = fsumf_(p, q, -E.e4)
525 e = E.e4 * q
526 if e or r > 0:
527 # Avoid possible division by zero when r = 0 by multiplying
528 # equations for s and t by r^3 and r, respectively.
529 s = d = e * p / _4_0 # s = r^3 * s
530 u = r = r / _6_0
531 r2 = r**2
532 r3 = r2 * r
533 t3 = r3 + s
534 d *= t3 + r3
535 if d < 0:
536 # t is complex, but the way u is defined, the result is real.
537 # There are three possible cube roots. We choose the root
538 # which avoids cancellation. Note, d < 0 implies r < 0.
539 u += cos(atan2(sqrt(-d), -t3) / _3_0) * r * _2_0
540 else:
541 # Pick the sign on the sqrt to maximize abs(t3). This
542 # minimizes loss of precision due to cancellation. The
543 # result is unchanged because of the way the t is used
544 # in definition of u.
545 if d > 0:
546 t3 += copysign0(sqrt(d), t3) # t3 = (r * t)^3
547 # N.B. cbrt always returns the real root, cbrt(-8) = -2.
548 t = cbrt(t3) # t = r * t
549 if t: # t can be zero; but then r2 / t -> 0.
550 u = fsumf_(u, t, r2 / t)
551 v = sqrt(e + u**2) # guaranteed positive
552 # Avoid loss of accuracy when u < 0. Underflow doesn't occur in
553 # E.e4 * q / (v - u) because u ~ e^4 when q is small and u < 0.
554 u = (e / (v - u)) if u < 0 else (u + v) # u+v, guaranteed positive
555 # Need to guard against w going negative due to roundoff in u - q.
556 w = E.e2abs * (u - q) / (_2_0 * v)
557 # Rearrange expression for k to avoid loss of accuracy due to
558 # subtraction. Division by 0 not possible because u > 0, w >= 0.
559 k1 = k2 = (u / (sqrt(u + w**2) + w)) if w > 0 else sqrt(u)
560 if f < 0:
561 k1 -= E.e2
562 else:
563 k2 += E.e2
564 sa, ca, h = _norm3(z / k1, R / k2)
565 h *= k1 - E.e21
566 C = 1
568 else: # e = E.e4 * q == 0 and r <= 0
569 # This leads to k = 0 (oblate, equatorial plane) and k + E.e^2 = 0
570 # (prolate, rotation axis) and the generation of 0/0 in the general
571 # formulas for phi and h, using the general formula and division
572 # by 0 in formula for h. Handle this case by taking the limits:
573 # f > 0: z -> 0, k -> E.e2 * sqrt(q) / sqrt(E.e4 - p)
574 # f < 0: r -> 0, k + E.e2 -> -E.e2 * sqrt(q) / sqrt(E.e4 - p)
575 q = E.e4 - p
576 if f < 0:
577 p, q = q, p
578 e = E.a
579 else:
580 e = E.b2_a
581 sa, ca, h = _norm3(sqrt(q * E._1_e21), sqrt(p))
582 if z < 0: # for tiny negative z, not for prolate
583 sa = neg(sa)
584 h *= neg(e / E.e2abs)
585 C = 3
587 else: # E.e4 == 0, spherical case
588 # Dealing with underflow in the general case with E.e2 = 0 is
589 # difficult. Origin maps to North pole, same as with ellipsoid.
590 sa, ca, _ = _norm3((z if h else _1_0), R)
591 h -= E.a
592 C = 2
594 # lon00 <https://GitHub.com/mrJean1/PyGeodesy/issues/77>
595 lon = self._polon(sb, cb, R, **lon00_name)
596 m = self._Matrix(sa, ca, sb, cb) if M else None
597 return Ecef9Tuple(x, y, z, atan1d(sa, ca), lon, h,
598 C, m, self.datum, name=self._name__(name))
601class EcefSudano(_EcefBase):
602 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
603 based on I{John J. Sudano}'s U{paper<https://www.ResearchGate.net/publication/3709199>}.
604 '''
605 _tol = EPS # DEPRECATED
607 def reverse(self, xyz, y=None, z=None, M=None, tol=EPS, **lon00_name): # PYCHOK unused M
608 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using
609 I{Sudano}'s U{iterative method<https://www.ResearchGate.net/publication/3709199>}.
611 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
612 coordinate (C{meter}).
613 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
614 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
615 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
616 @kwarg tol: Convergence tolerance for C{sin}(latitude) (C{scalar}).
617 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument
618 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude
619 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.
621 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with geodetic
622 coordinates C{(lat, lon, height)} for the given geocentric ones C{(x, y, z)},
623 case C{C}, C{M=None} always and C{datum} if available.
625 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
626 not C{scalar} for C{scalar} B{C{xyz}} or no convergence.
628 @see: Class L{EcefUPC}.
629 '''
630 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)
631 if lat is None:
632 E = self.ellipsoid
633 e = E.e2 * E.a
634 d = e - p
636 sa, ca, _ = _norm3(fabs(z), p * E.e21)
637 _S2 = Fsum(sa).fsum2f_
638 tol = Scalar_(tol=tol, low=self.tolerance, Error=EcefError)
639 # Sudano's Eq (A-6) and (A-7) refactored/reduced,
640 # replacing Rn from Eq (A-4) with n = E.a / ca:
641 # N = ca**2 * ((z + E.e2 * n * sa) * ca - p * sa)
642 # = ca**2 * (z * ca + E.e2 * E.a * sa - p * sa)
643 # = ca**2 * (z * ca + (E.e2 * E.a - p) * sa)
644 # D = ca**3 * (E.e2 * n / E.e2s2(sa)) - p
645 # = ca**2 * (E.e2 * E.a / E.e2s2(sa) - p / ca**2)
646 # N / D = (z * ca + (E.e2 * E.a - p) * sa) /
647 # (E.e2 * E.a / E.e2s2(sa) - p / ca**2)
648 for i in range(1, _TRIPS): # 6+ max
649 ca2 = _1_0 - sa**2
650 if ca2 < EPS02:
651 break
652 D = p / ca2 - e / E.e2s2(sa)
653 if fabs(D) < EPS0:
654 break
655 ca = sqrt(ca2)
656 sa, D = _S2(z * ca / D, d * sa / D)
657 if fabs(D) < tol:
658 break
659 else: # PYCHOK no cover
660 raise self._reversError(x, y, z, fabs(D), tol)
662 sa = copysign0(sa, z)
663 lat = atan1d(sa, ca)
664 # h = (fabs(z) + p - E.a * cos(a + E.e21) * sa / ca) / (ca + sa)
665 # Sudano's Eq (7) doesn't produce the correct height, ...
666 h = E._heightB(sa, ca, z, p) # ... use Veness' (Bowring eqn 7)
667 else:
668 i = 0
669 return Ecef9Tuple(x, y, z, lat, lon, h,
670 C, None, self.datum, # M=None
671 iteration=i, name=name)
673 @deprecated_property
674 def tolerance(self):
675 '''DEPRECATED on 2025.08.22, use keyword argument C{tol}.'''
676 return self._tol
678 @tolerance.setter # PYCHOK setter!
679 def tolerance(self, tol):
680 self._tol = Scalar_(tolerance=tol, low=EPS, Error=EcefError)
683class EcefUPC(_EcefBase):
684 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates based on
685 I{UPC}'s U{method<https://GSSC.ESA.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion>}.
686 '''
688 def reverse(self, xyz, y=None, z=None, M=None, tol=_TOL, **lon00_name): # PYCHOK unused M
689 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using I{UPC}'s
690 U{iterative method<https://GSSC.ESA.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf>}, page 186.
692 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
693 coordinate (C{meter}).
694 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
695 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
696 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
697 @kwarg tol: Convergence tolerance for the latitude (C{degrees}).
698 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument
699 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude
700 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.
702 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with geodetic
703 coordinates C{(lat, lon, height)} for the given geocentric ones C{(x, y, z)},
704 case C{C}, C{M=None} always and C{datum} if available.
706 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
707 not C{scalar} for C{scalar} B{C{xyz}} or no convergence.
709 @see: Class L{EcefSudano}.
710 '''
711 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)
712 if lat is None:
713 t = Degrees_(tol=tol, low=EPS, Error=EcefError).toRadians()
714 E = self.ellipsoid
715 a = E.a
716 e2 = E.e2 # signed
718 z_ = fabs(z)
719 ph_ = atan1(z_, E.e21 * p)
720 for i in range(1, _TRIPS): # 5..6 max
721 s, c = sincos2(ph_)
722 N = a / sqrt(_1_0 - s**2 * e2)
723 ph = atan1(z_, p - N * c * e2)
724 r = fabs(ph - ph_)
725 if r < t:
726 lat = copysign0(degrees(ph), z)
727 h = p / c - N
728 break
729 ph_ = ph
730 else: # PYCHOK no cover
731 raise self._reversError(x, y, z, degrees(r), tol)
732 else:
733 i = 0
734 return Ecef9Tuple(x, y, z, lat, lon, h,
735 C, None, self.datum, # M=None
736 iteration=i, name=name)
739class EcefVeness(_EcefBase):
740 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
741 transcoded from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal, Cartesian<https://
742 www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}.
744 @see: U{I{A Guide to Coordinate Systems in Great Britain}<https://www.OrdnanceSurvey.co.UK/
745 documents/resources/guide-coordinate-systems-great-britain.pdf>}, section I{B) Converting
746 between 3D Cartesian and ellipsoidal latitude, longitude and height coordinates}.
747 '''
749 def reverse(self, xyz, y=None, z=None, M=None, **lon00_name): # PYCHOK unused M
750 '''Conversion from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}
751 transcoded from I{Chris Veness}' U{JavaScript<https://www.Movable-Type.co.UK/
752 scripts/geodesy/docs/latlon-ellipsoidal.js.html>}.
754 Uses B. R. Bowring’s formulation for μm precision in concise form U{I{The accuracy
755 of geodetic latitude and height equations}<https://www.ResearchGate.net/publication/
756 233668213>}, Survey Review, Vol 28, 218, Oct 1985.
758 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
759 coordinate (C{meter}).
760 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
761 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
762 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
763 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument
764 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude
765 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.
767 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
768 geodetic coordinates C{(lat, lon, height)} for the given geocentric
769 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum} if available.
771 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
772 not C{scalar} for C{scalar} B{C{xyz}}.
774 @see: Toms, Ralph M. U{I{An Efficient Algorithm for Geocentric to Geodetic
775 Coordinate Conversion}<https://www.OSTI.gov/scitech/biblio/110235>},
776 Sept 1995 and U{I{An Improved Algorithm for Geocentric to Geodetic
777 Coordinate Conversion}<https://www.OSTI.gov/scitech/servlets/purl/231228>},
778 Apr 1996, both from Lawrence Livermore National Laboratory (LLNL).
779 '''
780 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)
781 if lat is None:
782 E = self.ellipsoid
783 a = E.a
784 B = E.b * E.e22
785 # parametric latitude (Bowring eqn 17, replaced)
786 t = (E.b * z) / (a * p) * (_1_0 + B / h) # h = hypot(z, p)
787 c = _1_0 / hypot1(t)
788 s = c * t
789 # geodetic latitude (Bowring eqn 18)
790 s, c, _ = _norm3(z + s**3 * B,
791 p - c**3 * a * E.e2)
792 lat = atan1d(s, c)
793 h = E._heightB(s, c, z, p) # height (Bowring eqn 7)
795 return Ecef9Tuple(x, y, z, lat, lon, h,
796 C, None, self.datum, name=name)
799class EcefYou(_EcefBase):
800 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
801 using I{Rey-Jer You}'s U{transformation<https://www.ResearchGate.net/publication/240359424>}
802 for I{non-prolate} ellipsoids.
804 @see: Featherstone, W.E., Claessens, S.J. U{I{Closed-form transformation between geodetic and
805 ellipsoidal coordinates}<https://Espace.Curtin.edu.AU/bitstream/handle/20.500.11937/
806 11589/115114_9021_geod2ellip_final.pdf>} Studia Geophysica et Geodaetica, 2008, 52,
807 pages 1-18 and U{PyMap3D <https://PyPI.org/project/pymap3d>}.
808 '''
809 _isYou = True
811 def __init__(self, a_ellipsoid=_EWGS84, f=None, **lon00_name): # PYCHOK signature
812 _EcefBase.__init__(self, a_ellipsoid, f=f, **lon00_name) # inherited documentation
814 E = self.ellipsoid
815 e2 = E.a2 - E.b2
816 if e2 < 0 or E.f < 0:
817 raise EcefError(ellipsoid=E, txt=_prolate_)
818 self._ee2 = sqrt0(e2), e2
820 def reverse(self, xyz, y=None, z=None, M=None, **lon00_name): # PYCHOK unused M
821 '''Convert geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}
822 using I{Rey-Jer You}'s transformation.
824 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
825 coordinate (C{meter}).
826 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
827 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
828 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
829 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument
830 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude
831 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.
833 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
834 geodetic coordinates C{(lat, lon, height)} for the given geocentric
835 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum} if available.
837 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or
838 B{C{z}} not C{scalar} for C{scalar} B{C{xyz}} or the
839 ellipsoid is I{prolate}.
840 '''
841 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)
842 if lat is None:
843 E = self.ellipsoid
844 a, b = E.a, E.b
845 e, e2 = self._ee2
847 u = hypot2_(x, y, z) - e2
848 u += hypot(u, e * z * _2_0)
849 u *= _0_5
850 if u > EPS02:
851 u = sqrt(u)
852 q = hypot(u, e)
853 B = atan1(q * z, u * p) # beta0 = atan(q / u * z / p)
854 sB, cB = sincos2(B)
855 if cB and sB:
856 q *= a
857 d = (q / cB - e2 * cB) / sB
858 if isnon0(d):
859 B += fsumf_(u * b, -q, e2) / d
860 sB, cB = sincos2(B)
861 elif u < (-EPS02):
862 raise EcefError(x=x, y=y, z=z, u=u, txt=_singular_)
863 else: # near polar # PYCHOK no cover
864 sB, cB, C = _copysign_1_0(z), _0_0, 2
866 lat = atan1d( a * sB, b * cB) # atan(E.a_b * tan(B))
867 h = hypot(p - a * cB, z - b * sB)
868 if hypot2_(x, y, z * E.a_b) < E.a2: # or lat < 0 or z < 0
869 h = neg(h) # inside ellipsoid
871 return Ecef9Tuple(x, y, z, lat, lon, h,
872 C, None, self.datum, name=name)
875class EcefMatrix(_NamedTuple):
876 '''A rotation matrix known as I{East-North-Up (ENU) to ECEF}.
878 @see: U{From ENU to ECEF<https://WikiPedia.org/wiki/
879 Geographic_coordinate_conversion#From_ECEF_to_ENU>} and
880 U{Issue #74<https://Github.com/mrJean1/PyGeodesy/issues/74>}.
881 '''
882 _Names_ = ('_0_0_', '_0_1_', '_0_2_', # row-order
883 '_1_0_', '_1_1_', '_1_2_',
884 '_2_0_', '_2_1_', '_2_2_')
885 _Units_ = (Scalar,) * len(_Names_)
887 def _validate(self, **unused): # PYCHOK unused
888 '''(INTERNAL) Allow C{_Names_} with leading underscore.
889 '''
890 _NamedTuple._validate(self, underOK=True)
892 def __new__(cls, sa, ca, sb, cb, *_more, **name):
893 '''New L{EcefMatrix} matrix.
895 @arg sa: C{sin(phi)} (C{float}).
896 @arg ca: C{cos(phi)} (C{float}).
897 @arg sb: C{sin(lambda)} (C{float}).
898 @arg cb: C{cos(lambda)} (C{float}).
899 @arg _more: (INTERNAL) from C{.multiply}.
901 @raise EcefError: If B{C{sa}}, B{C{ca}}, B{C{sb}} or
902 B{C{cb}} outside M{[-1.0, +1.0]}.
903 '''
904 t = sa, ca, sb, cb
905 if _more: # all 9 matrix elements ...
906 t += _more # ... from .multiply
908 elif max(map(fabs, t)) > _1_0:
909 raise EcefError(unstr(EcefMatrix, *t))
911 else: # build matrix from the following quaternion operations
912 # qrot(lam, [0,0,1]) * qrot(phi, [0,-1,0]) * [1,1,1,1]/2
913 # or
914 # qrot(pi/2 + lam, [0,0,1]) * qrot(-pi/2 + phi, [-1,0,0])
915 # where
916 # qrot(t,v) = [cos(t/2), sin(t/2)*v[1], sin(t/2)*v[2], sin(t/2)*v[3]]
918 # Local X axis (East) in geocentric coords
919 # M[0] = -slam; M[3] = clam; M[6] = 0;
920 # Local Y axis (North) in geocentric coords
921 # M[1] = -clam * sphi; M[4] = -slam * sphi; M[7] = cphi;
922 # Local Z axis (Up) in geocentric coords
923 # M[2] = clam * cphi; M[5] = slam * cphi; M[8] = sphi;
924 t = (-sb, -cb * sa, cb * ca,
925 cb, -sb * sa, sb * ca,
926 _0_0, ca, sa)
928 return _NamedTuple.__new__(cls, *t, **name)
930 def column(self, column):
931 '''Get this matrix' B{C{column}} 0, 1 or 2 as C{3-tuple}.
932 '''
933 if 0 <= column < 3:
934 return self[column::3]
935 raise _IndexError(column=column)
937 @property_RO
938 def _columns(self):
939 for c in range(3):
940 yield self[c::3]
942 def copy(self, **unused): # PYCHOK signature
943 '''Make a shallow or deep copy of this instance.
945 @return: The copy (C{This class} or subclass thereof).
946 '''
947 return self.classof(*self)
949 __copy__ = __deepcopy__ = copy
951 @Property_RO
952 def matrix3(self):
953 '''Get this matrix' rows (C{3-tuple} of 3 C{3-tuple}s).
954 '''
955 return tuple(self._rows)
957 @Property_RO
958 def matrixTransposed3(self):
959 '''Get this matrix' I{Transposed} rows (C{3-tuple} of 3 C{3-tuple}s).
960 '''
961 return tuple(self._columns)
963 def multiply(self, other):
964 '''Matrix multiply M{M0' ⋅ M} this matrix I{Transposed} with an other matrix.
966 @arg other: The other matrix (L{EcefMatrix}).
968 @return: The matrix product (L{EcefMatrix}).
970 @raise TypeError: If B{C{other}} is not an L{EcefMatrix}.
971 '''
972 _xinstanceof(EcefMatrix, other=other)
973 # like LocalCartesian.MatrixMultiply, C{self.matrixTransposed3 X other.matrix3}
974 # <https://GeographicLib.SourceForge.io/C++/doc/LocalCartesian_8cpp_source.html>
975 X = (_fdotf(t, *c) for t in self._columns for c in other._columns)
976 return _xnamed(EcefMatrix(*X), typename(EcefMatrix.multiply))
978 def rotate(self, xyz, *xyz0):
979 '''Forward rotation M{M0' ⋅ ([x, y, z] - [x0, y0, z0])'}.
981 @arg xyz: Local C{(x, y, z)} coordinates (C{3-tuple}).
982 @arg xyz0: Optional, local C{(x0, y0, z0)} origin (C{3-tuple}).
984 @return: Rotated C{(x, y, z)} location (C{3-tuple}).
986 @raise LenError: Unequal C{len(B{xyz})} and C{len(B{xyz0})}.
987 '''
988 if xyz0:
989 if len(xyz0) != len(xyz):
990 raise LenError(self.rotate, xyz0=len(xyz0), xyz=len(xyz))
991 xyz = tuple(s - s0 for s, s0 in zip(xyz, xyz0))
993 # x' = M[0] * x + M[3] * y + M[6] * z
994 # y' = M[1] * x + M[4] * y + M[7] * z
995 # z' = M[2] * x + M[5] * y + M[8] * z
996 return tuple(_fdotf(xyz, *c) for c in self._columns)
998 def row(self, row):
999 '''Get this matrix' B{C{row}} 0, 1 or 2 as C{3-tuple}.
1000 '''
1001 if 0 <= row < 3:
1002 r = row * 3
1003 return self[r:r+3]
1004 raise _IndexError(row=row)
1006 @property_RO
1007 def _rows(self):
1008 for r in (0, 3, 6):
1009 yield self[r:r+3]
1011 def unrotate(self, xyz, *xyz0):
1012 '''Inverse rotation M{[x0, y0, z0] + M0 ⋅ [x,y,z]'}.
1014 @arg xyz: Local C{(x, y, z)} coordinates (C{3-tuple}).
1015 @arg xyz0: Optional, local C{(x0, y0, z0)} origin (C{3-tuple}).
1017 @return: Unrotated C{(x, y, z)} location (C{3-tuple}).
1019 @raise LenError: Unequal C{len(B{xyz})} and C{len(B{xyz0})}.
1020 '''
1021 if xyz0:
1022 if len(xyz0) != len(xyz):
1023 raise LenError(self.unrotate, xyz0=len(xyz0), xyz=len(xyz))
1024 _xyz = _1_0_1T + xyz
1025 # x' = x0 + M[0] * x + M[1] * y + M[2] * z
1026 # y' = y0 + M[3] * x + M[4] * y + M[5] * z
1027 # z' = z0 + M[6] * x + M[7] * y + M[8] * z
1028 xyz_ = (_fdotf(_xyz, s, *r) for s, r in zip(xyz0, self._rows))
1029 else:
1030 # x' = M[0] * x + M[1] * y + M[2] * z
1031 # y' = M[3] * x + M[4] * y + M[5] * z
1032 # z' = M[6] * x + M[7] * y + M[8] * z
1033 xyz_ = (_fdotf(xyz, *r) for r in self._rows)
1034 return tuple(xyz_)
1037class Ecef9Tuple(_NamedTuple, _EcefLocal):
1038 '''9-Tuple C{(x, y, z, lat, lon, height, C, M, datum)} with I{geocentric} C{x},
1039 C{y} and C{z} plus I{geodetic} C{lat}, C{lon} and C{height}, case C{C} and
1040 optionally, rotation matrix C{M} (L{EcefMatrix}) and C{datum}, with C{lat}
1041 and C{lon} in C{degrees} and C{x}, C{y}, C{z} and C{height} in C{meter},
1042 conventionally. Case C{C=1} means normal, C{C=2} near polar and C{C=3}
1043 equatorial latitude and C{C=4} height exceeds C{heightMax}.
1044 '''
1045 _Names_ = (_x_, _y_, _z_, _lat_, _lon_, _height_, _C_, _M_, _datum_)
1046 _Units_ = ( Meter, Meter, Meter, Lat, Lon, Height, Int, _Pass, _Pass)
1048 @property_ROver
1049 def _CartesianBase(self):
1050 '''(INTERNAL) Get class C{CartesianBase}, I{once}.
1051 '''
1052 return _MODS.cartesianBase.CartesianBase # overwrite property_ROver
1054 @deprecated_method
1055 def convertDatum(self, datum2): # for backward compatibility
1056 '''DEPRECATED, use method L{toDatum}.'''
1057 return self.toDatum(datum2)
1059 @property_RO
1060 def _ecef9(self): # in ._EcefLocal._Ltp_ecef2local
1061 return self
1063 @property_RO
1064 def ellipsoid(self):
1065 '''Get the ellipsoid (L{Ellipsoid}).
1066 '''
1067 return (self.datum or _WGS84).ellipsoid
1069 @Property_RO
1070 def lam(self):
1071 '''Get the longitude in C{radians} (C{float}).
1072 '''
1073 return self.philam.lam
1075 @Property_RO
1076 def lamVermeille(self):
1077 '''Get the longitude in C{radians} M{[-PI*3/2..+PI*3/2]} after U{Vermeille
1078 <https://Search.ProQuest.com/docview/639493848>} (2004), page 95.
1080 @see: U{Karney<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>},
1081 U{Vermeille<https://Search.ProQuest.com/docview/847292978>} 2011, pp 112-113, 116
1082 and U{Featherstone, et.al.<https://Search.ProQuest.com/docview/872827242>}, page 7.
1083 '''
1084 x, y = self.x, self.y
1085 a = fabs(y)
1086 if a > EPS0:
1087 r = PI_2 - atan2(x, hypot(x, a) + a) * _2_0
1088 if y < 0:
1089 r = -r
1090 else: # y == 0
1091 r = PI if x < 0 else _0_0
1092 return Lam(Vermeille=r)
1094 @Property_RO
1095 def latlon(self):
1096 '''Get the lat-, longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).
1097 '''
1098 return LatLon2Tuple(self.lat, self.lon, name=self.name)
1100 @Property_RO
1101 def latlonheight(self):
1102 '''Get the lat-, longitude in C{degrees} and height (L{LatLon3Tuple}C{(lat, lon, height)}).
1103 '''
1104 return self.latlon.to3Tuple(self.height)
1106 @Property_RO
1107 def latlonheightdatum(self):
1108 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}).
1109 '''
1110 return self.latlonheight.to4Tuple(self.datum)
1112 @Property_RO
1113 def latlonVermeille(self):
1114 '''Get the latitude and I{Vermeille} longitude in C{degrees [-225..+225]} (L{LatLon2Tuple}C{(lat, lon)}).
1116 @see: Property C{lonVermeille}.
1117 '''
1118 return LatLon2Tuple(self.lat, self.lonVermeille, name=self.name)
1120 @Property_RO
1121 def lonVermeille(self):
1122 '''Get the longitude in C{degrees [-225..+225]} after U{Vermeille
1123 <https://Search.ProQuest.com/docview/639493848>} 2004, page 95.
1125 @see: Property C{lamVermeille}.
1126 '''
1127 return Lon(Vermeille=degrees(self.lamVermeille))
1129 @Property_RO
1130 def Mx(self):
1131 '''Compute rotation matrix (L{EcefMatrix}), seperate from C{M}.
1132 '''
1133 sa, ca, sb, cb, _, _, _ = _norm7(self.y, self.x, self.z, self.ellipsoid)
1134 return EcefMatrix(sa, ca, sb, cb, name=self.name)
1136 @Property_RO
1137 def phi(self):
1138 '''Get the latitude in C{radians} (C{float}).
1139 '''
1140 return self.philam.phi
1142 @Property_RO
1143 def philam(self):
1144 '''Get the lat-, longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).
1145 '''
1146 return PhiLam2Tuple(radians(self.lat), radians(self.lon), name=self.name)
1148 @Property_RO
1149 def philamheight(self):
1150 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
1151 '''
1152 return self.philam.to3Tuple(self.height)
1154 @Property_RO
1155 def philamheightdatum(self):
1156 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLam4Tuple}C{(phi, lam, height, datum)}).
1157 '''
1158 return self.philamheight.to4Tuple(self.datum)
1160 @Property_RO
1161 def philamVermeille(self):
1162 '''Get the latitude and I{Vermeille} longitude in C{radians [-PI*3/2..+PI*3/2]} (L{PhiLam2Tuple}C{(phi, lam)}).
1164 @see: Property C{lamVermeille}.
1165 '''
1166 return PhiLam2Tuple(radians(self.lat), self.lamVermeille, name=self.name)
1168 phiVermeille = phi
1170 def toCartesian(self, Cartesian=None, **Cartesian_kwds):
1171 '''Return the geocentric C{(x, y, z)} coordinates as an ellipsoidal or spherical
1172 C{Cartesian}.
1174 @kwarg Cartesian: Optional class to return C{(x, y, z)} (L{ellipsoidalKarney.Cartesian},
1175 L{ellipsoidalNvector.Cartesian}, L{ellipsoidalVincenty.Cartesian},
1176 L{sphericalNvector.Cartesian} or L{sphericalTrigonometry.Cartesian})
1177 or C{None}.
1178 @kwarg Cartesian_kwds: Optionally, additional B{C{Cartesian}} keyword arguments, ignored
1179 if C{B{Cartesian} is None}.
1181 @return: A B{C{Cartesian}} instance or a L{Vector4Tuple}C{(x, y, z, h)} if C{B{Cartesian}
1182 is None}.
1184 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}} item.
1185 '''
1186 if _isin(Cartesian, None, Vector4Tuple):
1187 r = self.xyzh
1188 elif Cartesian is Vector3Tuple:
1189 r = self.xyz
1190 else:
1191 _xsubclassof(self._CartesianBase, Cartesian=Cartesian)
1192 r = Cartesian(self, **_name1__(Cartesian_kwds, _or_nameof=self))
1193 return r
1195 def toDatum(self, datum2, **name):
1196 '''Convert this C{Ecef9Tuple} to an other datum.
1198 @arg datum2: Datum to convert I{to} (L{Datum}).
1199 @kwarg name: Optional C{B{name}=NN} (C{str}).
1201 @return: The converted 9-Tuple (C{Ecef9Tuple}).
1203 @raise TypeError: The B{C{datum2}} is not a L{Datum}.
1204 '''
1205 n = _name__(name, _or_nameof=self)
1206 if _isin(self.datum, None, datum2): # PYCHOK _Names_
1207 r = self.copy(name=n)
1208 else:
1209 c = self._CartesianBase(self, datum=self.datum, name=n) # PYCHOK _Names_
1210 # c.toLatLon converts datum, x, y, z, lat, lon, etc.
1211 # and returns another Ecef9Tuple iff LatLon is None
1212 r = c.toLatLon(datum=datum2, LatLon=None)
1213 return r
1215 def toLatLon(self, LatLon=None, **LatLon_kwds):
1216 '''Return the geodetic C{(lat, lon, height[, datum])} coordinates.
1218 @kwarg LatLon: Optional class to return C{(lat, lon, height[, datum])} or C{None}.
1219 @kwarg LatLon_kwds: Optional B{C{height}}, B{C{datum}} and other B{C{LatLon}}
1220 keyword arguments.
1222 @return: A B{C{LatLon}} instance or if C{B{LatLon} is None}, a L{LatLon4Tuple}C{(lat,
1223 lon, height, datum)} or L{LatLon3Tuple}C{(lat, lon, height)} if C{datum} is
1224 specified or not.
1226 @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}} item.
1227 '''
1228 lat, lon, D = self.lat, self.lon, self.datum # PYCHOK Ecef9Tuple
1229 kwds = _name1__(LatLon_kwds, _or_nameof=self)
1230 kwds = _xkwds(kwds, height=self.height, datum=D) # PYCHOK Ecef9Tuple
1231 d = kwds.get(_datum_, LatLon)
1232 if LatLon is None:
1233 r = LatLon3Tuple(lat, lon, kwds[_height_], name=kwds[_name_])
1234 if d is not None:
1235 # assert d is not LatLon
1236 r = r.to4Tuple(d) # checks type(d)
1237 else:
1238 if d is None:
1239 _ = kwds.pop(_datum_) # remove None datum
1240 r = LatLon(lat, lon, **kwds)
1241 _xdatum(_xattr(r, datum=D), D)
1242 return r
1244 def toVector(self, Vector=None, **Vector_kwds):
1245 '''Return these geocentric C{(x, y, z)} coordinates as vector.
1247 @kwarg Vector: Optional vector class to return C{(x, y, z)} or C{None}.
1248 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword arguments,
1249 ignored if C{B{Vector} is None}.
1251 @return: A B{C{Vector}} instance or a L{Vector3Tuple}C{(x, y, z)} if
1252 C{B{Vector} is None}.
1254 @raise TypeError: Invalid B{C{Vector}} or B{C{Vector_kwds}} item.
1256 @see: Propertes C{xyz} and C{xyzh}
1257 '''
1258 return self.xyz if Vector is None else Vector(
1259 *self.xyz, **_name1__(Vector_kwds, _or_nameof=self)) # PYCHOK Ecef9Tuple
1261# def _T_x_M(self, T):
1262# '''(INTERNAL) Update M{self.M = T.multiply(self.M)}.
1263# '''
1264# return self.dup(M=T.multiply(self.M))
1266 @Property_RO
1267 def xyz(self):
1268 '''Get the geocentric C{(x, y, z)} coordinates (L{Vector3Tuple}C{(x, y, z)}).
1269 '''
1270 return Vector3Tuple(self.x, self.y, self.z, name=self.name)
1272 @Property_RO
1273 def xyzh(self):
1274 '''Get the geocentric C{(x, y, z)} coordinates and C{height} (L{Vector4Tuple}C{(x, y, z, h)})
1275 '''
1276 return self.xyz.to4Tuple(self.height)
1279def _4Ecef(this, Ecef): # in .datums.Datum.ecef, .ellipsoids.Ellipsoid.ecef
1280 '''Return an ECEF converter for C{this} L{Datum} or L{Ellipsoid}.
1281 '''
1282 if Ecef is None:
1283 Ecef = EcefKarney
1284 else:
1285 _xinstanceof(*_Ecefs, Ecef=Ecef)
1286 return Ecef(this, name=this.name)
1289def _llhn4(latlonh, lon, height, suffix=NN, Error=EcefError, **name): # in .ltp
1290 '''(INTERNAL) Get a C{(lat, lon, h, name)} 4-tuple.
1291 '''
1292 try:
1293 lat, lon = latlonh.lat, latlonh.lon
1294 h = _xattr(latlonh, height=_xattr(latlonh, h=height))
1295 n = _name__(name, _or_nameof=latlonh) # == latlonh._name__(name)
1296 except AttributeError:
1297 lat, h, n = latlonh, height, _name__(**name)
1298 try:
1299 return Lat(lat), Lon(lon), Height(h), n
1300 except (TypeError, ValueError) as x:
1301 t = _lat_, _lon_, _height_
1302 if suffix:
1303 t = (_ + suffix for _ in t)
1304 d = dict(zip(t, (lat, lon, h)))
1305 raise Error(cause=x, **d)
1308def _norm3(y, x, eps=0):
1309 '''(INTERNAL) Return C{y, x, h} normalized.
1310 '''
1311 h = hypot(y, x) # EPS0, EPS_2
1312 return (y / h, x / h, h) if h > eps else (_0_0, _1_0, h)
1315def _norm7(y, x, z=0, E=_EWGS84):
1316 '''(INTERNAL) Return C{phi, lam, p, h, C}.
1317 '''
1318 sb, cb, p = _norm3(y, x) # lam, distance to polar axis
1319 sa, ca, h = _norm3(z, p) # phi, distance to earth center
1320 if h > E._heightMax:
1321 # We are really far away (> 12M light years). Treat the earth
1322 # as a point and h above as an acceptable approximation to the
1323 # height. This avoids overflow, e.g., in the computation of d
1324 # below. It's possible that h has overflowed to INF, that's OK.
1325 # Treat finite x, y, but R overflows to +INF by scaling by 2.
1326 sb, cb, p = _norm3(y * _0_5, x * _0_5)
1327 sa, ca, _ = _norm3(z * _0_5, p)
1328 C = 4
1329 else:
1330 C = 0
1331 return sa, ca, sb, cb, p, h, C
1334def _xEcef(Ecef): # PYCHOK .latlonBase
1335 '''(INTERNAL) Validate B{C{Ecef}} I{class}.
1336 '''
1337 if issubclassof(Ecef, _EcefBase):
1338 return Ecef
1339 raise _TypesError(_Ecef_, Ecef, *_Ecefs)
1342# kwd lon00 unused but will throw a TypeError if misspelled, etc.
1343def _xyzn4(xyz, y, z, Types, Error=EcefError, lon00=0, # PYCHOK unused
1344 _xyz_y_z_names=_xyz_y_z, **name): # in .ltp
1345 '''(INTERNAL) Get an C{(x, y, z, name)} 4-tuple.
1346 '''
1347 try:
1348 n = _name__(name, _or_nameof=xyz) # == xyz._name__(name)
1349 try:
1350 t = xyz.x, xyz.y, xyz.z, n
1351 if not isinstance(xyz, Types):
1352 raise _TypesError(_xyz_y_z_names[0], xyz, *Types)
1353 except AttributeError:
1354 t = map1(float, xyz, y, z) + (n,)
1355 except (TypeError, ValueError) as x:
1356 d = dict(zip(_xyz_y_z_names, (xyz, y, z)))
1357 raise Error(cause=x, **d)
1358 return t
1359# assert _xyz_y_z == _args_kwds_names(_xyzn4)[:3]
1362_Ecefs = tuple(_ for _ in locals().values()
1363 if issubclassof(_, _EcefBase) and
1364 _ is not _EcefBase)
1365__all__ += _ALL_DOCS(_EcefBase)
1367# **) MIT License
1368#
1369# Copyright (C) 2016-2026 -- mrJean1 at Gmail -- All Rights Reserved.
1370#
1371# Permission is hereby granted, free of charge, to any person obtaining a
1372# copy of this software and associated documentation files (the "Software"),
1373# to deal in the Software without restriction, including without limitation
1374# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1375# and/or sell copies of the Software, and to permit persons to whom the
1376# Software is furnished to do so, subject to the following conditions:
1377#
1378# The above copyright notice and this permission notice shall be included
1379# in all copies or substantial portions of the Software.
1380#
1381# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1382# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1383# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1384# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1385# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1386# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1387# OTHER DEALINGS IN THE SOFTWARE.