Coverage for pygeodesy/ecef.py: 95%

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, _isNAN, _isNAN0, 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.05.19'

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]

101

102

103class EcefError(_ValueError):

104 '''An ECEF or C{Ecef*} related issue.

105 '''

106 pass

107

108

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

116

117 def __init__(self, a_ellipsoid=_EWGS84, f=None, lon00=INT0, **name):

118 '''New C{Ecef*} converter.

119

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}).

129

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_

143

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

151

152 except (TypeError, ValueError) as x:

153 t = unstr(self.classname, a=a_ellipsoid, f=f)

154 raise EcefError(_SPACE_(t, _ellipsoid_), cause=x)

155

156 if name:

157 self.name = name

158 if lon00 is not INT0:

159 self.lon00 = lon00

160

161 def __eq__(self, other):

162 '''Compare this and an other Ecef.

163

164 @arg other: The other ecef (C{Ecef*}).

165

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)

170

171 @Property_RO

172 def datum(self):

173 '''Get the datum (L{Datum}).

174 '''

175 return self._datum

176

177 @Property_RO

178 def ellipsoid(self):

179 '''Get the ellipsoid (L{Ellipsoid} or L{Ellipsoid2}).

180 '''

181 return self._E

182

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

188

189 a = equatorialRadius = equatoradius # Karney property

190

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

197

198 f = flattening

199

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)

209

210 E = self.ellipsoid

211 n = E.roc1_(sa, ca) if self._isYou else E.roc1_(sa)

212 H = _isNAN0(h)

213 c = (H + n) * ca

214 x = cb * c

215 y = sb * c

216 z = (H + n * E.e21) * sa

217

218 m = self._Matrix(sa, ca, sb, cb) if M else None

219 return Ecef9Tuple(x, y, z, lat, lon, h,

220 0, m, self.datum, # C=0, always

221 name=self._name__(name))

222

223 def forward(self, latlonh, lon=None, height=0, M=False, **name):

224 '''Convert from geodetic C{(lat, lon, height)} to geocentric C{(x, y, z)}.

225

226 @arg latlonh: Either a C{LatLon}, an L{Ecef9Tuple} or C{scalar}

227 latitude (C{degrees}).

228 @kwarg lon: Optional C{scalar} longitude for C{scalar} B{C{latlonh}}

229 (C{degrees}).

230 @kwarg height: Optional height (C{meter}), vertically above (or below)

231 the surface of the ellipsoid.

232 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).

233 @kwarg name: Optional C{B{name}=NN} (C{str}).

234

235 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with

236 geocentric C{(x, y, z)} coordinates for the given geodetic ones

237 C{(lat, lon, height)}, case C{C} 0, optional C{M} (L{EcefMatrix})

238 and C{datum} if available.

239

240 @raise EcefError: If B{C{latlonh}} not C{LatLon}, L{Ecef9Tuple} or

241 C{scalar} or B{C{lon}} not C{scalar} for C{scalar}

242 B{C{latlonh}} or C{abs(lat)} exceeds 90°.

243

244 @note: Use method C{.forward_} to specify C{lat} and C{lon} in C{radians}

245 and avoid double angle conversions.

246 '''

247 llhn = _llhn4(latlonh, lon, height, **name)

248 return self._forward(*llhn, M=M)

249

250 def forward_(self, phi, lam, height=0, M=False, **name):

251 '''Like method C{.forward} except with geodetic lat- and longitude given

252 in I{radians}.

253

254 @arg phi: Latitude in I{radians} (C{scalar}).

255 @arg lam: Longitude in I{radians} (C{scalar}).

256 @kwarg height: Optional height (C{meter}), vertically above (or below)

257 the surface of the ellipsoid.

258 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).

259 @kwarg name: Optional C{B{name}=NN} (C{str}).

260

261 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}

262 with C{lat} set to C{degrees90(B{phi})} and C{lon} to

263 C{degrees180(B{lam})}.

264

265 @raise EcefError: If B{C{phi}} or B{C{lam}} invalid or not C{scalar}.

266 '''

267 try: # like function C{_llhn4} below

268 plhn = Phi(phi), Lam(lam), Height(height), _name__(name)

269 except (TypeError, ValueError) as x:

270 raise EcefError(phi=phi, lam=lam, height=height, cause=x)

271 return self._forward(*plhn, M=M, _philam=True)

272

273 @property_ROver

274 def _Geocentrics(self):

275 '''(INTERNAL) Get the valid geocentric classes. I{once}.

276 '''

277 return (Ecef9Tuple, # overwrite property_ROver

278 _MODS.vector3d.Vector3d) # _MODS.cartesianBase.CartesianBase

279

280 @property

281 def lon00(self):

282 '''Get the I{"polar"} longitude (C{degrees}), see method C{reverse}.

283 '''

284 return self._lon00

285

286 @lon00.setter # PYCHOK setter!

287 def lon00(self, lon00):

288 '''Set the I{"polar"} longitude (C{degrees}), see method C{reverse}.

289 '''

290 self._lon00 = Degrees(lon00=lon00)

291

292 def _Matrix(self, *sa_ca_sb_cb):

293 '''(INTERNAL) Create a rotation L{EcefMatrix}.

294 '''

295 return self._xnamed(EcefMatrix(*sa_ca_sb_cb))

296

297 def _polon(self, y, x, p, **lon00_name):

298 '''(INTERNAL) Handle I{"polar"} longitude.

299 '''

300 return atan2d(y, x) if p else _xkwds_get(lon00_name, lon00=self.lon00)

301

302 def reverse(self, xyz, y=None, z=None, M=False, **tol_lon00_name): # PYCHOK no cover

303 '''I{Must be overloaded}.'''

304 self._notOverloaded(xyz, y=y, z=z, M=M, **tol_lon00_name)

305

306 def _reversError(self, x, y, z, r, tol):

307 '''(INTERNAL) Convergence error.

308 '''

309 t = unstr(self.reverse, x=x, y=y, z=z)

310 return EcefError(t, txt=Fmt.no_convergence(r, tol))

311

312 def toStr(self, prec=9, **unused): # PYCHOK signature

313 '''Return this C{Ecef*} as a string.

314

315 @kwarg prec: Precision, number of decimal digits (0..9).

316

317 @return: This C{Ecef*} (C{str}).

318 '''

319 return self.attrs(_a_, _f_, _datum_, _name_, prec=prec) # _ellipsoid_

320

321 def _xyzllhCpn9(self, xyz, y, z, **lon00_name):

322 '''(INTERNAL) Get C{x, y, z} and determine case C{C}, C{lat}, C{lon}, etc.

323 '''

324 x, y, z, n = _xyzn4(xyz, y, z, self._Geocentrics, **lon00_name)

325

326 s, c, p = _norm3(y, x) # distance to polar axis

327 if p < EPS0: # polar

328 lat = copysign0(_90_0, z)

