Coverage for pygeodesy/ecef.py: 95%
447 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -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{EcefVeness} transcoded from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal,
12Cartesian<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}, class L{EcefYou}
13implementing I{Rey-Jer You}'s U{transformations <https://www.ResearchGate.net/publication/240359424>} and
14classes L{EcefFarrell22} and L{EcefFarrell22} from I{Jay A. Farrell}'s U{Table 2.1 and 2.2
15<https://Books.Google.com/books?id=fW4foWASY6wC>}, page 29-30.
17Following is a copy of I{Karney}'s U{Detailed Description
18<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Geocentric.html>}.
20Convert between geodetic coordinates C{lat}-, C{lon}gitude and height C{h} (measured vertically
21from the surface of the ellipsoid) to geocentric C{x}, C{y} and C{z} coordinates, also known as
22I{Earth-Centered, Earth-Fixed} (U{ECEF<https://WikiPedia.org/wiki/ECEF>}).
24The origin of geocentric coordinates is at the center of the earth. The C{z} axis goes thru
25the North pole, C{lat} = 90°. The C{x} axis goes thru C{lat} = 0°, C{lon} = 0°.
27The I{local (cartesian) origin} is at (C{lat0}, C{lon0}, C{height0}). The I{local} C{x} axis points
28East, the I{local} C{y} axis points North and the I{local} C{z} axis is normal to the ellipsoid. The
29plane C{z = -height0} is tangent to the ellipsoid, hence the alternate name I{local tangent plane}.
31Forward conversion from geodetic to geocentric (ECEF) coordinates is straightforward.
33For the reverse transformation we use Hugues Vermeille's U{I{Direct transformation from geocentric
34coordinates to geodetic coordinates}<https://DOI.org/10.1007/s00190-002-0273-6>}, J. Geodesy
35(2002) 76, page 451-454.
37Several changes have been made to ensure that the method returns accurate results for all finite
38inputs (even if h is infinite). The changes are described in Appendix B of C. F. F. Karney
39U{I{Geodesics on an ellipsoid of revolution}<https://ArXiv.org/abs/1102.1215v1>}, Feb. 2011, 85,
40105-117 (U{preprint<https://ArXiv.org/abs/1102.1215v1>}). Vermeille similarly updated his method
41in U{I{An analytical method to transform geocentric into geodetic coordinates}
42<https://DOI.org/10.1007/s00190-010-0419-x>}, J. Geodesy (2011) 85, page 105-117. See U{Geocentric
43coordinates<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>} for more information.
45The errors in these routines are close to round-off. Specifically, for points within 5,000 km of
46the surface of the ellipsoid (either inside or outside the ellipsoid), the error is bounded by 7
47nm (7 nanometers) for the WGS84 ellipsoid. See U{Geocentric coordinates
48<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>} for further information on the errors.
50@see: Module L{ltp} and class L{LocalCartesian}, a transcription of I{Charles Karney}'s C++ class
51U{LocalCartesian <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1LocalCartesian.html>},
52providing conversion to and from I{local} cartesian cordinates in a I{local tangent plane} as
53opposed to I{geocentric} (ECEF) ones.
54'''
56from pygeodesy.basics import copysign0, isscalar, issubclassof, neg, map1, \
57 _xinstanceof, _xsubclassof
58from pygeodesy.constants import EPS, EPS0, EPS02, EPS1, EPS2, EPS_2, PI, PI_2, \
59 _0_0, _0_0001, _0_01, _0_5, _1_0, _1_0_1T, _2_0, _3_0, _4_0, \
60 _6_0, _60_0, _90_0, _100_0, isnon0, \
61 _N_2_0 # PYCHOK used!
62from pygeodesy.datums import a_f2Tuple, _ellipsoidal_datum
63# from pygeodesy.ellipsoids import a_f2Tuple # from .datums
64from pygeodesy.errors import _IndexError, LenError, _ValueError, _TypesError, \
65 _xdatum, _xkwds
66from pygeodesy.fmath import cbrt, fdot, hypot, hypot1, hypot2_
67from pygeodesy.fsums import Fsum, fsum_
68from pygeodesy.interns import NN, _a_, _C_, _datum_, _ellipsoid_, _f_, _h_, \
69 _height_, _lat_, _lon_, _M_, _name_, _singular_, \
70 _SPACE_, _x_, _xyz_, _y_, _z_
71from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
72from pygeodesy.named import _NamedBase, _NamedTuple, notOverloaded, _Pass, _xnamed
73from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \
74 PhiLam2Tuple, Vector3Tuple, Vector4Tuple
75from pygeodesy.props import deprecated_method, Property_RO, property_RO, property_doc_
76from pygeodesy.streprs import Fmt, unstr
77from pygeodesy.units import Height, Int, Lam, Lat, Lon, Meter, Phi, Scalar, Scalar_
78from pygeodesy.utily import atan2d, degrees90, degrees180, sincos2, sincos2_, \
79 sincos2d_
81from math import asin, atan2, cos, degrees, fabs, radians, sqrt
83__all__ = _ALL_LAZY.ecef
84__version__ = '23.03.19'
86_Ecef_ = 'Ecef'
87_prolate_ = 'prolate'
88_TRIPS = 17 # 8..9 sufficient, EcefSudano.reverse
89_xyz_y_z = _xyz_, _y_, _z_ # _xargs_names(_xyzn4)[:3]
92class EcefError(_ValueError):
93 '''An ECEF or C{Ecef*} related issue.
94 '''
95 pass
98def _llhn4(latlonh, lon, height, suffix=NN, Error=EcefError, name=NN): # in .ltp.LocalCartesian.forward and -.reset
99 '''(INTERNAL) Get C{lat, lon, h, name} as C{4-tuple}.
100 '''
101 try:
102 lat, lon = latlonh.lat, latlonh.lon
103 h = getattr(latlonh, _height_,
104 getattr(latlonh, _h_, height))
105 n = getattr(latlonh, _name_, NN)
106 except AttributeError:
107 lat, h, n = latlonh, height, NN
109 try:
110 llhn = Lat(lat), Lon(lon), Height(h), (name or n)
111 except (TypeError, ValueError) as x:
112 t = _lat_, _lon_, _height_
113 if suffix:
114 t = (_ + suffix for _ in t)
115 d = dict(zip(t, (lat, lon, h)))
116 raise Error(cause=x, **d)
117 return llhn
120def _xyzn4(xyz, y, z, Types, Error=EcefError, name=NN, # in .ltp
121 _xyz_y_z_names=_xyz_y_z):
122 '''(INTERNAL) Get an C{(x, y, z, name)} 4-tuple.
123 '''
124 try:
125 try:
126 t = xyz.x, xyz.y, xyz.z, getattr(xyz, _name_, name)
127 if not isinstance(xyz, Types):
128 raise _TypesError(_xyz_y_z_names[0], xyz, *Types)
129 except AttributeError:
130 t = map1(float, xyz, y, z) + (name,)
132 except (TypeError, ValueError) as x:
133 d = dict(zip(_xyz_y_z_names, (xyz, y, z)))
134 raise Error(cause=x, **d)
135 return t
137# assert _xyz_y_z == _xargs_names(_xyzn4)[:3]
140class _EcefBase(_NamedBase):
141 '''(INTERNAL) Base class for L{EcefFarrell21}, L{EcefFarrell22}, L{EcefKarney},
142 L{EcefSudano}, L{EcefVeness} and L{EcefYou}.
143 '''
144 _datum = None
145 _E = None
147 def __init__(self, a_ellipsoid, f=None, name=NN):
148 '''New C{Ecef*} converter.
150 @arg a_ellipsoid: A (non-prolate) ellipsoid (L{Ellipsoid}, L{Ellipsoid2},
151 L{Datum} or L{a_f2Tuple}) or C{scalar} ellipsoid's
152 equatorial radius (C{meter}).
153 @kwarg f: C{None} or the ellipsoid flattening (C{scalar}), required
154 for C{scalar} B{C{a_ellipsoid}}, C{B{f}=0} represents a
155 sphere, negative B{C{f}} a prolate ellipsoid.
156 @kwarg name: Optional name (C{str}).
158 @raise EcefError: If B{C{a_ellipsoid}} not L{Ellipsoid}, L{Ellipsoid2},
159 L{Datum} or L{a_f2Tuple} or C{scalar} or B{C{f}} not
160 C{scalar} or if C{scalar} B{C{a_ellipsoid}} not positive
161 or B{C{f}} not less than 1.0.
162 '''
163 if name:
164 self.name = name
165 try:
166 E = a_ellipsoid
167 if f is None:
168 pass
169 elif isscalar(E) and isscalar(f):
170 E = a_f2Tuple(E, f)
171 else:
172 raise ValueError # _invalid_
174 d = _ellipsoidal_datum(E, name=name)
175 E = d.ellipsoid
176 if E.a < EPS or E.f > EPS1:
177 raise ValueError # _invalid_
179 except (TypeError, ValueError) as x:
180 t = unstr(self.classname, a=a_ellipsoid, f=f)
181 raise EcefError(_SPACE_(t, _ellipsoid_), cause=x)
183 self._datum = d
184 self._E = E
186 def __eq__(self, other):
187 '''Compare this and an other Ecef.
189 @arg other: The other ecef (C{Ecef*}).
191 @return: C{True} if equal, C{False} otherwise.
192 '''
193 return other is self or (isinstance(other, self.__class__) and
194 other.ellipsoid == self.ellipsoid)
196 @Property_RO
197 def equatoradius(self):
198 '''Get the I{equatorial} radius, semi-axis (C{meter}).
199 '''
200 return self.ellipsoid.a
202 a = equatorialRadius = equatoradius # Karney property
204 @Property_RO
205 def datum(self):
206 '''Get the datum (L{Datum}).
207 '''
208 return self._datum
210 @Property_RO
211 def ellipsoid(self):
212 '''Get the ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).
213 '''
214 return self._E
216 @Property_RO
217 def flattening(self): # Karney property
218 '''Get the I{flattening} (C{float}), M{(a - b) / a}, positive for
219 I{oblate}, negative for I{prolate} or C{0} for I{near-spherical}.
220 '''
221 return self.ellipsoid.f
223 f = flattening
225 def _forward(self, lat, lon, h, name, M=False, _philam=False): # in .ltp.LocalCartesian.forward and -.reset
226 '''(INTERNAL) Common for all C{Ecef*}.
227 '''
228 E = self.ellipsoid
230 if _philam:
231 sa, ca, sb, cb = sincos2_(lat, lon)
232 lat = Lat(degrees90( lat))
233 lon = Lon(degrees180(lon))
234 else:
235 sa, ca, sb, cb = sincos2d_(lat, lon)
237 n = E.roc1_(sa, ca) if self._isYou else E.roc1_(sa)
238 z = (h + n * E.e21) * sa
239 x = (h + n) * ca
241 m = self._Matrix(sa, ca, sb, cb) if M else None
242 return Ecef9Tuple(x * cb, x * sb, z, lat, lon, h,
243 0, m, self.datum,
244 name=name or self.name)
246 def forward(self, latlonh, lon=None, height=0, M=False, name=NN):
247 '''Convert from geodetic C{(lat, lon, height)} to geocentric C{(x, y, z)}.
249 @arg latlonh: Either a C{LatLon}, an L{Ecef9Tuple} or C{scalar}
250 latitude (C{degrees}).
251 @kwarg lon: Optional C{scalar} longitude for C{scalar} B{C{latlonh}}
252 (C{degrees}).
253 @kwarg height: Optional height (C{meter}), vertically above (or below)
254 the surface of the ellipsoid.
255 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
256 @kwarg name: Optional name (C{str}).
258 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
259 geocentric C{(x, y, z)} coordinates for the given geodetic ones
260 C{(lat, lon, height)}, case C{C} 0, optional C{M} (L{EcefMatrix})
261 and C{datum} if available.
263 @raise EcefError: If B{C{latlonh}} not C{LatLon}, L{Ecef9Tuple} or
264 C{scalar} or B{C{lon}} not C{scalar} for C{scalar}
265 B{C{latlonh}} or C{abs(lat)} exceeds 90°.
267 @note: Use method C{.forward_} to specify C{lat} and C{lon} in C{radians}
268 and avoid double angle conversions.
269 '''
270 llhn = _llhn4(latlonh, lon, height, name=name)
271 return _EcefBase._forward(self, *llhn, M=M)
273 def forward_(self, phi, lam, height=0, M=False, name=NN):
274 '''Like method C{.forward} except with geodetic lat- and longitude given
275 in I{radians}.
277 @arg phi: Latitude in I{radians} (C{scalar}).
278 @arg lam: Longitude in I{radians} (C{scalar}).
279 @kwarg height: Optional height (C{meter}), vertically above (or below)
280 the surface of the ellipsoid.
281 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
282 @kwarg name: Optional name (C{str}).
284 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
285 with C{lat} set to C{degrees90(B{phi})} and C{lon} to
286 C{degrees180(B{lam})}.
288 @raise EcefError: If B{C{phi}} or B{C{lam}} invalid or not C{scalar}.
289 '''
290 try: # like function C{_llhn4} above
291 plhn = Phi(phi), Lam(lam), Height(height), name
292 except (TypeError, ValueError) as x:
293 raise EcefError(phi=phi, lam=lam, height=height, cause=x)
294 return self._forward(*plhn, M=M, _philam=True)
296 @property_RO
297 def _Geocentrics(self):
298 '''(INTERNAL) Valid geocentric classes.
299 '''
300 t = Ecef9Tuple, _MODS.cartesianBase.CartesianBase
301 _EcefBase._Geocentrics = t # overwrite the property
302 return t
304 @Property_RO
305 def _isYou(self):
306 '''(INTERNAL) Is this an C{EcefYou}?.
307 '''
308 return isinstance(self, EcefYou)
310 def _Matrix(self, sa, ca, sb, cb):
311 '''Creation a rotation matrix.
313 @arg sa: C{sin(phi)} (C{float}).
314 @arg ca: C{cos(phi)} (C{float}).
315 @arg sb: C{sin(lambda)} (C{float}).
316 @arg cb: C{cos(lambda)} (C{float}).
318 @return: An L{EcefMatrix}.
319 '''
320 return self._xnamed(EcefMatrix(sa, ca, sb, cb))
322 def reverse(self, xyz, y=None, z=None, M=False, name=NN): # PYCHOK no cover
323 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
324 '''
325 notOverloaded(self, xyz, y=y, z=z, M=M, name=name)
327 def toStr(self, prec=9, **unused): # PYCHOK signature
328 '''Return this C{Ecef*} as a string.
330 @kwarg prec: Precision, number of decimal digits (0..9).
332 @return: This C{Ecef*} (C{str}).
333 '''
334 return self.attrs(_a_, _f_, _datum_, _name_, prec=prec) # _ellipsoid_
337class EcefFarrell21(_EcefBase):
338 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)
339 coordinates based on I{Jay A. Farrell}'s U{Table 2.1<https://Books.Google.com/
340 books?id=fW4foWASY6wC>}, page 29.
341 '''
343 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M
344 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using
345 I{Farrell}'s U{Table 2.1<https://Books.Google.com/books?id=fW4foWASY6wC>},
346 page 29.
348 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
349 coordinate (C{meter}).
350 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
351 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
352 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
353 @kwarg name: Optional name (C{str}).
355 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
356 geodetic coordinates C{(lat, lon, height)} for the given geocentric
357 ones C{(x, y, z)}, case C{C=1}, C{M=None} always and C{datum}
358 if available.
360 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
361 not C{scalar} for C{scalar} B{C{xyz}} or C{sqrt} domain or
362 zero division error.
364 @see: L{EcefFarrell22} and L{EcefVeness}.
365 '''
366 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name)
368 E = self.ellipsoid
369 a = E.a
370 a2 = E.a2
371 b2 = E.b2
372 e_ = E.a_b * E.e # 0.0820944... WGS84
373 e2 = E.e2
374 e4 = E.e4
376 try: # names as page 29
377 z2 = z**2
378 ez = (_1_0 - e2) * z2 # E.e2s2(z)
380 p = hypot(x, y)
381 p2 = p**2
382 G = p2 + ez - e2 * (a2 - b2) # p2 + ez - e4 * a2
383 F = b2 * z2 * 54
384 c = e4 * p2 * F / G**3
385 s = cbrt(_1_0 + c + sqrt(c**2 + c * 2))
386 P = F / (_3_0 * (fsum_(_1_0, s, _1_0 / s) * G)**2)
387 Q = sqrt(_1_0 + _2_0 * e4 * P)
388 Q1 = Q + _1_0
389 r0 = P * e2 * p / Q1 - sqrt(fsum_(a2 * (Q1 / Q) * _0_5,
390 -P * ez / (Q * Q1),
391 -P * p2 * _0_5))
392 r = p + e2 * r0
393 v = b2 / (a * sqrt(r**2 + ez))
395 h = hypot(r, z) * (_1_0 - v)
396 t = atan2(z + e_**2 * v * z, p)
397 # note, phi and lam are swapped on page 29
399 except (ValueError, ZeroDivisionError) as e:
400 raise EcefError(x=x, y=y, z=z, cause=e)
402 return Ecef9Tuple(x, y, z, degrees90(t), atan2d(y, x), h,
403 1, None, self.datum,
404 name=name or self.name)
407class EcefFarrell22(_EcefBase):
408 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)
409 coordinates based on I{Jay A. Farrell}'s U{Table 2.2<https://Books.Google.com/
410 books?id=fW4foWASY6wC>}, page 30.
411 '''
413 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M
414 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using
415 I{Farrell}'s U{Table 2.2<https://Books.Google.com/books?id=fW4foWASY6wC>},
416 page 30.
418 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
419 coordinate (C{meter}).
420 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
421 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
422 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
423 @kwarg name: Optional name (C{str}).
425 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
426 geodetic coordinates C{(lat, lon, height)} for the given geocentric
427 ones C{(x, y, z)}, case C{C=1}, C{M=None} always and C{datum}
428 if available.
430 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
431 not C{scalar} for C{scalar} B{C{xyz}} or C{sqrt} domain or
432 zero division error.
434 @see: L{EcefFarrell21} and L{EcefVeness}.
435 '''
436 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name)
438 E = self.ellipsoid
439 a = E.a
440 b = E.b
442 try: # see EcefVeness.reverse
443 p = hypot(x, y)
444 s, c = sincos2(atan2(z * a, p * b))
446 t = atan2(z + E.e22 * b * s**3,
447 p - E.e2 * a * c**3)
449 s, c = sincos2(t)
450 if c: # E.roc1_(s) = E.a / sqrt(1 - E.e2 * s**2)
451 h = p / c - (E.roc1_(s) if s else a)
452 else: # polar
453 h = fabs(z) - b
454 # note, phi and lam are swapped on page 30
456 except (ValueError, ZeroDivisionError) as e:
457 raise EcefError(x=x, y=y, z=z, cause=e)
459 return Ecef9Tuple(x, y, z, degrees90(t), atan2d(y, x), h,
460 1, None, self.datum,
461 name=name or self.name)
464class EcefKarney(_EcefBase):
465 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)
466 coordinates transcoded from I{Karney}'s C++ U{Geocentric<https://GeographicLib.SourceForge.io/
467 C++/doc/classGeographicLib_1_1Geocentric.html>} methods.
469 @note: On methods C{.forward} and C{.forwar_}, let C{v} be a unit vector located
470 at C{(lat, lon, h)}. We can express C{v} as column vectors in one of two
471 ways, C{v1} in east, north, up coordinates (where the components are
472 relative to a local coordinate system at C{C(lat0, lon0, h0)}) or as C{v0}
473 in geocentric C{x, y, z} coordinates. Then, M{v0 = M ⋅ v1} where C{M} is
474 the rotation matrix.
475 '''
477 @Property_RO
478 def hmax(self):
479 '''Get the distance or height limit (C{meter}, conventionally).
480 '''
481 return self.equatoradius / EPS_2 # self.equatoradius * _2_EPS, 12M lighyears
483 def reverse(self, xyz, y=None, z=None, M=False, name=NN):
484 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}.
486 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
487 coordinate (C{meter}).
488 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
489 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
490 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).
491 @kwarg name: Optional name (C{str}).
493 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
494 geodetic coordinates C{(lat, lon, height)} for the given geocentric
495 ones C{(x, y, z)}, case C{C}, optional C{M} (L{EcefMatrix}) and
496 C{datum} if available.
498 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
499 not C{scalar} for C{scalar} B{C{xyz}}.
501 @note: In general, there are multiple solutions and the result which minimizes
502 C{height} is returned, i.e., C{(lat, lon)} corresponds to the closest
503 point on the ellipsoid. If there are still multiple solutions with
504 different latitudes (applies only if C{z} = 0), then the solution with
505 C{lat} > 0 is returned. If there are still multiple solutions with
506 different longitudes (applies only if C{x} = C{y} = 0) then C{lon} = 0
507 is returned. The returned C{height} value is not below M{−E.a * (1 −
508 E.e2) / sqrt(1 − E.e2 * sin(lat)**2)}. The returned C{lon} is in the
509 range [−180°, 180°]. Like C{forward} above, M{v1 = Transpose(M) ⋅ v0}.
510 '''
511 def norm3(y, x):
512 h = hypot(y, x) # EPS0, EPS_2
513 return (y / h, x / h, h) if h > 0 else (_0_0, _1_0, h)
515 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name)
517 E = self.ellipsoid
519 sb, cb, R = norm3(y, x)
520 h = hypot(R, z) # distance to earth center
521 if h > self.hmax: # PYCHOK no cover
522 # We are really far away (> 12M light years). Treat the earth
523 # as a point and h, above as an acceptable approximation to the
524 # height. This avoids overflow, e.g., in the computation of disc
525 # below. It's possible that h has overflowed to INF, that's OK.
526 # Treat finite x, y, but R overflows to +INF by scaling by 2.
527 sb, cb, R = norm3(y * _0_5, x * _0_5)
528 sa, ca, _ = norm3(z * _0_5, R)
529 C = 1
531 elif E.e4: # E.isEllipsoidal
532 # Treat prolate spheroids by swapping R and Z here and by
533 # switching the arguments to phi = atan2(...) at the end.
534 p = (R / E.a)**2
535 q = E.e21 * (z / E.a)**2
536 if E.isProlate:
537 p, q = q, p
538 r = p + q - E.e4
539 e = E.e4 * q
540 if e or r > 0:
541 # Avoid possible division by zero when r = 0 by multiplying
542 # equations for s and t by r^3 and r, respectively.
543 s = e * p / _4_0 # s = r^3 * s
544 u = r = r / _6_0
545 r2 = r**2
546 r3 = r * r2
547 t3 = s + r3
548 disc = s * (r3 + t3)
549 if disc < 0:
550 # t is complex, but the way u is defined, the result is real.
551 # There are three possible cube roots. We choose the root
552 # which avoids cancellation. Note, disc < 0 implies r < 0.
553 u += _2_0 * r * cos(atan2(sqrt(-disc), -t3) / _3_0)
554 else:
555 # Pick the sign on the sqrt to maximize abs(T3). This
556 # minimizes loss of precision due to cancellation. The
557 # result is unchanged because of the way the t is used
558 # in definition of u.
559 if disc > 0:
560 t3 += copysign0(sqrt(disc), t3) # t3 = (r * t)^3
561 # N.B. cbrt always returns the real root, cbrt(-8) = -2.
562 t = cbrt(t3) # t = r * t
563 # t can be zero; but then r2 / t -> 0.
564 if t:
565 u = fsum_(u, t, r2 / t)
566 v = sqrt(e + u**2) # guaranteed positive
567 # Avoid loss of accuracy when u < 0. Underflow doesn't occur in
568 # E.e4 * q / (v - u) because u ~ e^4 when q is small and u < 0.
569 uv = (e / (v - u)) if u < 0 else (u + v) # u+v, guaranteed positive
570 # Need to guard against w going negative due to roundoff in uv - q.
571 w = max(_0_0, E.e2abs * (uv - q) / (_2_0 * v))
572 # Rearrange expression for k to avoid loss of accuracy due to
573 # subtraction. Division by 0 not possible because uv > 0, w >= 0.
574 k1 = k2 = uv / (sqrt(uv + w**2) + w)
575 if E.isProlate:
576 k1 -= E.e2
577 else:
578 k2 += E.e2
579 sa, ca, h = norm3(z / k1, R / k2)
580 h *= k1 - E.e21
581 C = 2
583 else: # e = E.e4 * q == 0 and r <= 0
584 # This leads to k = 0 (oblate, equatorial plane) and k + E.e^2 = 0
585 # (prolate, rotation axis) and the generation of 0/0 in the general
586 # formulas for phi and h, using the general formula and division
587 # by 0 in formula for h. Handle this case by taking the limits:
588 # f > 0: z -> 0, k -> E.e2 * sqrt(q) / sqrt(E.e4 - p)
589 # f < 0: r -> 0, k + E.e2 -> -E.e2 * sqrt(q) / sqrt(E.e4 - p)
590 q = E.e4 - p
591 if E.isProlate:
592 p, q = q, p
593 e = E.a
594 else:
595 e = E.b2_a
596 sa, ca, h = norm3(sqrt(q * E._1_e21), sqrt(p))
597 if z < 0:
598 sa = neg(sa) # for tiny negative z, not for prolate
599 h *= neg(e / E.e2abs)
600 C = 3
602 else: # E.e4 == 0, spherical case
603 # Dealing with underflow in the general case with E.e2 = 0 is
604 # difficult. Origin maps to North pole, same as with ellipsoid.
605 sa, ca, _ = norm3((z if h else _1_0), R)
606 h -= E.a
607 C = 4
609 m = self._Matrix(sa, ca, sb, cb) if M else None
610 return Ecef9Tuple(x, y, z, atan2d(sa, ca),
611 atan2d(sb, cb), h,
612 C, m, self.datum,
613 name=name or self.name)
616class EcefSudano(_EcefBase):
617 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
618 based on I{John J. Sudano}'s U{paper<https://www.ResearchGate.net/publication/
619 3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}.
620 '''
621 _tol = EPS2
623 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M
624 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using
625 I{Sudano}'s U{iterative method<https://www.ResearchGate.net/publication/
626 3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}.
628 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
629 coordinate (C{meter}).
630 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
631 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
632 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
633 @kwarg name: Optional name (C{str}).
635 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with geodetic
636 coordinates C{(lat, lon, height)} for the given geocentric ones C{(x, y, z)},
637 iteration C{C}, C{M=None} always and C{datum} if available.
639 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
640 not C{scalar} for C{scalar} B{C{xyz}} or no convergence.
641 '''
642 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name)
644 E = self.ellipsoid
645 e = E.e2 * E.a
646 h = hypot(x, y) # Rh
647 d = e - h
649 a = atan2(z, h * E.e21)
650 sa, ca = sincos2(fabs(a))
651 # Sudano's Eq (A-6) and (A-7) refactored/reduced,
652 # replacing Rn from Eq (A-4) with n = E.a / ca:
653 # N = ca**2 * ((z + E.e2 * n * sa) * ca - h * sa)
654 # = ca**2 * (z * ca + E.e2 * E.a * sa - h * sa)
655 # = ca**2 * (z * ca + (E.e2 * E.a - h) * sa)
656 # D = ca**3 * (E.e2 * n / E.e2s2(sa)) - h
657 # = ca**2 * (E.e2 * E.a / E.e2s2(sa) - h / ca**2)
658 # N / D = (z * ca + (E.e2 * E.a - h) * sa) /
659 # (E.e2 * E.a / E.e2s2(sa) - h / ca**2)
660 tol = self.tolerance
661 _S2_ = Fsum(sa).fsum2_
662 for C in range(1, _TRIPS):
663 ca2 = _1_0 - sa**2
664 if ca2 < EPS_2: # PYCHOK no cover
665 ca = _0_0
666 break
667 ca = sqrt(ca2)
668 r = e / E.e2s2(sa) - h / ca2
669 if fabs(r) < EPS_2:
670 break
671 a = None
672 sa, r = _S2_(-z * ca / r, -d * sa / r)
673 if fabs(r) < tol:
674 break
675 else:
676 t = unstr(self.reverse, x=x, y=y, z=z)
677 raise EcefError(Fmt.no_convergence(r, tol), txt=t)
679 if a is None:
680 a = copysign0(asin(sa), z)
681 h = fsum_(h * ca, fabs(z * sa), -E.a * E.e2s(sa)) # use Veness',
682 # since Sudano's Eq (7) doesn't provide the correct height
683 # h = (fabs(z) + h - E.a * cos(a + E.e21) * sa / ca) / (ca + sa)
685 r = Ecef9Tuple(x, y, z, degrees90(a), atan2d(y, x), h,
686 C, None, self.datum,
687 name=name or self.name)
688 r._iteration = C
689 return r
691 @property_doc_(''' the convergence tolerance (C{float}).''')
692 def tolerance(self):
693 '''Get the convergence tolerance (C{scalar}).
694 '''
695 return self._tol
697 @tolerance.setter # PYCHOK setter!
698 def tolerance(self, tol):
699 '''Set the convergence tolerance (C{scalar}).
701 @raise EcefError: Non-scalar or invalid B{C{tol}}.
702 '''
703 self._tol = Scalar_(tolerance=tol, low=EPS, Error=EcefError)
706class EcefVeness(_EcefBase):
707 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
708 transcoded from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal, Cartesian<https://
709 www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}.
711 @see: U{I{A Guide to Coordinate Systems in Great Britain}<https://
712 www.OrdnanceSurvey.co.UK/documents/resources/guide-coordinate-systems-great-britain.pdf>},
713 section I{B) Converting between 3D Cartesian and ellipsoidal
714 latitude, longitude and height coordinates}.
715 '''
717 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M
718 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}
719 transcoded from I{Chris Veness}' U{JavaScript<https://www.Movable-Type.co.UK/
720 scripts/geodesy/docs/latlon-ellipsoidal.js.html>}.
722 Uses B. R. Bowring’s formulation for μm precision in concise form U{I{The accuracy
723 of geodetic latitude and height equations}<https://www.ResearchGate.net/publication/
724 233668213_The_Accuracy_of_Geodetic_Latitude_and_Height_Equations>}, Survey Review,
725 Vol 28, 218, Oct 1985.
727 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
728 coordinate (C{meter}).
729 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
730 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
731 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
732 @kwarg name: Optional name (C{str}).
734 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
735 geodetic coordinates C{(lat, lon, height)} for the given geocentric
736 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum} if available.
738 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}
739 not C{scalar} for C{scalar} B{C{xyz}}.
741 @see: Toms, Ralph M. U{I{An Efficient Algorithm for Geocentric to Geodetic
742 Coordinate Conversion}<https://www.OSTI.gov/scitech/biblio/110235>},
743 Sept 1995 and U{I{An Improved Algorithm for Geocentric to Geodetic
744 Coordinate Conversion}<https://www.OSTI.gov/scitech/servlets/purl/231228>},
745 Apr 1996, both from Lawrence Livermore National Laboratory (LLNL) and
746 Sudano, John J, U{I{An exact conversion from an Earth-centered coordinate
747 system to latitude longitude and altitude}<https://www.ResearchGate.net/
748 publication/3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}.
749 '''
750 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name)
752 E = self.ellipsoid
754 p = hypot(x, y) # distance from minor axis
755 r = hypot(p, z) # polar radius
756 if min(p, r) > EPS0:
757 # parametric latitude (Bowring eqn 17, replaced)
758 t = (E.b * z) / (E.a * p) * (_1_0 + E.e22 * E.b / r)
759 c = _1_0 / hypot1(t)
760 s = t * c
762 # geodetic latitude (Bowring eqn 18)
763 a = atan2(z + E.e22 * E.b * s**3,
764 p - E.e2 * E.a * c**3)
766 # height above ellipsoid (Bowring eqn 7)
767 sa, ca = sincos2(a)
768# r = E.a / E.e2s(sa) # length of normal terminated by minor axis
769# h = p * ca + z * sa - (E.a * E.a / r)
770 h = fsum_(p * ca, z * sa, -E.a * E.e2s(sa))
772 C, lat, lon = 1, degrees90(a), atan2d(y, x)
774 # see <https://GIS.StackExchange.com/questions/28446>
775 elif p > EPS: # lat arbitrarily zero
776 C, lat, lon, h = 2, _0_0, atan2d(y, x), p - E.a
778 else: # polar lat, lon arbitrarily zero
779 C, lat, lon, h = 3, copysign0(_90_0, z), _0_0, fabs(z) - E.b
781 return Ecef9Tuple(x, y, z, lat, lon, h,
782 C, None, # M=None
783 self.datum, name=name or self.name)
786class EcefYou(_EcefBase):
787 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates
788 using I{Rey-Jer You}'s U{transformation<https://www.ResearchGate.net/publication/240359424>}.
790 @see: Featherstone, W.E., Claessens, S.J. U{I{Closed-form transformation between geodetic and
791 ellipsoidal coordinates}<https://Espace.Curtin.edu.AU/bitstream/handle/20.500.11937/
792 11589/115114_9021_geod2ellip_final.pdf>} Studia Geophysica et Geodaetica, 2008, 52,
793 pages 1-18 and U{PyMap3D <https://PyPI.org/project/pymap3d>}.
794 '''
796 def __init__(self, a_ellipsoid, f=None, name=NN):
797 _EcefBase.__init__(self, a_ellipsoid, f=f, name=name) # inherited documentation
798 E = self.ellipsoid
799 if E.isProlate or (E.a2 - E.b2) < 0:
800 raise EcefError(ellipsoid=E, txt=_prolate_)
802 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M
803 '''Convert geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}
804 using I{Rey-Jer You}'s transformation.
806 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}
807 coordinate (C{meter}).
808 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).
809 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).
810 @kwarg M: I{Ignored}, rotation matrix C{M} not available.
811 @kwarg name: Optional name (C{str}).
813 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with
814 geodetic coordinates C{(lat, lon, height)} for the given geocentric
815 ones C{(x, y, z)}, case C{C=1}, C{M=None} always and C{datum} if
816 available.
818 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or
819 B{C{z}} not C{scalar} for C{scalar} B{C{xyz}}.
820 '''
821 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name)
823 r2 = hypot2_(x, y, z)
825 E = self.ellipsoid
826 e2 = E.a2 - E.b2 # == E.e2 * E.a2
827 if e2 < 0:
828 raise EcefError(ellipsoid=E, txt=_prolate_)
829 e = sqrt(e2) # XXX sqrt0(e2)?
831 q = hypot(x, y)
832 u = fsum_(r2, -e2, hypot(r2 - e2, 2 * e * z)) * _0_5
833 if u > EPS02:
834 u = sqrt(u)
835 p = hypot(u, e)
836 B = atan2(p * z, u * q) # beta0 = atan(p / u * z / q)
837 sB, cB = sincos2(B)
838 if cB and sB:
839 p *= E.a
840 d = (p / cB - e2 * cB) / sB
841 if isnon0(d):
842 B += fsum_(u * E.b, -p, e2) / d
843 sB, cB = sincos2(B)
844 elif u < 0:
845 raise EcefError(x=x, y=y, z=z, txt=_singular_)
846 else:
847 sB, cB = copysign0(_1_0, z), _0_0
849 h = hypot(z - E.b * sB, q - E.a * cB)
850 if hypot2_(x, y, z * E.a_b) < E.a2:
851 h = neg(h) # inside ellipsoid
853 return Ecef9Tuple(x, y, z, atan2d(E.a * sB, E.b * cB), # atan(E.a_b * tan(B))
854 atan2d(y, x), h,
855 1, None, # C=1, M=None
856 self.datum, name=name or self.name)
859class EcefMatrix(_NamedTuple):
860 '''A rotation matrix.
861 '''
862 _Names_ = ('_0_0_', '_0_1_', '_0_2_', # row-order
863 '_1_0_', '_1_1_', '_1_2_',
864 '_2_0_', '_2_1_', '_2_2_')
865 _Units_ = (Scalar,) * len(_Names_)
867 def _validate(self, **_OK): # PYCHOK unused
868 '''(INTERNAL) Allow C{_Names_} with leading underscore.
869 '''
870 _NamedTuple._validate(self, _OK=True)
872 def __new__(cls, sa, ca, sb, cb, *_more):
873 '''New L{EcefMatrix} matrix.
875 @arg sa: C{sin(phi)} (C{float}).
876 @arg ca: C{cos(phi)} (C{float}).
877 @arg sb: C{sin(lambda)} (C{float}).
878 @arg cb: C{cos(lambda)} (C{float}).
879 @arg _more: (INTERNAL) from C{.multiply}.
881 @raise EcefError: If B{C{sa}}, B{C{ca}}, B{C{sb}} or
882 B{C{cb}} outside M{[-1.0, +1.0]}.
883 '''
884 t = sa, ca, sb, cb
885 if _more: # all 9 matrix elements ...
886 t += _more # ... from .multiply
888 elif max(map(fabs, t)) > _1_0:
889 raise EcefError(unstr(EcefMatrix.__name__, *t))
891 else: # build matrix from the following quaternion operations
892 # qrot(lam, [0,0,1]) * qrot(phi, [0,-1,0]) * [1,1,1,1]/2
893 # or
894 # qrot(pi/2 + lam, [0,0,1]) * qrot(-pi/2 + phi, [-1,0,0])
895 # where
896 # qrot(t,v) = [cos(t/2), sin(t/2)*v[1], sin(t/2)*v[2], sin(t/2)*v[3]]
898 # Local X axis (east) in geocentric coords
899 # M[0] = -slam; M[3] = clam; M[6] = 0;
900 # Local Y axis (north) in geocentric coords
901 # M[1] = -clam * sphi; M[4] = -slam * sphi; M[7] = cphi;
902 # Local Z axis (up) in geocentric coords
903 # M[2] = clam * cphi; M[5] = slam * cphi; M[8] = sphi;
904 t = (-sb, -cb * sa, cb * ca,
905 cb, -sb * sa, sb * ca,
906 _0_0, ca, sa)
908 return _NamedTuple.__new__(cls, *t)
910 def column(self, column):
911 '''Get matrix B{C{column}} as 3-tuple.
912 '''
913 if 0 <= column < 3:
914 return self[column::3]
915 raise _IndexError(column=column)
917 @Property_RO
918 def _column_0(self):
919 return self.column(0)
921 @Property_RO
922 def _column_1(self):
923 return self.column(1)
925 @Property_RO
926 def _column_2(self):
927 return self.column(2)
929 def copy(self, **unused): # PYCHOK signature
930 '''Make a shallow or deep copy of this instance.
932 @return: The copy (C{This class} or subclass thereof).
933 '''
934 return self.classof(*self)
936 __copy__ = __deepcopy__ = copy
938 def multiply(self, other):
939 '''Matrix multiply M{M0' ⋅ M} this matrix transposed with
940 an other matrix.
942 @arg other: The other matrix (L{EcefMatrix}).
944 @return: The matrix product (L{EcefMatrix}).
946 @raise TypeError: If B{C{other}} is not L{EcefMatrix}.
947 '''
948 _xinstanceof(EcefMatrix, other=other)
950 # like LocalCartesian.MatrixMultiply, transposed(self) X other
951 # <https://GeographicLib.SourceForge.io/C++/doc/LocalCartesian_8cpp_source.html>
952 M = (fdot(self[r::3], *other[c::3]) for r in range(3) for c in range(3))
953 return _xnamed(EcefMatrix(*M), EcefMatrix.multiply.__name__)
955 def rotate(self, xyz, *xyz0):
956 '''Forward rotation M{M0' ⋅ ([x, y, z] - [x0, y0, z0])'}.
958 @arg xyz: Local C{(x, y, z)} coordinates (C{3-tuple}).
959 @arg xyz0: Optional, local C{(x0, y0, z0)} origin (C{3-tuple}).
961 @return: Rotated C{(x, y, z)} location (C{3-tuple}).
963 @raise LenError: Unequal C{len(B{xyz})} and C{len(B{xyz0})}.
964 '''
965 if xyz0:
966 if len(xyz0) != len(xyz):
967 raise LenError(self.rotate, xyz0=len(xyz0), xyz=len(xyz))
969 xyz = tuple(c - c0 for c, c0 in zip(xyz, xyz0))
971 # x' = M[0] * x + M[3] * y + M[6] * z
972 # y' = M[1] * x + M[4] * y + M[7] * z
973 # z' = M[2] * x + M[5] * y + M[8] * z
974 return (fdot(xyz, *self._column_0),
975 fdot(xyz, *self._column_1),
976 fdot(xyz, *self._column_2))
978 def row(self, row):
979 '''Get matrix B{C{row}} as 3-tuple.
980 '''
981 if 0 <= row < 3:
982 r = row * 3
983 return self[r:r+3]
984 raise _IndexError(row=row)
986 @Property_RO
987 def _row_0(self):
988 return self.row(0)
990 @Property_RO
991 def _row_1(self):
992 return self.row(1)
994 @Property_RO
995 def _row_2(self):
996 return self.row(2)
998 def unrotate(self, xyz, *xyz0):
999 '''Inverse rotation M{[x0, y0, z0] + M0 ⋅ [x,y,z]'}.
1001 @arg xyz: Local C{(x, y, z)} coordinates (C{3-tuple}).
1002 @arg xyz0: Optional, local C{(x0, y0, z0)} origin (C{3-tuple}).
1004 @return: Unrotated C{(x, y, z)} location (C{3-tuple}).
1006 @raise LenError: Unequal C{len(B{xyz})} and C{len(B{xyz0})}.
1007 '''
1008 if xyz0:
1009 if len(xyz0) != len(xyz):
1010 raise LenError(self.unrotate, xyz0=len(xyz0), xyz=len(xyz))
1012 _xyz = _1_0_1T + xyz
1013 # x' = x0 + M[0] * x + M[1] * y + M[2] * z
1014 # y' = y0 + M[3] * x + M[4] * y + M[5] * z
1015 # z' = z0 + M[6] * x + M[7] * y + M[8] * z
1016 xyz_ = (fdot(_xyz, xyz0[0], *self._row_0),
1017 fdot(_xyz, xyz0[1], *self._row_1),
1018 fdot(_xyz, xyz0[2], *self._row_2))
1019 else:
1020 # x' = M[0] * x + M[1] * y + M[2] * z
1021 # y' = M[3] * x + M[4] * y + M[5] * z
1022 # z' = M[6] * x + M[7] * y + M[8] * z
1023 xyz_ = (fdot(xyz, *self._row_0),
1024 fdot(xyz, *self._row_1),
1025 fdot(xyz, *self._row_2))
1026 return xyz_
1029class Ecef9Tuple(_NamedTuple):
1030 '''9-Tuple C{(x, y, z, lat, lon, height, C, M, datum)} with I{geocentric}
1031 C{x}, C{y} and C{z} plus I{geodetic} C{lat}, C{lon} and C{height}, case
1032 C{C} (see the C{Ecef*.reverse} methods) and optionally, the rotation
1033 matrix C{M} (L{EcefMatrix}) and C{datum}, with C{lat} and C{lon} in
1034 C{degrees} and C{x}, C{y}, C{z} and C{height} in C{meter}, conventionally.
1035 '''
1036 _Names_ = (_x_, _y_, _z_, _lat_, _lon_, _height_, _C_, _M_, _datum_)
1037 _Units_ = ( Meter, Meter, Meter, Lat, Lon, Height, Int, _Pass, _Pass)
1039 _IndexM = _Names_.index(_M_) # for ._M_x_M
1041 @property_RO
1042 def _CartesianBase(self):
1043 '''(INTERNAL) Get/cache class C{CartesianBase}.
1044 '''
1045 Ecef9Tuple._CartesianBase = C = _MODS.cartesianBase.CartesianBase # overwrite property
1046 return C
1048 @deprecated_method
1049 def convertDatum(self, datum2): # for backward compatibility
1050 '''DEPRECATED, use method L{toDatum}.'''
1051 return self.toDatum(datum2)
1053 @Property_RO
1054 def lam(self):
1055 '''Get the longitude in C{radians} (C{float}).
1056 '''
1057 return self.philam.lam
1059 @Property_RO
1060 def lamVermeille(self):
1061 '''Get the longitude in C{radians [-PI*3/2..+PI*3/2]} after U{Vermeille
1062 <https://Search.ProQuest.com/docview/639493848>} (2004), page 95.
1064 @see: U{Karney<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>},
1065 U{Vermeille<https://Search.ProQuest.com/docview/847292978>} 2011, pp 112-113, 116
1066 and U{Featherstone, et.al.<https://Search.ProQuest.com/docview/872827242>}, page 7.
1067 '''
1068 x, y = self.x, self.y
1069 if y > EPS0:
1070 r = _N_2_0 * atan2(x, hypot(y, x) + y) + PI_2
1071 elif y < -EPS0:
1072 r = _2_0 * atan2(x, hypot(y, x) - y) - PI_2
1073 else: # y == 0
1074 r = PI if x < 0 else _0_0
1075 return Lam(Vermeille=r)
1077 @Property_RO
1078 def latlon(self):
1079 '''Get the lat-, longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).
1080 '''
1081 return LatLon2Tuple(self.lat, self.lon, name=self.name)
1083 @Property_RO
1084 def latlonheight(self):
1085 '''Get the lat-, longitude in C{degrees} and height (L{LatLon3Tuple}C{(lat, lon, height)}).
1086 '''
1087 return self.latlon.to3Tuple(self.height)
1089 @Property_RO
1090 def latlonheightdatum(self):
1091 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}).
1092 '''
1093 return self.latlonheight.to4Tuple(self.datum)
1095 @Property_RO
1096 def latlonVermeille(self):
1097 '''Get the latitude and I{Vermeille} longitude in C{degrees [-225..+225]} (L{LatLon2Tuple}C{(lat, lon)}).
1099 @see: Property C{lonVermeille}.
1100 '''
1101 return LatLon2Tuple(self.lat, self.lonVermeille, name=self.name)
1103 @Property_RO
1104 def lonVermeille(self):
1105 '''Get the longitude in C{degrees [-225..+225]} after U{Vermeille
1106 <https://Search.ProQuest.com/docview/639493848>} (2004), p 95.
1108 @see: Property C{lamVermeille}.
1109 '''
1110 return Lon(Vermeille=degrees(self.lamVermeille))
1112 def _T_x_M(self, T):
1113 '''(INTERNAL) Update M{self.M = T.multiply(self.M)}.
1114 '''
1115 t = list(self)
1116 M = self._IndexM
1117 t[M] = T.multiply(t[M])
1118 return self.classof(*t)
1120 @Property_RO
1121 def phi(self):
1122 '''Get the latitude in C{radians} (C{float}).
1123 '''
1124 return self.philam.phi
1126 @Property_RO
1127 def philam(self):
1128 '''Get the lat-, longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).
1129 '''
1130 return PhiLam2Tuple(radians(self.lat), radians(self.lon), name=self.name)
1132 @Property_RO
1133 def philamheight(self):
1134 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
1135 '''
1136 return self.philam.to3Tuple(self.height)
1138 @Property_RO
1139 def philamheightdatum(self):
1140 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLam4Tuple}C{(phi, lam, height, datum)}).
1141 '''
1142 return self.philamheight.to4Tuple(self.datum)
1144 @Property_RO
1145 def philamVermeille(self):
1146 '''Get the latitude and I{Vermeille} longitude in C{radians [-PI*3/2..+PI*3/2]} (L{PhiLam2Tuple}C{(phi, lam)}).
1148 @see: Property C{lamVermeille}.
1149 '''
1150 return PhiLam2Tuple(radians(self.lat), self.lamVermeille, name=self.name)
1152 def toCartesian(self, Cartesian=None, **Cartesian_kwds):
1153 '''Return the geocentric C{(x, y, z)} coordinates as an ellipsoidal or spherical
1154 C{Cartesian}.
1156 @kwarg Cartesian: Optional class to return C{(x, y, z)} (L{ellipsoidalKarney.Cartesian},
1157 L{ellipsoidalNvector.Cartesian}, L{ellipsoidalVincenty.Cartesian},
1158 L{sphericalNvector.Cartesian} or L{sphericalTrigonometry.Cartesian})
1159 or C{None}.
1160 @kwarg Cartesian_kwds: Optional, additional B{C{Cartesian}} keyword arguments, ignored
1161 if C{B{Cartesian} is None}.
1163 @return: A C{B{Cartesian}(x, y, z, **B{Cartesian_kwds})} instance or
1164 a L{Vector4Tuple}C{(x, y, z, h)} if C{B{Cartesian} is None}.
1166 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}}.
1167 '''
1168 if Cartesian in (None, Vector4Tuple):
1169 r = self.xyzh
1170 elif Cartesian is Vector3Tuple:
1171 r = self.xyz
1172 else:
1173 _xsubclassof(self._CartesianBase, Cartesian=Cartesian)
1174 r = Cartesian(self, **_xkwds(Cartesian_kwds, name=self.name))
1175 return r
1177 def toDatum(self, datum2):
1178 '''Convert this C{Ecef9Tuple} to an other datum.
1180 @arg datum2: Datum to convert I{to} (L{Datum}).
1182 @return: The converted 9-Tuple (C{Ecef9Tuple}).
1184 @raise TypeError: The B{C{datum2}} is not a L{Datum}.
1185 '''
1186 if self.datum in (None, datum2): # PYCHOK _Names_
1187 r = self.copy()
1188 else:
1189 c = self._CartesianBase(self, datum=self.datum, name=self.name) # PYCHOK _Names_
1190 # c.toLatLon converts datum, x, y, z, lat, lon, etc.
1191 # and returns another Ecef9Tuple iff LatLon is None
1192 r = c.toLatLon(datum=datum2, LatLon=None)
1193 return r
1195 def toLatLon(self, LatLon=None, **LatLon_kwds):
1196 '''Return the geodetic C{(lat, lon, height[, datum])} coordinates.
1198 @kwarg LatLon: Optional class to return C{(lat, lon, height[, datum])}
1199 or C{None}.
1200 @kwarg LatLon_kwds: Optional B{C{height}}, B{C{datum}} and other
1201 B{C{LatLon}} keyword arguments.
1203 @return: An instance of C{B{LatLon}(lat, lon, **B{LatLon_kwds})}
1204 or if B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon,
1205 height)} respectively L{LatLon4Tuple}C{(lat, lon, height,
1206 datum)} depending on whether C{datum} is un-/specified.
1208 @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}}.
1209 '''
1210 kwds = _xkwds(LatLon_kwds, height=self.height, datum=self.datum, name=self.name) # PYCHOK Ecef9Tuple
1211 d = kwds[_datum_]
1212 if LatLon is None:
1213 r = LatLon3Tuple(self.lat, self.lon, kwds[_height_], name=kwds[_name_]) # PYCHOK Ecef9Tuple
1214 if d:
1215 r = r.to4Tuple(d) # checks type(d)
1216 else:
1217 if d is None: # remove None datum
1218 _ = kwds.pop[_datum_]
1219 r = LatLon(self.lat, self.lon, **kwds) # PYCHOK Ecef9Tuple
1220 _xdatum(getattr(r, _datum_, self.datum), self.datum) # PYCHOK Ecef9Tuple
1221 return r
1223 def toLocal(self, ltp, Xyz=None, **Xyz_kwds):
1224 '''Convert this geocentric to I{local} C{x}, C{y} and C{z}.
1226 @kwarg ltp: The I{local tangent plane} (LTP) to use (L{Ltp}).
1227 @kwarg Xyz: Optional class to return C{x}, C{y} and C{z}
1228 (L{XyzLocal}, L{Enu}, L{Ned}) or C{None}.
1229 @kwarg Xyz_kwds: Optional, additional B{C{Xyz}} keyword
1230 arguments, ignored if C{B{Xyz} is None}.
1232 @return: An B{C{Xyz}} instance or if C{B{Xyz} is None},
1233 a L{Local9Tuple}C{(x, y, z, lat, lon, height,
1234 ltp, ecef, M)} with C{M=None}, always.
1236 @raise TypeError: Invalid B{C{ltp}}.
1237 '''
1238 return _MODS.ltp._xLtp(ltp)._ecef2local(self, Xyz, Xyz_kwds)
1240 def toVector(self, Vector=None, **Vector_kwds):
1241 '''Return the geocentric C{(x, y, z)} coordinates as vector.
1243 @kwarg Vector: Optional vector class to return C{(x, y, z)} or
1244 C{None}.
1245 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword
1246 arguments, ignored if C{B{Vector} is None}.
1248 @return: A C{Vector}C{(x, y, z, **Vector_kwds)} instance or a
1249 L{Vector3Tuple}C{(x, y, z)} if B{C{Vector}} is C{None}.
1251 @see: Propertes C{xyz} and C{xyzh}
1252 '''
1253 return self.xyz if Vector is None else self._xnamed(
1254 Vector(self.x, self.y, self.z, **Vector_kwds)) # PYCHOK Ecef9Tuple
1256 @Property_RO
1257 def xyz(self):
1258 '''Get the geocentric C{(x, y, z)} coordinates (L{Vector3Tuple}C{(x, y, z)}).
1259 '''
1260 return Vector3Tuple(self.x, self.y, self.z, name=self.name)
1262 @Property_RO
1263 def xyzh(self):
1264 '''Get the geocentric C{(x, y, z)} coordinates and C{height} (L{Vector4Tuple}C{(x, y, z, h)})
1265 '''
1266 return self.xyz.to4Tuple(self.height)
1269def _4Ecef(this, Ecef): # in .datums.Datum.ecef, .ellipsoids.Ellipsoid.ecef
1270 '''Return an ECEF converter for C{this} L{Datum} or L{Ellipsoid}.
1271 '''
1272 if Ecef is None:
1273 Ecef = EcefKarney
1274 else:
1275 _xinstanceof(*_Ecefs, Ecef=Ecef)
1276 return Ecef(this, name=this.name)
1279def _xEcef(Ecef): # PYCHOK .latlonBase.py
1280 '''(INTERNAL) Validate B{C{Ecef}} I{class}.
1281 '''
1282 if issubclassof(Ecef, _EcefBase):
1283 return Ecef
1284 raise _TypesError(_Ecef_, Ecef, *_Ecefs)
1287_Ecefs = (EcefKarney, EcefSudano, EcefVeness, EcefYou,
1288 EcefFarrell21, EcefFarrell22)
1290__all__ += _ALL_DOCS(_EcefBase)
1292# **) MIT License
1293#
1294# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1295#
1296# Permission is hereby granted, free of charge, to any person obtaining a
1297# copy of this software and associated documentation files (the "Software"),
1298# to deal in the Software without restriction, including without limitation
1299# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1300# and/or sell copies of the Software, and to permit persons to whom the
1301# Software is furnished to do so, subject to the following conditions:
1302#
1303# The above copyright notice and this permission notice shall be included
1304# in all copies or substantial portions of the Software.
1305#
1306# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1307# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1308# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1309# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1310# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1311# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1312# OTHER DEALINGS IN THE SOFTWARE.