329 h = fabs(z) - self.ellipsoid.b

330 C = 2

331 p = 0 # force lon00

332 elif fabs(z) < EPS0: # equatorial

333 lat = _0_0

334 h = p - self.ellipsoid.a

335 C = 3

336 else: # see _norm7 below, EcefKarney

337 h = hypot(z, p) # distance to earth center

338 if h > self.ellipsoid._heightMax:

339 lat = atan1d(z / h, p / h)

340 C = 4 # too high

341# p *= _0_5

342 else: # pass h for EcefVeness

343 lat = None # -90..90

344 C = 1 # normal

345 lon = self._polon(s, c, p, **lon00_name)

346 return x, y, z, lat, lon, h, C, p, self._name__(n)

347

348

349class EcefFarrell21(_EcefBase):

350 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)

351 coordinates based on I{Jay A. Farrell}'s U{Table 2.1<https://Books.Google.com/

352 books?id=fW4foWASY6wC>}, page 29.

353 '''

354

355 def reverse(self, xyz, y=None, z=None, M=None, **lon00_name): # PYCHOK unused M

356 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using

357 I{Farrell}'s U{Table 2.1<https://Books.Google.com/books?id=fW4foWASY6wC>},

358 page 29, aka the I{Heikkinen application} of U{Ferrari's solution

359 <https://WikiPedia.org/wiki/Geographic_coordinate_conversion>}.

360

361 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}

362 coordinate (C{meter}).

363 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).

364 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).

365 @kwarg M: I{Ignored}, rotation matrix C{M} not available.

366 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument

367 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude

368 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.

369

370 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with

371 geodetic coordinates C{(lat, lon, height)} for the given geocentric

372 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum}.

373

374 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}

375 not C{scalar} for C{scalar} B{C{xyz}} or C{sqrt} domain or

376 zero division error.

377

378 @see: L{EcefFarrell22} and L{EcefVeness}.

379 '''

380 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)

381 if lat is None:

382 E = self.ellipsoid

383 a = E.a

384 a2 = E.a2

385 b2 = E.b2

386 e2 = E.e2

387 e2_ = E.e2abs * E.a2_b2 # (E.e * E.a_b)**2 = 0.0820944... WGS84

388 e4 = E.e4

389

390 try: # names as page 29

391 z2 = z**2

392 ez = z2 * (_1_0 - e2) # E.e2s2(z)

393

394 p2 = p**2

395 G = p2 + ez - e2 * (a2 - b2) # p2 + ez - e4 * a2

396 F = b2 * z2 * 54

397 c = e4 * p2 * F / G**3

398 s = sqrt(c * (c + _2_0))

399 c = cbrt(s + c + _1_0)

400 G *= fsumf_(c, _1_0, _1_0 / c) # k

401 P = F / (G**2 * _3_0)

402 Q = sqrt(_2_0 * e4 * P + _1_0)

403 Q1 = Q + _1_0

404 s = fsumf_(a2 * (Q1 / Q) * _0_5,

405 -P * ez / (Q * Q1),

406 -P * p2 * _0_5)

407 r = p * P * e2 / Q1 - sqrt(s)

408 r = p + r * e2

409 v = b2 / (sqrt(r**2 + ez) * a) # z0 / z

410

411 h = hypot(r, z) * (_1_0 - v)

412 lat = atan1d((e2_ * v + _1_0) * z, p)

413 # note, phi and lam are swapped on page 29

414

415 except (ValueError, ZeroDivisionError) as X:

416 raise EcefError(x=x, y=y, z=z, cause=X)

417

418 return Ecef9Tuple(x, y, z, lat, lon, h,

419 C, None, self.datum, name=name)

420

421

422class EcefFarrell22(_EcefBase):

423 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)

424 coordinates based on I{Jay A. Farrell}'s U{Table 2.2<https://Books.Google.com/

425 books?id=fW4foWASY6wC>}, page 30.

426 '''

427

428 def reverse(self, xyz, y=None, z=None, M=None, **lon00_name): # PYCHOK unused M

429 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using

430 I{Farrell}'s U{Table 2.2<https://Books.Google.com/books?id=fW4foWASY6wC>},

431 page 30.

432

433 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}

434 coordinate (C{meter}).

435 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).

436 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).

437 @kwarg M: I{Ignored}, rotation matrix C{M} not available.

438 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument

439 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude

440 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.

441

442 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with

443 geodetic coordinates C{(lat, lon, height)} for the given geocentric

444 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum}.

445

446 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}

447 not C{scalar} for C{scalar} B{C{xyz}} or C{sqrt} domain or

448 zero division error.

449

450 @see: L{EcefFarrell21} and L{EcefVeness}.

451 '''

452 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)

453 if lat is None:

454 E = self.ellipsoid

455 a, b = E.a, E.b

456 s, c, _ = _norm3(z * a, p * b) # Bowring

457 s, c, _ = _norm3(z + s**3 * b * E.e22,

458 p - c**3 * a * E.e2)

459 lat = atan1d(s, c)

460 h = (p / fabs(c) - (E.roc1_(s) if s else a)) if c else (fabs(z) - b)

461 # note, phi and lam are swapped on page 30

462

463 return Ecef9Tuple(x, y, z, lat, lon, h,

464 C, None, self.datum, name=name)

465

466

467class EcefKarney(_EcefBase):

468 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF)

469 coordinates transcoded from I{Karney}'s C++ U{Geocentric<https://GeographicLib.SourceForge.io/

470 C++/doc/classGeographicLib_1_1Geocentric.html>} methods.

471

472 @note: On methods C{.forward} and C{.forwar_}, let C{v} be a unit vector located

473 at C{(lat, lon, h)}. We can express C{v} as column vectors in one of two

474 ways, C{v1} in East, North, Up (ENU) coordinates (where the components are

475 relative to a local coordinate system at C{C(lat0, lon0, h0)}) or as C{v0}

476 in geocentric C{x, y, z} coordinates. Then, M{v0 = M ⋅ v1} where C{M} is

477 the rotation matrix.

478 '''

479

480 def reverse(self, xyz, y=None, z=None, M=False, **lon00_name):

481 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}.

482

483 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}

484 coordinate (C{meter}).

485 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).

486 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).

487 @kwarg M: Optionally, return the rotation L{EcefMatrix} (C{bool}).

488 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument

489 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude

490 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.

491

492 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with

493 geodetic coordinates C{(lat, lon, height)} for the given geocentric

494 ones C{(x, y, z)}, case C{C}, optional C{M} (L{EcefMatrix}) and

495 C{datum} if available.

496

497 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}

498 not C{scalar} for C{scalar} B{C{xyz}}.

499

500 @note: In general, there are multiple solutions and the result which minimizes

501 C{height} is returned, i.e., the C{(lat, lon)} corresponding to the

502 closest point on the ellipsoid. If there are still multiple solutions

503 with different latitudes (applies only if C{z} = 0), then the solution

504 with C{lat} > 0 is returned. If there are still multiple solutions with

505 different longitudes (applies only if C{x} = C{y} = 0), then C{lon00} is

506 returned. The returned C{lon} is in the range [−180°, 180°] and C{height}

507 is not below M{−E.a * (1 − E.e2) / sqrt(1 − E.e2 * sin(lat)**2)}. Like

508 C{forward} above, M{v1 = Transpose(M) ⋅ v0}.

509 '''

510 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, **lon00_name)

511

512 E = self.ellipsoid

513 f = E.f

514

515 sa, ca, sb, cb, R, h, C = _norm7(y, x, z, E)

516 if C: # PYCHOK no cover

517 pass # too high, far

518 elif E.e4: # E.isEllipsoidal

519 # Treat prolate spheroids by swapping R and Z here and by

520 # switching the arguments to phi = atan2(...) at the end.

521 p = (R / E.a)**2

522 q = (z / E.a)**2 * E.e21

523 if f < 0:

524 p, q = q, p

525 r = fsumf_(p, q, -E.e4)

526 e = E.e4 * q

527 if e or r > 0:

528 # Avoid possible division by zero when r = 0 by multiplying

529 # equations for s and t by r^3 and r, respectively.

530 s = d = e * p / _4_0 # s = r^3 * s

531 u = r = r / _6_0

532 r2 = r**2

533 r3 = r2 * r

534 t3 = r3 + s

535 d *= t3 + r3

536 if d < 0:

537 # t is complex, but the way u is defined, the result is real.

538 # There are three possible cube roots. We choose the root

539 # which avoids cancellation. Note, d < 0 implies r < 0.

540 u += cos(atan2(sqrt(-d), -t3) / _3_0) * r * _2_0

541 else:

542 # Pick the sign on the sqrt to maximize abs(t3). This

543 # minimizes loss of precision due to cancellation. The

544 # result is unchanged because of the way the t is used

545 # in definition of u.

546 if d > 0:

547 t3 += copysign0(sqrt(d), t3) # t3 = (r * t)^3

548 # N.B. cbrt always returns the real root, cbrt(-8) = -2.

549 t = cbrt(t3) # t = r * t

550 if t: # t can be zero; but then r2 / t -> 0.

551 u = fsumf_(u, t, r2 / t)

552 v = sqrt(e + u**2) # guaranteed positive

553 # Avoid loss of accuracy when u < 0. Underflow doesn't occur in

554 # E.e4 * q / (v - u) because u ~ e^4 when q is small and u < 0.

555 u = (e / (v - u)) if u < 0 else (u + v) # u+v, guaranteed positive

556 # Need to guard against w going negative due to roundoff in u - q.

557 w = E.e2abs * (u - q) / (_2_0 * v)

558 # Rearrange expression for k to avoid loss of accuracy due to

559 # subtraction. Division by 0 not possible because u > 0, w >= 0.

560 k1 = k2 = (u / (sqrt(u + w**2) + w)) if w > 0 else sqrt(u)

561 if f < 0:

562 k1 -= E.e2

563 else:

564 k2 += E.e2

565 sa, ca, h = _norm3(z / k1, R / k2)

566 h *= k1 - E.e21

567 C = 1

568

569 else: # e = E.e4 * q == 0 and r <= 0

570 # This leads to k = 0 (oblate, equatorial plane) and k + E.e^2 = 0

571 # (prolate, rotation axis) and the generation of 0/0 in the general

572 # formulas for phi and h, using the general formula and division

573 # by 0 in formula for h. Handle this case by taking the limits:

574 # f > 0: z -> 0, k -> E.e2 * sqrt(q) / sqrt(E.e4 - p)

575 # f < 0: r -> 0, k + E.e2 -> -E.e2 * sqrt(q) / sqrt(E.e4 - p)

576 q = E.e4 - p

577 if f < 0:

578 p, q = q, p

579 e = E.a

580 else:

581 e = E.b2_a

582 sa, ca, h = _norm3(sqrt(q * E._1_e21), sqrt(p))

583 if z < 0: # for tiny negative z, not for prolate

584 sa = neg(sa)

585 h *= neg(e / E.e2abs)

586 C = 3

587

588 else: # E.e4 == 0, spherical case

589 # Dealing with underflow in the general case with E.e2 = 0 is

590 # difficult. Origin maps to North pole, same as with ellipsoid.

591 sa, ca, _ = _norm3((z if h else _1_0), R)

592 h -= E.a

593 C = 2

594

595 # lon00 <https://GitHub.com/mrJean1/PyGeodesy/issues/77>

596 lon = self._polon(sb, cb, R, **lon00_name)

597 m = self._Matrix(sa, ca, sb, cb) if M else None

598 return Ecef9Tuple(x, y, z, atan1d(sa, ca), lon, h,

599 C, m, self.datum, name=self._name__(name))

600

601

602class EcefSudano(_EcefBase):

603 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates

604 based on I{John J. Sudano}'s U{paper<https://www.ResearchGate.net/publication/3709199>}.

605 '''

606 _tol = EPS # DEPRECATED

607

608 def reverse(self, xyz, y=None, z=None, M=None, tol=EPS, **lon00_name): # PYCHOK unused M

609 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using

610 I{Sudano}'s U{iterative method<https://www.ResearchGate.net/publication/3709199>}.

611

612 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}

613 coordinate (C{meter}).

614 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).

615 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).

616 @kwarg M: I{Ignored}, rotation matrix C{M} not available.

617 @kwarg tol: Convergence tolerance for C{sin}(latitude) (C{scalar}).

618 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument

619 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude

620 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.

621

622 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with geodetic

623 coordinates C{(lat, lon, height)} for the given geocentric ones C{(x, y, z)},

624 case C{C}, C{M=None} always and C{datum} if available.

625

626 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}

627 not C{scalar} for C{scalar} B{C{xyz}} or no convergence.

628

629 @see: Class L{EcefUPC}.

630 '''

631 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)

632 if lat is None:

633 E = self.ellipsoid

634 e = E.e2 * E.a

635 d = e - p

636

637 sa, ca, _ = _norm3(fabs(z), p * E.e21)

638 _S2 = Fsum(sa).fsum2f_

639 tol = Scalar_(tol=tol, low=self.tolerance, Error=EcefError)

640 # Sudano's Eq (A-6) and (A-7) refactored/reduced,

641 # replacing Rn from Eq (A-4) with n = E.a / ca:

642 # N = ca**2 * ((z + E.e2 * n * sa) * ca - p * sa)

643 # = ca**2 * (z * ca + E.e2 * E.a * sa - p * sa)

644 # = ca**2 * (z * ca + (E.e2 * E.a - p) * sa)

645 # D = ca**3 * (E.e2 * n / E.e2s2(sa)) - p

646 # = ca**2 * (E.e2 * E.a / E.e2s2(sa) - p / ca**2)

647 # N / D = (z * ca + (E.e2 * E.a - p) * sa) /

648 # (E.e2 * E.a / E.e2s2(sa) - p / ca**2)

649 for i in range(1, _TRIPS): # 6+ max

650 ca2 = _1_0 - sa**2

651 if ca2 < EPS02:

652 break

653 D = p / ca2 - e / E.e2s2(sa)

654 if fabs(D) < EPS0:

655 break

656 ca = sqrt(ca2)

657 sa, D = _S2(z * ca / D, d * sa / D)

658 if fabs(D) < tol:

659 break

660 else: # PYCHOK no cover

661 raise self._reversError(x, y, z, fabs(D), tol)

662

663 sa = copysign0(sa, z)

664 lat = atan1d(sa, ca)

665 # h = (fabs(z) + p - E.a * cos(a + E.e21) * sa / ca) / (ca + sa)

666 # Sudano's Eq (7) doesn't produce the correct height, ...

667 h = E._heightB(sa, ca, z, p) # ... use Veness' (Bowring eqn 7)

668 else:

669 i = 0

670 return Ecef9Tuple(x, y, z, lat, lon, h,

671 C, None, self.datum, # M=None

672 iteration=i, name=name)

673

674 @deprecated_property

675 def tolerance(self):

676 '''DEPRECATED on 2025.08.22, use keyword argument C{tol}.'''

677 return self._tol

678

679 @tolerance.setter # PYCHOK setter!

680 def tolerance(self, tol):

681 self._tol = Scalar_(tolerance=tol, low=EPS, Error=EcefError)

682

683

684class EcefUPC(_EcefBase):

685 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates based on

686 I{UPC}'s U{method<https://GSSC.ESA.int/navipedia/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion>}.

687 '''

688

689 def reverse(self, xyz, y=None, z=None, M=None, tol=_TOL, **lon00_name): # PYCHOK unused M

690 '''Convert from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)} using I{UPC}'s

691 U{iterative method<https://GSSC.ESA.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf>}, page 186.

692

693 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}

694 coordinate (C{meter}).

695 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).

696 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).

697 @kwarg M: I{Ignored}, rotation matrix C{M} not available.

698 @kwarg tol: Convergence tolerance for the latitude (C{degrees}).

699 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument

700 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude

701 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.

702

703 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with geodetic

704 coordinates C{(lat, lon, height)} for the given geocentric ones C{(x, y, z)},

705 case C{C}, C{M=None} always and C{datum} if available.

706

707 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}

708 not C{scalar} for C{scalar} B{C{xyz}} or no convergence.

709

710 @see: Class L{EcefSudano}.

711 '''

712 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)

713 if lat is None:

714 t = Degrees_(tol=tol, low=EPS, Error=EcefError).toRadians()

715 E = self.ellipsoid

716 a = E.a

717 e2 = E.e2 # signed

718

719 z_ = fabs(z)

720 ph_ = atan1(z_, E.e21 * p)

721 for i in range(1, _TRIPS): # 5..6 max

722 s, c = sincos2(ph_)

723 N = a / sqrt(_1_0 - s**2 * e2)

724 ph = atan1(z_, p - N * c * e2)

725 r = fabs(ph - ph_)

726 if r < t:

727 lat = copysign0(degrees(ph), z)

728 h = p / c - N

729 break

730 ph_ = ph

731 else: # PYCHOK no cover

732 raise self._reversError(x, y, z, degrees(r), tol)

733 else:

734 i = 0

735 return Ecef9Tuple(x, y, z, lat, lon, h,

736 C, None, self.datum, # M=None

737 iteration=i, name=name)

738

739

740class EcefVeness(_EcefBase):

741 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates

742 transcoded from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal, Cartesian<https://

743 www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}.

744

745 @see: U{I{A Guide to Coordinate Systems in Great Britain}<https://www.OrdnanceSurvey.co.UK/

746 documents/resources/guide-coordinate-systems-great-britain.pdf>}, section I{B) Converting

747 between 3D Cartesian and ellipsoidal latitude, longitude and height coordinates}.

748 '''

749

750 def reverse(self, xyz, y=None, z=None, M=None, **lon00_name): # PYCHOK unused M

751 '''Conversion from geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}

752 transcoded from I{Chris Veness}' U{JavaScript<https://www.Movable-Type.co.UK/

753 scripts/geodesy/docs/latlon-ellipsoidal.js.html>}.

754

755 Uses B. R. Bowring’s formulation for μm precision in concise form U{I{The accuracy

756 of geodetic latitude and height equations}<https://www.ResearchGate.net/publication/

757 233668213>}, Survey Review, Vol 28, 218, Oct 1985.

758

759 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}

760 coordinate (C{meter}).

761 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).

762 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).

763 @kwarg M: I{Ignored}, rotation matrix C{M} not available.

764 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument

765 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude

766 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.

767

768 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with

769 geodetic coordinates C{(lat, lon, height)} for the given geocentric

770 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum} if available.

771

772 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or B{C{z}}

773 not C{scalar} for C{scalar} B{C{xyz}}.

774

775 @see: Toms, Ralph M. U{I{An Efficient Algorithm for Geocentric to Geodetic

776 Coordinate Conversion}<https://www.OSTI.gov/scitech/biblio/110235>},

777 Sept 1995 and U{I{An Improved Algorithm for Geocentric to Geodetic

778 Coordinate Conversion}<https://www.OSTI.gov/scitech/servlets/purl/231228>},

779 Apr 1996, both from Lawrence Livermore National Laboratory (LLNL).

780 '''

781 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)

782 if lat is None:

783 E = self.ellipsoid

784 a = E.a

785 B = E.b * E.e22

786 # parametric latitude (Bowring eqn 17, replaced)

787 t = (E.b * z) / (a * p) * (_1_0 + B / h) # h = hypot(z, p)

788 c = _1_0 / hypot1(t)

789 s = c * t

790 # geodetic latitude (Bowring eqn 18)

791 s, c, _ = _norm3(z + s**3 * B,

792 p - c**3 * a * E.e2)

793 lat = atan1d(s, c)

794 h = E._heightB(s, c, z, p) # height (Bowring eqn 7)

795

796 return Ecef9Tuple(x, y, z, lat, lon, h,

797 C, None, self.datum, name=name)

798

799

800class EcefYou(_EcefBase):

801 '''Conversion between geodetic and geocentric, I{Earth-Centered, Earth-Fixed} (ECEF) coordinates

802 using I{Rey-Jer You}'s U{transformation<https://www.ResearchGate.net/publication/240359424>}

803 for I{non-prolate} ellipsoids.

804

805 @see: Featherstone, W.E., Claessens, S.J. U{I{Closed-form transformation between geodetic and

806 ellipsoidal coordinates}<https://Espace.Curtin.edu.AU/bitstream/handle/20.500.11937/

807 11589/115114_9021_geod2ellip_final.pdf>} Studia Geophysica et Geodaetica, 2008, 52,

808 pages 1-18 and U{PyMap3D <https://PyPI.org/project/pymap3d>}.

809 '''

810 _isYou = True

811

812 def __init__(self, a_ellipsoid=_EWGS84, f=None, **lon00_name): # PYCHOK signature

813 _EcefBase.__init__(self, a_ellipsoid, f=f, **lon00_name) # inherited documentation

814

815 E = self.ellipsoid

816 e2 = E.a2 - E.b2

817 if e2 < 0 or E.f < 0:

818 raise EcefError(ellipsoid=E, txt=_prolate_)

819 self._ee2 = sqrt0(e2), e2

820

821 def reverse(self, xyz, y=None, z=None, M=None, **lon00_name): # PYCHOK unused M

822 '''Convert geocentric C{(x, y, z)} to geodetic C{(lat, lon, height)}

823 using I{Rey-Jer You}'s transformation.

824

825 @arg xyz: A geocentric (C{Cartesian}, L{Ecef9Tuple}) or C{scalar} ECEF C{x}

826 coordinate (C{meter}).

827 @kwarg y: ECEF C{y} coordinate for C{scalar} B{C{xyz}} and B{C{z}} (C{meter}).

828 @kwarg z: ECEF C{z} coordinate for C{scalar} B{C{xyz}} and B{C{y}} (C{meter}).

829 @kwarg M: I{Ignored}, rotation matrix C{M} not available.

830 @kwarg lon00_name: Optional C{B{name}=NN} (C{str}) and optional keyword argument

831 C{B{lon00}=INT0} (C{degrees}), an arbitrary I{"polar"} longitude

832 returned if C{B{x}=0} and C{B{y}=0}, see property C{lon00}.

833

834 @return: An L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with

835 geodetic coordinates C{(lat, lon, height)} for the given geocentric

836 ones C{(x, y, z)}, case C{C}, C{M=None} always and C{datum} if available.

837

838 @raise EcefError: Invalid B{C{xyz}} or C{scalar} C{x} or B{C{y}} and/or

839 B{C{z}} not C{scalar} for C{scalar} B{C{xyz}} or the

840 ellipsoid is I{prolate}.

841 '''

842 x, y, z, lat, lon, h, C, p, name = self._xyzllhCpn9(xyz, y, z, **lon00_name)

843 if lat is None:

844 E = self.ellipsoid

845 a, b = E.a, E.b

846 e, e2 = self._ee2

847

848 u = hypot2_(x, y, z) - e2

849 u += hypot(u, e * z * _2_0)

850 u *= _0_5

851 if u > EPS02:

852 u = sqrt(u)

853 q = hypot(u, e)

854 B = atan1(q * z, u * p) # beta0 = atan(q / u * z / p)

855 sB, cB = sincos2(B)

856 if cB and sB:

857 q *= a

858 d = (q / cB - e2 * cB) / sB

859 if isnon0(d):

860 B += fsumf_(u * b, -q, e2) / d

861 sB, cB = sincos2(B)

862 elif u < (-EPS02):

863 raise EcefError(x=x, y=y, z=z, u=u, txt=_singular_)

864 else: # near polar # PYCHOK no cover

865 sB, cB, C = _copysign_1_0(z), _0_0, 2

866

867 lat = atan1d( a * sB, b * cB) # atan(E.a_b * tan(B))

868 h = hypot(p - a * cB, z - b * sB)

869 if hypot2_(x, y, z * E.a_b) < E.a2: # or lat < 0 or z < 0

870 h = neg(h) # inside ellipsoid

871

872 return Ecef9Tuple(x, y, z, lat, lon, h,

873 C, None, self.datum, name=name)

874

875

876class EcefMatrix(_NamedTuple):

877 '''A rotation matrix known as I{East-North-Up (ENU) to ECEF}.

878

879 @see: U{From ENU to ECEF<https://WikiPedia.org/wiki/

880 Geographic_coordinate_conversion#From_ECEF_to_ENU>} and

881 U{Issue #74<https://Github.com/mrJean1/PyGeodesy/issues/74>}.

882 '''

883 _Names_ = ('_0_0_', '_0_1_', '_0_2_', # row-order

884 '_1_0_', '_1_1_', '_1_2_',

885 '_2_0_', '_2_1_', '_2_2_')

886 _Units_ = (Scalar,) * len(_Names_)

887

888 def _validate(self, **unused): # PYCHOK unused

889 '''(INTERNAL) Allow C{_Names_} with leading underscore.

890 '''

891 _NamedTuple._validate(self, underOK=True)

892

893 def __new__(cls, sa, ca, sb, cb, *_more, **name):

894 '''New L{EcefMatrix} matrix.

895

896 @arg sa: C{sin(phi)} (C{float}).

897 @arg ca: C{cos(phi)} (C{float}).

898 @arg sb: C{sin(lambda)} (C{float}).

899 @arg cb: C{cos(lambda)} (C{float}).

900 @arg _more: (INTERNAL) from C{.multiply}.

901

902 @raise EcefError: If B{C{sa}}, B{C{ca}}, B{C{sb}} or

903 B{C{cb}} outside M{[-1.0, +1.0]}.

904 '''

905 t = sa, ca, sb, cb

906 if _more: # all 9 matrix elements ...

907 t += _more # ... from .multiply

908

909 elif max(map(fabs, t)) > _1_0:

910 raise EcefError(unstr(EcefMatrix, *t))

911

912 else: # build matrix from the following quaternion operations

913 # qrot(lam, [0,0,1]) * qrot(phi, [0,-1,0]) * [1,1,1,1]/2

914 # or

915 # qrot(pi/2 + lam, [0,0,1]) * qrot(-pi/2 + phi, [-1,0,0])

916 # where

917 # qrot(t,v) = [cos(t/2), sin(t/2)*v[1], sin(t/2)*v[2], sin(t/2)*v[3]]

918

919 # Local X axis (East) in geocentric coords

920 # M[0] = -slam; M[3] = clam; M[6] = 0;

921 # Local Y axis (North) in geocentric coords

922 # M[1] = -clam * sphi; M[4] = -slam * sphi; M[7] = cphi;

923 # Local Z axis (Up) in geocentric coords

924 # M[2] = clam * cphi; M[5] = slam * cphi; M[8] = sphi;

925 t = (-sb, -cb * sa, cb * ca,

926 cb, -sb * sa, sb * ca,

927 _0_0, ca, sa)

928

929 return _NamedTuple.__new__(cls, *t, **name)

930

931 def column(self, column):

932 '''Get this matrix' B{C{column}} 0, 1 or 2 as C{3-tuple}.

933 '''

934 if 0 <= column < 3:

935 return self[column::3]

936 raise _IndexError(column=column)

937

938 @property_RO

939 def _columns(self):

940 for c in range(3):

941 yield self[c::3]

942

943 def copy(self, **unused): # PYCHOK signature

944 '''Make a shallow or deep copy of this instance.

945

946 @return: The copy (C{This class} or subclass thereof).

947 '''

948 return self.classof(*self)

949

950 __copy__ = __deepcopy__ = copy

951

952 @Property_RO

953 def matrix3(self):

954 '''Get this matrix' rows (C{3-tuple} of 3 C{3-tuple}s).

955 '''

956 return tuple(self._rows)

957

958 @Property_RO

959 def matrixTransposed3(self):

960 '''Get this matrix' I{Transposed} rows (C{3-tuple} of 3 C{3-tuple}s).

961 '''

962 return tuple(self._columns)

963

964 def multiply(self, other):

965 '''Matrix multiply M{M0' ⋅ M} this matrix I{Transposed} with an other matrix.

966

967 @arg other: The other matrix (L{EcefMatrix}).

968

969 @return: The matrix product (L{EcefMatrix}).

970

971 @raise TypeError: If B{C{other}} is not an L{EcefMatrix}.

972 '''

973 _xinstanceof(EcefMatrix, other=other)

974 # like LocalCartesian.MatrixMultiply, C{self.matrixTransposed3 X other.matrix3}

975 # <https://GeographicLib.SourceForge.io/C++/doc/LocalCartesian_8cpp_source.html>

976 X = (_fdotf(t, *c) for t in self._columns for c in other._columns)

977 return _xnamed(EcefMatrix(*X), typename(EcefMatrix.multiply))

978

979 def rotate(self, xyz, *xyz0):

980 '''Forward rotation M{M0' ⋅ ([x, y, z] - [x0, y0, z0])'}.

981

982 @arg xyz: Local C{(x, y, z)} coordinates (C{3-tuple}).

983 @arg xyz0: Optional, local C{(x0, y0, z0)} origin (C{3-tuple}).

984

985 @return: Rotated C{(x, y, z)} location (C{3-tuple}).

986

987 @raise LenError: Unequal C{len(B{xyz})} and C{len(B{xyz0})}.

988 '''

989 if xyz0:

990 if len(xyz0) != len(xyz):

991 raise LenError(self.rotate, xyz0=len(xyz0), xyz=len(xyz))

992 xyz = tuple(s - s0 for s, s0 in zip(xyz, xyz0))

993

994 # x' = M[0] * x + M[3] * y + M[6] * z

995 # y' = M[1] * x + M[4] * y + M[7] * z

996 # z' = M[2] * x + M[5] * y + M[8] * z

997 return tuple(_fdotf(xyz, *c) for c in self._columns)

998

999 def row(self, row):

1000 '''Get this matrix' B{C{row}} 0, 1 or 2 as C{3-tuple}.

1001 '''

1002 if 0 <= row < 3:

1003 r = row * 3

1004 return self[r:r+3]

1005 raise _IndexError(row=row)

1006

1007 @property_RO

1008 def _rows(self):

1009 for r in (0, 3, 6):

1010 yield self[r:r+3]

1011

1012 def unrotate(self, xyz, *xyz0):

1013 '''Inverse rotation M{[x0, y0, z0] + M0 ⋅ [x,y,z]'}.

1014

1015 @arg xyz: Local C{(x, y, z)} coordinates (C{3-tuple}).

1016 @arg xyz0: Optional, local C{(x0, y0, z0)} origin (C{3-tuple}).

1017

1018 @return: Unrotated C{(x, y, z)} location (C{3-tuple}).

1019

1020 @raise LenError: Unequal C{len(B{xyz})} and C{len(B{xyz0})}.

1021 '''

1022 if xyz0:

1023 if len(xyz0) != len(xyz):

1024 raise LenError(self.unrotate, xyz0=len(xyz0), xyz=len(xyz))

1025 _xyz = _1_0_1T + xyz

1026 # x' = x0 + M[0] * x + M[1] * y + M[2] * z

1027 # y' = y0 + M[3] * x + M[4] * y + M[5] * z

1028 # z' = z0 + M[6] * x + M[7] * y + M[8] * z

1029 xyz_ = (_fdotf(_xyz, s, *r) for s, r in zip(xyz0, self._rows))

1030 else:

1031 # x' = M[0] * x + M[1] * y + M[2] * z

1032 # y' = M[3] * x + M[4] * y + M[5] * z

1033 # z' = M[6] * x + M[7] * y + M[8] * z

1034 xyz_ = (_fdotf(xyz, *r) for r in self._rows)

1035 return tuple(xyz_)

1036

1037

1038class Ecef9Tuple(_NamedTuple, _EcefLocal):

1039 '''9-Tuple C{(x, y, z, lat, lon, height, C, M, datum)} with I{geocentric} C{x},

1040 C{y} and C{z} plus I{geodetic} C{lat}, C{lon} and C{height}, case C{C} and

1041 optionally, rotation matrix C{M} (L{EcefMatrix}) and C{datum}, with C{lat}

1042 and C{lon} in C{degrees} and C{x}, C{y}, C{z} and C{height} in C{meter},

1043 conventionally. Case C{C=1} means normal, C{C=2} near polar and C{C=3}

1044 equatorial latitude and C{C=4} height exceeds C{heightMax}.

1045 '''

1046 _Names_ = (_x_, _y_, _z_, _lat_, _lon_, _height_, _C_, _M_, _datum_)

1047 _Units_ = ( Meter, Meter, Meter, Lat, Lon, Height, Int, _Pass, _Pass)

1048

1049 @property_ROver

1050 def _CartesianBase(self):

1051 '''(INTERNAL) Get class C{CartesianBase}, I{once}.

1052 '''

1053 return _MODS.cartesianBase.CartesianBase # overwrite property_ROver

1054

1055 @deprecated_method

1056 def convertDatum(self, datum2): # for backward compatibility

1057 '''DEPRECATED, use method L{toDatum}.'''

1058 return self.toDatum(datum2)

1059

1060 @property_RO

1061 def _ecef9(self): # in ._EcefLocal._Ltp_ecef2local

1062 return self

1063

1064 @property_RO

1065 def ellipsoid(self):

1066 '''Get the ellipsoid (L{Ellipsoid}).

1067 '''

1068 return (self.datum or _WGS84).ellipsoid

1069

1070 @Property_RO

1071 def lam(self):

1072 '''Get the longitude in C{radians} (C{float}).

1073 '''

1074 return self.philam.lam

1075

1076 @Property_RO

1077 def lamVermeille(self):

1078 '''Get the longitude in C{radians} M{[-PI*3/2..+PI*3/2]} after U{Vermeille

1079 <https://Search.ProQuest.com/docview/639493848>} (2004), page 95.

1080

1081 @see: U{Karney<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>},

1082 U{Vermeille<https://Search.ProQuest.com/docview/847292978>} 2011, pp 112-113, 116

1083 and U{Featherstone, et.al.<https://Search.ProQuest.com/docview/872827242>}, page 7.

1084 '''

1085 x, y = self.x, self.y

1086 a = fabs(y)

1087 if a > EPS0:

1088 r = PI_2 - atan2(x, hypot(x, a) + a) * _2_0

1089 if y < 0:

1090 r = -r

1091 else: # y == 0

1092 r = PI if x < 0 else _0_0

1093 return Lam(Vermeille=r)

1094

1095 @Property_RO

1096 def latlon(self):

1097 '''Get the lat-, longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).

1098 '''

1099 return LatLon2Tuple(self.lat, self.lon, name=self.name)

1100

1101 @Property_RO

1102 def latlonheight(self):

1103 '''Get the lat-, longitude in C{degrees} and height (L{LatLon3Tuple}C{(lat, lon, height)}).

1104 '''

1105 return self.latlon.to3Tuple(self.height)

1106

1107 @Property_RO

1108 def latlonheightdatum(self):

1109 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}).

1110 '''

1111 return self.latlonheight.to4Tuple(self.datum)

1112

1113 @Property_RO

1114 def latlonVermeille(self):

1115 '''Get the latitude and I{Vermeille} longitude in C{degrees [-225..+225]} (L{LatLon2Tuple}C{(lat, lon)}).

1116

1117 @see: Property C{lonVermeille}.

1118 '''

1119 return LatLon2Tuple(self.lat, self.lonVermeille, name=self.name)

1120

1121 @Property_RO

1122 def lonVermeille(self):

1123 '''Get the longitude in C{degrees [-225..+225]} after U{Vermeille

1124 <https://Search.ProQuest.com/docview/639493848>} 2004, page 95.

1125

1126 @see: Property C{lamVermeille}.

1127 '''

1128 return Lon(Vermeille=degrees(self.lamVermeille))

1129

1130 @Property_RO

1131 def Mx(self):

1132 '''Compute rotation matrix (L{EcefMatrix}), seperate from C{M}.

1133 '''

1134 sa, ca, sb, cb, _, _, _ = _norm7(self.y, self.x, self.z, self.ellipsoid)

1135 return EcefMatrix(sa, ca, sb, cb, name=self.name)

1136

1137 @Property_RO

1138 def phi(self):

1139 '''Get the latitude in C{radians} (C{float}).

1140 '''

1141 return self.philam.phi

1142

1143 @Property_RO

1144 def philam(self):

1145 '''Get the lat-, longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).

1146 '''

1147 return PhiLam2Tuple(radians(self.lat), radians(self.lon), name=self.name)

1148

1149 @Property_RO

1150 def philamheight(self):

1151 '''Get the lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).

1152 '''

1153 return self.philam.to3Tuple(self.height)

1154

1155 @Property_RO

1156 def philamheightdatum(self):

1157 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLam4Tuple}C{(phi, lam, height, datum)}).

1158 '''

1159 return self.philamheight.to4Tuple(self.datum)

1160

1161 @Property_RO

1162 def philamVermeille(self):

1163 '''Get the latitude and I{Vermeille} longitude in C{radians [-PI*3/2..+PI*3/2]} (L{PhiLam2Tuple}C{(phi, lam)}).

1164

1165 @see: Property C{lamVermeille}.

1166 '''

1167 return PhiLam2Tuple(radians(self.lat), self.lamVermeille, name=self.name)

1168

1169 phiVermeille = phi

1170

1171 def toCartesian(self, Cartesian=None, **Cartesian_kwds):

1172 '''Return the geocentric C{(x, y, z)} coordinates as an ellipsoidal or spherical

1173 C{Cartesian}.

1174

1175 @kwarg Cartesian: Optional class to return C{(x, y, z)} (L{ellipsoidalKarney.Cartesian},

1176 L{ellipsoidalNvector.Cartesian}, L{ellipsoidalVincenty.Cartesian},

1177 L{sphericalNvector.Cartesian} or L{sphericalTrigonometry.Cartesian})

1178 or C{None}.

1179 @kwarg Cartesian_kwds: Optionally, additional B{C{Cartesian}} keyword arguments, ignored

1180 if C{B{Cartesian} is None}.

1181

1182 @return: A B{C{Cartesian}} instance or a L{Vector4Tuple}C{(x, y, z, h)} if C{B{Cartesian}

1183 is None}.

1184

1185 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}} item.

1186 '''

1187 if _isin(Cartesian, None, Vector4Tuple):

1188 r = self.xyzh

1189 elif Cartesian is Vector3Tuple:

1190 r = self.xyz

1191 else:

1192 _xsubclassof(self._CartesianBase, Cartesian=Cartesian)

1193 r = Cartesian(self, **_name1__(Cartesian_kwds, _or_nameof=self))

1194 return r

1195

1196 def toDatum(self, datum2, **name):

1197 '''Convert this C{Ecef9Tuple} to an other datum.

1198

1199 @arg datum2: Datum to convert I{to} (L{Datum}).

1200 @kwarg name: Optional C{B{name}=NN} (C{str}).

1201

1202 @return: The converted 9-Tuple (C{Ecef9Tuple}).

1203

1204 @raise TypeError: The B{C{datum2}} is not a L{Datum}.

1205 '''

1206 n = _name__(name, _or_nameof=self)

1207 if _isin(self.datum, None, datum2): # PYCHOK _Names_

1208 r = self.copy(name=n)

1209 else:

1210 c = self._CartesianBase(self, datum=self.datum, name=n) # PYCHOK _Names_

1211 # c.toLatLon converts datum, x, y, z, lat, lon, etc.

1212 # and returns another Ecef9Tuple iff LatLon is None

1213 r = c.toLatLon(datum=datum2, LatLon=None)

1214 return r

1215

1216 def toLatLon(self, LatLon=None, **LatLon_kwds):

1217 '''Return the geodetic C{(lat, lon, height[, datum])} coordinates.

1218

1219 @kwarg LatLon: Optional class to return C{(lat, lon, height[, datum])} or C{None}.

1220 @kwarg LatLon_kwds: Optional B{C{height}}, B{C{datum}} and other B{C{LatLon}}

1221 keyword arguments.

1222

1223 @return: A B{C{LatLon}} instance or if C{B{LatLon} is None}, a L{LatLon4Tuple}C{(lat,

1224 lon, height, datum)} or L{LatLon3Tuple}C{(lat, lon, height)} if C{datum} is

1225 specified or not.

1226

1227 @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}} item.

1228 '''

1229 lat, lon, D = self.lat, self.lon, self.datum # PYCHOK Ecef9Tuple

1230 kwds = _name1__(LatLon_kwds, _or_nameof=self)

1231 kwds = _xkwds(kwds, height=self.height, datum=D) # PYCHOK Ecef9Tuple

1232 d = kwds.get(_datum_, LatLon)

1233 if LatLon is None:

1234 r = LatLon3Tuple(lat, lon, kwds[_height_], name=kwds[_name_])

1235 if d is not None:

1236 # assert d is not LatLon

1237 r = r.to4Tuple(d) # checks type(d)

1238 else:

1239 if d is None:

1240 _ = kwds.pop(_datum_) # remove None datum

1241 r = LatLon(lat, lon, **kwds)

1242 _xdatum(_xattr(r, datum=D), D)

1243 return r

1244

1245 def toVector(self, Vector=None, **Vector_kwds):

1246 '''Return these geocentric C{(x, y, z)} coordinates as vector.

1247

1248 @kwarg Vector: Optional vector class to return C{(x, y, z)} or C{None}.

1249 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword arguments,

1250 ignored if C{B{Vector} is None}.

1251

1252 @return: A B{C{Vector}} instance or a L{Vector3Tuple}C{(x, y, z)} if

1253 C{B{Vector} is None}.

1254

1255 @raise TypeError: Invalid B{C{Vector}} or B{C{Vector_kwds}} item.

1256

1257 @see: Propertes C{xyz} and C{xyzh}

1258 '''

1259 return self.xyz if Vector is None else Vector(

1260 *self.xyz, **_name1__(Vector_kwds, _or_nameof=self)) # PYCHOK Ecef9Tuple

1261

1262# def _T_x_M(self, T):

1263# '''(INTERNAL) Update M{self.M = T.multiply(self.M)}.

1264# '''

1265# return self.dup(M=T.multiply(self.M))

1266

1267 @Property_RO

1268 def xyz(self):

1269 '''Get the geocentric C{(x, y, z)} coordinates (L{Vector3Tuple}C{(x, y, z)}).

1270 '''

1271 return Vector3Tuple(self.x, self.y, self.z, name=self.name)

1272

1273 @Property_RO

1274 def xyzh(self):

1275 '''Get the geocentric C{(x, y, z)} coordinates and C{height} (L{Vector4Tuple}C{(x, y, z, h)})

1276 '''

1277 return self.xyz.to4Tuple(self.height)

1278

1279

1280def _4Ecef(this, Ecef): # in .datums.Datum.ecef, .ellipsoids.Ellipsoid.ecef

1281 '''Return an ECEF converter for C{this} L{Datum} or L{Ellipsoid}.

1282 '''

1283 if Ecef is None:

1284 Ecef = EcefKarney

1285 else:

1286 _xinstanceof(*_Ecefs, Ecef=Ecef)

1287 return Ecef(this, name=this.name)

1288

1289

1290def _llhn4(latlonh, lon, height, suffix=NN, Error=EcefError, **name): # in .ltp

1291 '''(INTERNAL) Get a C{(lat, lon, h, name)} 4-tuple.

1292 '''

1293 try:

1294 lat, lon = latlonh.lat, latlonh.lon

1295 h = _xattr(latlonh, height=_xattr(latlonh, h=height))

1296 n = _name__(name, _or_nameof=latlonh) # == latlonh._name__(name)

1297 except AttributeError:

1298 lat, h, n = latlonh, height, _name__(**name)

1299 try:

1300 return Lat(lat), Lon(lon), Height(h), n

1301 except (TypeError, ValueError) as x:

1302 t = _lat_, _lon_, _height_

1303 if suffix:

1304 t = (_ + suffix for _ in t)

1305 d = dict(zip(t, (lat, lon, h)))

1306 raise Error(cause=x, **d)

1307

1308

1309def _norm3(y, x, eps=0):

1310 '''(INTERNAL) Return C{y, x, h} normalized.

1311 '''

1312 h = hypot(y, x) # EPS0, EPS_2

1313 return (y / h, x / h, h) if h > eps else (_0_0, _1_0, h)

1314

1315

1316def _norm7(y, x, z=0, E=_EWGS84):

1317 '''(INTERNAL) Return C{phi, lam, p, h, C}.

1318 '''

1319 sb, cb, p = _norm3(y, x) # lam, distance to polar axis

1320 sa, ca, h = _norm3(z, p) # phi, distance to earth center

1321 if h > E._heightMax:

1322 # We are really far away (> 12M light years). Treat the earth

1323 # as a point and h above as an acceptable approximation to the

1324 # height. This avoids overflow, e.g., in the computation of d

1325 # below. It's possible that h has overflowed to INF, that's OK.

1326 # Treat finite x, y, but R overflows to +INF by scaling by 2.

1327 sb, cb, p = _norm3(y * _0_5, x * _0_5)

1328 sa, ca, _ = _norm3(z * _0_5, p)

1329 C = 4

1330 else:

1331 C = 0

1332 return sa, ca, sb, cb, p, h, C

1333

1334

1335def _xEcef(Ecef): # PYCHOK .latlonBase

1336 '''(INTERNAL) Validate B{C{Ecef}} I{class}.

1337 '''

1338 if issubclassof(Ecef, _EcefBase):

1339 return Ecef

1340 raise _TypesError(_Ecef_, Ecef, *_Ecefs)

1341

1342

1343# kwd lon00 unused but will throw a TypeError if misspelled, etc.

1344def _xyzn4(xyz, y, z, Types, Error=EcefError, lon00=0, # PYCHOK unused

1345 _xyz_y_z_names=_xyz_y_z, **name): # in .ltp

1346 '''(INTERNAL) Get an C{(x, y, z, name)} 4-tuple.

1347 '''

1348 try:

1349 n = _name__(name, _or_nameof=xyz) # == xyz._name__(name)

1350 try:

1351 t = xyz.x, xyz.y, xyz.z, n

1352 if not isinstance(xyz, Types):

1353 raise _TypesError(_xyz_y_z_names[0], xyz, *Types)

1354 except AttributeError:

1355 t = map1(float, xyz, y, z) + (n,)

1356 except (TypeError, ValueError) as x:

1357 d = dict(zip(_xyz_y_z_names, (xyz, y, z)))

1358 raise Error(cause=x, **d)

1359 return t

1360# assert _xyz_y_z == _args_kwds_names(_xyzn4)[:3]

1361

1362

1363_Ecefs = tuple(_ for _ in locals().values()

1364 if issubclassof(_, _EcefBase) and

1365 _ is not _EcefBase)

1366__all__ += _ALL_DOCS(_EcefBase)

1367

1368# **) MIT License

1369#

1371#

1372# Permission is hereby granted, free of charge, to any person obtaining a

1373# copy of this software and associated documentation files (the "Software"),

1374# to deal in the Software without restriction, including without limitation

1375# the rights to use, copy, modify, merge, publish, distribute, sublicense,

1376# and/or sell copies of the Software, and to permit persons to whom the

1377# Software is furnished to do so, subject to the following conditions:

1378#

1379# The above copyright notice and this permission notice shall be included

1380# in all copies or substantial portions of the Software.

1381#

1382# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

1383# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

1384# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

1385# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

1386# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

1387# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

1388# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pygeodesy / ecef.py: 95%

508 statements