Coverage for pygeodesy/ecef.py: 96%

432 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-06-07 08:37 -0400

1 

2# -*- coding: utf-8 -*- 

3 

4u'''I{Geocentric} Earth-Centered, Earth-Fixed (ECEF) coordinates. 

5 

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. 

16 

17Following is a copy of I{Karney}'s U{Detailed Description 

18<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Geocentric.html>}. 

19 

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

23 

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°. 

26 

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

30 

31Forward conversion from geodetic to geocentric (ECEF) coordinates is straightforward. 

32 

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. 

36 

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. 

44 

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. 

49 

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>}, 

52for conversion between geodetic and I{local cartesian} cordinates in a I{local tangent plane} as 

53opposed to I{geocentric} (ECEF) ones. 

54''' 

55 

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, \ 

60 _3_0, _4_0, _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 _xattr, _xdatum, _xkwds 

66from pygeodesy.fmath import cbrt, fdot, hypot, hypot1, hypot2_ 

67from pygeodesy.fsums import Fsum, fsumf_ 

68from pygeodesy.interns import NN, _a_, _C_, _datum_, _ellipsoid_, _f_, _height_, \ 

69 _lat_, _lon_, _M_, _name_, _singular_, _SPACE_, \ 

70 _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_ 

80 

81from math import asin, atan2, cos, degrees, fabs, radians, sqrt 

82 

83__all__ = _ALL_LAZY.ecef 

84__version__ = '23.05.23' 

85 

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] 

90 

91 

92class EcefError(_ValueError): 

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

94 ''' 

95 pass 

96 

97 

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 = latlonh.lat 

103 lon = latlonh.lon 

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

105 n = _xattr(latlonh, name=NN) 

106 except AttributeError: 

107 lat, h, n = latlonh, height, NN 

108 

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 

118 

119 

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, _xattr(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,) 

131 

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 

136 

137# assert _xyz_y_z == _xargs_names(_xyzn4)[:3] 

138 

139 

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 

146 

147 def __init__(self, a_ellipsoid, f=None, name=NN): 

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

149 

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

157 

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_ 

173 

174 d = _ellipsoidal_datum(E, name=name) 

175 E = d.ellipsoid 

176 if E.a < EPS or E.f > EPS1: 

177 raise ValueError # _invalid_ 

178 

179 except (TypeError, ValueError) as x: 

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

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

182 

183 self._datum = d 

184 self._E = E 

185 

186 def __eq__(self, other): 

187 '''Compare this and an other Ecef. 

188 

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

190 

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) 

195 

196 @Property_RO 

197 def equatoradius(self): 

198 '''Get the I{equatorial} radius, semi-axis (C{meter}). 

199 ''' 

200 return self.ellipsoid.a 

201 

202 a = equatorialRadius = equatoradius # Karney property 

203 

204 @Property_RO 

205 def datum(self): 

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

207 ''' 

208 return self._datum 

209 

210 @Property_RO 

211 def ellipsoid(self): 

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

213 ''' 

214 return self._E 

215 

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 

222 

223 f = flattening 

224 

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 if _philam: # lat, lon in radians 

229 sa, ca, sb, cb = sincos2_(lat, lon) 

230 lat = Lat(degrees90( lat)) 

231 lon = Lon(degrees180(lon)) 

232 else: 

233 sa, ca, sb, cb = sincos2d_(lat, lon) 

234 

235 E = self.ellipsoid 

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

237 z = (h + n * E.e21) * sa 

238 x = (h + n) * ca 

239 

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

241 return Ecef9Tuple(x * cb, x * sb, z, lat, lon, h, 

242 0, m, self.datum, 

243 name=name or self.name) 

244 

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

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

247 

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

249 latitude (C{degrees}). 

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

251 (C{degrees}). 

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

253 the surface of the ellipsoid. 

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

255 @kwarg name: Optional name (C{str}). 

256 

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

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

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

260 and C{datum} if available. 

261 

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

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

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

265 

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

267 and avoid double angle conversions. 

268 ''' 

269 llhn = _llhn4(latlonh, lon, height, name=name) 

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

271 

272 def forward_(self, phi, lam, height=0, M=False, name=NN): 

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

274 in I{radians}. 

275 

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

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

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

279 the surface of the ellipsoid. 

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

281 @kwarg name: Optional name (C{str}). 

282 

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

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

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

286 

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

288 ''' 

289 try: # like function C{_llhn4} above 

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

291 except (TypeError, ValueError) as x: 

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

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

294 

295 @property_RO 

296 def _Geocentrics(self): 

297 '''(INTERNAL) Valid geocentric classes. 

298 ''' 

299 t = Ecef9Tuple, _MODS.cartesianBase.CartesianBase 

300 _EcefBase._Geocentrics = t # overwrite the property 

301 return t 

302 

303 @Property_RO 

304 def _isYou(self): 

305 '''(INTERNAL) Is this an C{EcefYou}?. 

306 ''' 

307 return isinstance(self, EcefYou) 

308 

309 def _Matrix(self, sa, ca, sb, cb): 

310 '''Creation a rotation matrix. 

311 

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

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

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

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

316 

317 @return: An L{EcefMatrix}. 

318 ''' 

319 return self._xnamed(EcefMatrix(sa, ca, sb, cb)) 

320 

321 def reverse(self, xyz, y=None, z=None, M=False, name=NN): # PYCHOK no cover 

322 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 

323 ''' 

324 notOverloaded(self, xyz, y=y, z=z, M=M, name=name) 

325 

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

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

328 

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

330 

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

332 ''' 

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

334 

335 

336class EcefFarrell21(_EcefBase): 

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

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

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

340 ''' 

341 

342 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M 

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

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

345 page 29. 

346 

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

348 coordinate (C{meter}). 

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

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

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

352 @kwarg name: Optional name (C{str}). 

353 

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

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

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

357 if available. 

358 

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

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

361 zero division error. 

362 

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

364 ''' 

365 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name) 

366 

367 E = self.ellipsoid 

368 a = E.a 

369 a2 = E.a2 

370 b2 = E.b2 

371 e2 = E.e2 

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

373 e4 = E.e4 

374 

375 try: # names as page 29 

376 z2 = z**2 

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

378 

379 p = hypot(x, y) 

380 p2 = p**2 

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

382 F = b2 * z2 * 54 

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

384 s = cbrt(_1_0 + c + sqrt(c**2 + c * 2)) 

385 P = F / (_3_0 * (fsumf_(_1_0, s, _1_0 / s) * G)**2) 

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

387 Q1 = Q + _1_0 

388 r0 = P * e2 * p / Q1 - sqrt(fsumf_(a2 * (Q1 / Q) * _0_5, 

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

390 -P * p2 * _0_5)) 

391 r = p + e2 * r0 

392 v = b2 / (a * sqrt(r**2 + ez)) 

393 

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

395 t = atan2((e2_ * v + _1_0) * z, p) 

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

397 

398 except (ValueError, ZeroDivisionError) as e: 

399 raise EcefError(x=x, y=y, z=z, cause=e) 

400 

401 return Ecef9Tuple(x, y, z, degrees90(t), atan2d(y, x), h, 

402 1, None, self.datum, 

403 name=name or self.name) 

404 

405 

406class EcefFarrell22(_EcefBase): 

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

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

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

410 ''' 

411 

412 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M 

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

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

415 page 30. 

416 

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

418 coordinate (C{meter}). 

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

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

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

422 @kwarg name: Optional name (C{str}). 

423 

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

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

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

427 if available. 

428 

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

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

431 zero division error. 

432 

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

434 ''' 

435 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name) 

436 

437 E = self.ellipsoid 

438 a = E.a 

439 b = E.b 

440 

441 try: # see EcefVeness.reverse 

442 p = hypot(x, y) 

443 s, c = sincos2(atan2(z * a, p * b)) 

444 

445 t = atan2(z + E.e22 * b * s**3, 

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

447 s, c = sincos2(t) 

448 if c: # E.roc1_(s) = E.a / sqrt(1 - E.e2 * s**2) 

449 h = p / c - (E.roc1_(s) if s else a) 

450 else: # polar 

451 h = fabs(z) - b 

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

453 

454 except (ValueError, ZeroDivisionError) as e: 

455 raise EcefError(x=x, y=y, z=z, cause=e) 

456 

457 return Ecef9Tuple(x, y, z, degrees90(t), atan2d(y, x), h, 

458 1, None, self.datum, 

459 name=name or self.name) 

460 

461 

462class EcefKarney(_EcefBase): 

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

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

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

466 

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

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

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

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

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

472 the rotation matrix. 

473 ''' 

474 

475 @Property_RO 

476 def hmax(self): 

477 '''Get the distance or height limit (C{meter}, conventionally). 

478 ''' 

479 return self.equatoradius / EPS_2 # self.equatoradius * _2_EPS, 12M lighyears 

480 

481 def reverse(self, xyz, y=None, z=None, M=False, name=NN): 

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

483 

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

485 coordinate (C{meter}). 

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

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

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

489 @kwarg name: Optional name (C{str}). 

490 

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

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

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

494 C{datum} if available. 

495 

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

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

498 

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

500 C{height} is returned, i.e., C{(lat, lon)} corresponds to the closest 

501 point on the ellipsoid. If there are still multiple solutions with 

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

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

504 different longitudes (applies only if C{x} = C{y} = 0) then C{lon} = 0 

505 is returned. The returned C{height} value is not below M{−E.a * (1 − 

506 E.e2) / sqrt(1 − E.e2 * sin(lat)**2)}. The returned C{lon} is in the 

507 range [−180°, 180°]. Like C{forward} above, M{v1 = Transpose(M) ⋅ v0}. 

508 ''' 

509 def _norm3(y, x): 

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

511 return (y / h, x / h, h) if h > 0 else (_0_0, _1_0, h) 

512 

513 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name) 

514 

515 E = self.ellipsoid 

516 

517 sb, cb, R = _norm3(y, x) 

518 h = hypot(R, z) # distance to earth center 

519 if h > self.hmax: # PYCHOK no cover 

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

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

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

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

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

525 sb, cb, R = _norm3(y * _0_5, x * _0_5) 

526 sa, ca, _ = _norm3(z * _0_5, R) 

527 C = 1 

528 

529 elif E.e4: # E.isEllipsoidal 

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

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

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

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

534 if E.isProlate: 

535 p, q = q, p 

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

537 e = E.e4 * q 

538 if e or r > 0: 

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

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

541 s = e * p / _4_0 # s = r^3 * s 

542 u = r = r / _6_0 

543 r2 = r**2 

544 r3 = r * r2 

545 t3 = s + r3 

546 d = s * (r3 + t3) 

547 if d < 0: 

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

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

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

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

552 else: 

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

554 # minimizes loss of precision due to cancellation. The 

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

556 # in definition of u. 

557 if d > 0: 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

572 if E.isProlate: 

573 k1 -= E.e2 

574 else: 

575 k2 += E.e2 

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

577 h *= k1 - E.e21 

578 C = 2 

579 

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

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

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

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

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

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

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

587 q = E.e4 - p 

588 if E.isProlate: 

589 p, q = q, p 

590 e = E.a 

591 else: 

592 e = E.b2_a 

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

594 if z < 0: 

595 sa = neg(sa) # for tiny negative z, not for prolate 

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

597 C = 3 

598 

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

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

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

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

603 h -= E.a 

604 C = 4 

605 

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

607 return Ecef9Tuple(x, y, z, atan2d(sa, ca), 

608 atan2d(sb, cb), h, 

609 C, m, self.datum, 

610 name=name or self.name) 

611 

612 

613class EcefSudano(_EcefBase): 

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

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

616 3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}. 

617 ''' 

618 _tol = EPS2 

619 

620 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M 

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

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

623 3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}. 

624 

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

626 coordinate (C{meter}). 

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

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

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

630 @kwarg name: Optional name (C{str}). 

631 

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

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

634 iteration C{C}, C{M=None} always and C{datum} if available. 

635 

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

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

638 ''' 

639 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name) 

640 

641 E = self.ellipsoid 

642 e = E.e2 * E.a 

643 h = hypot(x, y) # Rh 

644 d = e - h 

645 

646 a = atan2(z, h * E.e21) 

647 sa, ca = sincos2(fabs(a)) 

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

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

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

651 # = ca**2 * (z * ca + E.e2 * E.a * sa - h * sa) 

652 # = ca**2 * (z * ca + (E.e2 * E.a - h) * sa) 

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

654 # = ca**2 * (E.e2 * E.a / E.e2s2(sa) - h / ca**2) 

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

656 # (E.e2 * E.a / E.e2s2(sa) - h / ca**2) 

657 tol = self.tolerance 

658 _S2_ = Fsum(sa).fsum2_ 

659 for C in range(1, _TRIPS): 

660 ca2 = _1_0 - sa**2 

661 if ca2 < EPS_2: # PYCHOK no cover 

662 ca = _0_0 

663 break 

664 ca = sqrt(ca2) 

665 r = e / E.e2s2(sa) - h / ca2 

666 if fabs(r) < EPS_2: 

667 break 

668 a = None 

669 sa, r = _S2_(-z * ca / r, -d * sa / r) 

670 if fabs(r) < tol: 

671 break 

672 else: 

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

674 raise EcefError(Fmt.no_convergence(r, tol), txt=t) 

675 

676 if a is None: 

677 a = copysign0(asin(sa), z) 

678 h = fsumf_(h * ca, fabs(z * sa), -E.a * E.e2s(sa)) # use Veness', 

679 # since Sudano's Eq (7) doesn't provide the correct height 

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

681 

682 r = Ecef9Tuple(x, y, z, degrees90(a), atan2d(y, x), h, 

683 C, None, self.datum, 

684 name=name or self.name) 

685 r._iteration = C 

686 return r 

687 

688 @property_doc_(''' the convergence tolerance (C{float}).''') 

689 def tolerance(self): 

690 '''Get the convergence tolerance (C{scalar}). 

691 ''' 

692 return self._tol 

693 

694 @tolerance.setter # PYCHOK setter! 

695 def tolerance(self, tol): 

696 '''Set the convergence tolerance (C{scalar}). 

697 

698 @raise EcefError: Non-scalar or invalid B{C{tol}}. 

699 ''' 

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

701 

702 

703class EcefVeness(_EcefBase): 

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

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

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

707 

708 @see: U{I{A Guide to Coordinate Systems in Great Britain}<https:// 

709 www.OrdnanceSurvey.co.UK/documents/resources/guide-coordinate-systems-great-britain.pdf>}, 

710 section I{B) Converting between 3D Cartesian and ellipsoidal 

711 latitude, longitude and height coordinates}. 

712 ''' 

713 

714 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M 

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

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

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

718 

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

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

721 233668213_The_Accuracy_of_Geodetic_Latitude_and_Height_Equations>}, Survey Review, 

722 Vol 28, 218, Oct 1985. 

723 

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

725 coordinate (C{meter}). 

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

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

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

729 @kwarg name: Optional name (C{str}). 

730 

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

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

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

734 

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

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

737 

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

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

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

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

742 Apr 1996, both from Lawrence Livermore National Laboratory (LLNL) and 

743 Sudano, John J, U{I{An exact conversion from an Earth-centered coordinate 

744 system to latitude longitude and altitude}<https://www.ResearchGate.net/ 

745 publication/3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}. 

746 ''' 

747 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name) 

748 

749 E = self.ellipsoid 

750 p = hypot(x, y) # distance from minor axis 

751 r = hypot(p, z) # polar radius 

752 if min(p, r) > EPS0: 

753 b = E.b * E.e22 

754 # parametric latitude (Bowring eqn 17, replaced) 

755 t = (E.b * z) / (E.a * p) * (_1_0 + b / r) 

756 c = _1_0 / hypot1(t) 

757 s = t * c 

758 

759 # geodetic latitude (Bowring eqn 18) 

760 a = atan2(z + b * s**3, 

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

762 

763 # height above ellipsoid (Bowring eqn 7) 

764 sa, ca = sincos2(a) 

765# r = E.a / E.e2s(sa) # length of normal terminated by minor axis 

766# h = p * ca + z * sa - (E.a * E.a / r) 

767 h = fsumf_(p * ca, z * sa, -E.a * E.e2s(sa)) 

768 

769 C, lat, lon = 1, degrees90(a), atan2d(y, x) 

770 

771 # see <https://GIS.StackExchange.com/questions/28446> 

772 elif p > EPS: # lat arbitrarily zero 

773 C, lat, lon, h = 2, _0_0, atan2d(y, x), p - E.a 

774 

775 else: # polar lat, lon arbitrarily zero 

776 C, lat, lon, h = 3, copysign0(_90_0, z), _0_0, fabs(z) - E.b 

777 

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

779 C, None, # M=None 

780 self.datum, name=name or self.name) 

781 

782 

783class EcefYou(_EcefBase): 

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

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

786 

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

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

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

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

791 ''' 

792 

793 def __init__(self, a_ellipsoid, f=None, name=NN): 

794 _EcefBase.__init__(self, a_ellipsoid, f=f, name=name) # inherited documentation 

795 E = self.ellipsoid 

796 if E.isProlate or (E.a2 - E.b2) < 0: 

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

798 

799 def reverse(self, xyz, y=None, z=None, M=None, name=NN): # PYCHOK unused M 

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

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

802 

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

804 coordinate (C{meter}). 

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

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

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

808 @kwarg name: Optional name (C{str}). 

809 

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

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

812 ones C{(x, y, z)}, case C{C=1}, C{M=None} always and C{datum} if 

813 available. 

814 

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

816 B{C{z}} not C{scalar} for C{scalar} B{C{xyz}}. 

817 ''' 

818 x, y, z, name = _xyzn4(xyz, y, z, self._Geocentrics, name=name) 

819 

820 r2 = hypot2_(x, y, z) 

821 

822 E = self.ellipsoid 

823 e2 = E.a2 - E.b2 # == E.e2 * E.a2 

824 if e2 < 0: 

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

826 e = sqrt(e2) # XXX sqrt0(e2)? 

827 

828 q = hypot(x, y) 

829 u = fsumf_(r2, -e2, hypot(r2 - e2, 2 * e * z)) * _0_5 

830 if u > EPS02: 

831 u = sqrt(u) 

832 p = hypot(u, e) 

833 B = atan2(p * z, u * q) # beta0 = atan(p / u * z / q) 

834 sB, cB = sincos2(B) 

835 if cB and sB: 

836 p *= E.a 

837 d = (p / cB - e2 * cB) / sB 

838 if isnon0(d): 

839 B += fsumf_(u * E.b, -p, e2) / d 

840 sB, cB = sincos2(B) 

841 elif u < 0: 

842 raise EcefError(x=x, y=y, z=z, txt=_singular_) 

843 else: 

844 sB, cB = copysign0(_1_0, z), _0_0 

845 

846 h = hypot(z - E.b * sB, q - E.a * cB) 

847 if hypot2_(x, y, z * E.a_b) < E.a2: 

848 h = neg(h) # inside ellipsoid 

849 

850 return Ecef9Tuple(x, y, z, atan2d(E.a * sB, E.b * cB), # atan(E.a_b * tan(B)) 

851 atan2d(y, x), h, 

852 1, None, # C=1, M=None 

853 self.datum, name=name or self.name) 

854 

855 

856class EcefMatrix(_NamedTuple): 

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

858 

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

860 Geographic_coordinate_conversion#From_ECEF_to_ENU>} and 

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

862 ''' 

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

864 '_1_0_', '_1_1_', '_1_2_', 

865 '_2_0_', '_2_1_', '_2_2_') 

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

867 

868 def _validate(self, **_OK): # PYCHOK unused 

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

870 ''' 

871 _NamedTuple._validate(self, _OK=True) 

872 

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

874 '''New L{EcefMatrix} matrix. 

875 

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

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

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

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

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

881 

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

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

884 ''' 

885 t = sa, ca, sb, cb 

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

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

888 

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

890 raise EcefError(unstr(EcefMatrix.__name__, *t)) 

891 

892 else: # build matrix from the following quaternion operations 

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

894 # or 

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

896 # where 

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

898 

899 # Local X axis (East) in geocentric coords 

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

901 # Local Y axis (North) in geocentric coords 

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

903 # Local Z axis (Up) in geocentric coords 

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

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

906 cb, -sb * sa, sb * ca, 

907 _0_0, ca, sa) 

908 

909 return _NamedTuple.__new__(cls, *t) 

910 

911 def column(self, column): 

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

913 ''' 

914 if 0 <= column < 3: 

915 return self[column::3] 

916 raise _IndexError(column=column) 

917 

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

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

920 

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

922 ''' 

923 return self.classof(*self) 

924 

925 __copy__ = __deepcopy__ = copy 

926 

927 @Property_RO 

928 def matrix3(self): 

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

930 ''' 

931 return tuple(map(self.row, range(3))) 

932 

933 @Property_RO 

934 def matrixTransposed3(self): 

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

936 ''' 

937 return tuple(map(self.column, range(3))) 

938 

939 def multiply(self, other): 

940 '''Matrix multiply M{M0' ⋅ M} this matrix I{Transposed} 

941 with an other matrix. 

942 

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

944 

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

946 

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

948 ''' 

949 _xinstanceof(EcefMatrix, other=other) 

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

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

952 # X = (fdot(self.column(r), *other.column(c)) for r in (0,1,2) for c in (0,1,2)) 

953 X = (fdot(self[r::3], *other[c::3]) for r in range(3) for c in range(3)) 

954 return _xnamed(EcefMatrix(*X), EcefMatrix.multiply.__name__) 

955 

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

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

958 

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

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

961 

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

963 

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

965 ''' 

966 if xyz0: 

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

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

969 xyz = tuple(c - c0 for c, c0 in zip(xyz, xyz0)) 

970 

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[0::3]), # .column(0) 

975 fdot(xyz, *self[1::3]), # .column(1) 

976 fdot(xyz, *self[2::3])) # .column(2) 

977 

978 def row(self, row): 

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

980 ''' 

981 if 0 <= row < 3: 

982 r = row * 3 

983 return self[r:r+3] 

984 raise _IndexError(row=row) 

985 

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

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

988 

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

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

991 

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

993 

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

995 ''' 

996 if xyz0: 

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

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

999 _xyz = _1_0_1T + xyz 

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

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

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

1003 xyz_ = (fdot(_xyz, xyz0[0], *self[0:3]), # .row(0) 

1004 fdot(_xyz, xyz0[1], *self[3:6]), # .row(1) 

1005 fdot(_xyz, xyz0[2], *self[6:9])) # .row(2) 

1006 else: 

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

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

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

1010 xyz_ = (fdot(xyz, *self[0:3]), # .row(0) 

1011 fdot(xyz, *self[3:6]), # .row(1) 

1012 fdot(xyz, *self[6:9])) # .row(2) 

1013 return xyz_ 

1014 

1015 

1016class Ecef9Tuple(_NamedTuple): 

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

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

1019 C{C} (see the C{Ecef*.reverse} methods) and optionally, the rotation 

1020 matrix C{M} (L{EcefMatrix}) and C{datum}, with C{lat} and C{lon} in 

1021 C{degrees} and C{x}, C{y}, C{z} and C{height} in C{meter}, conventionally. 

1022 ''' 

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

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

1025 

1026 @property_RO 

1027 def _CartesianBase(self): 

1028 '''(INTERNAL) Get/cache class C{CartesianBase}. 

1029 ''' 

1030 Ecef9Tuple._CartesianBase = C = _MODS.cartesianBase.CartesianBase # overwrite property 

1031 return C 

1032 

1033 @deprecated_method 

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

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

1036 return self.toDatum(datum2) 

1037 

1038 @Property_RO 

1039 def lam(self): 

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

1041 ''' 

1042 return self.philam.lam 

1043 

1044 @Property_RO 

1045 def lamVermeille(self): 

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

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

1048 

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

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

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

1052 ''' 

1053 x, y = self.x, self.y 

1054 if y > EPS0: 

1055 r = _N_2_0 * atan2(x, hypot(y, x) + y) + PI_2 

1056 elif y < -EPS0: 

1057 r = _2_0 * atan2(x, hypot(y, x) - y) - PI_2 

1058 else: # y == 0 

1059 r = PI if x < 0 else _0_0 

1060 return Lam(Vermeille=r) 

1061 

1062 @Property_RO 

1063 def latlon(self): 

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

1065 ''' 

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

1067 

1068 @Property_RO 

1069 def latlonheight(self): 

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

1071 ''' 

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

1073 

1074 @Property_RO 

1075 def latlonheightdatum(self): 

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

1077 ''' 

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

1079 

1080 @Property_RO 

1081 def latlonVermeille(self): 

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

1083 

1084 @see: Property C{lonVermeille}. 

1085 ''' 

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

1087 

1088 @Property_RO 

1089 def lonVermeille(self): 

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

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

1092 

1093 @see: Property C{lamVermeille}. 

1094 ''' 

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

1096 

1097 @Property_RO 

1098 def phi(self): 

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

1100 ''' 

1101 return self.philam.phi 

1102 

1103 @Property_RO 

1104 def philam(self): 

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

1106 ''' 

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

1108 

1109 @Property_RO 

1110 def philamheight(self): 

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

1112 ''' 

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

1114 

1115 @Property_RO 

1116 def philamheightdatum(self): 

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

1118 ''' 

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

1120 

1121 @Property_RO 

1122 def philamVermeille(self): 

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

1124 

1125 @see: Property C{lamVermeille}. 

1126 ''' 

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

1128 

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

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

1131 C{Cartesian}. 

1132 

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

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

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

1136 or C{None}. 

1137 @kwarg Cartesian_kwds: Optional, additional B{C{Cartesian}} keyword arguments, ignored 

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

1139 

1140 @return: A C{B{Cartesian}(x, y, z, **B{Cartesian_kwds})} instance or 

1141 a L{Vector4Tuple}C{(x, y, z, h)} if C{B{Cartesian} is None}. 

1142 

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

1144 ''' 

1145 if Cartesian in (None, Vector4Tuple): 

1146 r = self.xyzh 

1147 elif Cartesian is Vector3Tuple: 

1148 r = self.xyz 

1149 else: 

1150 _xsubclassof(self._CartesianBase, Cartesian=Cartesian) 

1151 r = Cartesian(self, **_xkwds(Cartesian_kwds, name=self.name)) 

1152 return r 

1153 

1154 def toDatum(self, datum2): 

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

1156 

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

1158 

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

1160 

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

1162 ''' 

1163 if self.datum in (None, datum2): # PYCHOK _Names_ 

1164 r = self.copy() 

1165 else: 

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

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

1168 # and returns another Ecef9Tuple iff LatLon is None 

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

1170 return r 

1171 

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

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

1174 

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

1176 or C{None}. 

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

1178 B{C{LatLon}} keyword arguments. 

1179 

1180 @return: An instance of C{B{LatLon}(lat, lon, **B{LatLon_kwds})} 

1181 or if B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, 

1182 height)} respectively L{LatLon4Tuple}C{(lat, lon, height, 

1183 datum)} depending on whether C{datum} is un-/specified. 

1184 

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

1186 ''' 

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

1188 kwds = _xkwds(LatLon_kwds, height=self.height, datum=D, name=self.name) # PYCHOK Ecef9Tuple 

1189 d = kwds.get(_datum_, LatLon) 

1190 if LatLon is None: 

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

1192 if d is not None: 

1193 # assert d is not LatLon 

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

1195 else: 

1196 if d is None: 

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

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

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

1200 return r 

1201 

1202 def toLocal(self, ltp, Xyz=None, **Xyz_kwds): 

1203 '''Convert this geocentric to I{local} C{x}, C{y} and C{z}. 

1204 

1205 @kwarg ltp: The I{local tangent plane} (LTP) to use (L{Ltp}). 

1206 @kwarg Xyz: Optional class to return C{x}, C{y} and C{z} 

1207 (L{XyzLocal}, L{Enu}, L{Ned}) or C{None}. 

1208 @kwarg Xyz_kwds: Optional, additional B{C{Xyz}} keyword 

1209 arguments, ignored if C{B{Xyz} is None}. 

1210 

1211 @return: An B{C{Xyz}} instance or if C{B{Xyz} is None}, 

1212 a L{Local9Tuple}C{(x, y, z, lat, lon, height, 

1213 ltp, ecef, M)} with C{M=None}, always. 

1214 

1215 @raise TypeError: Invalid B{C{ltp}}. 

1216 ''' 

1217 return _MODS.ltp._xLtp(ltp)._ecef2local(self, Xyz, Xyz_kwds) 

1218 

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

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

1221 

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

1223 C{None}. 

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

1225 arguments, ignored if C{B{Vector} is None}. 

1226 

1227 @return: A C{Vector}C{(x, y, z, **Vector_kwds)} instance or a 

1228 L{Vector3Tuple}C{(x, y, z)} if B{C{Vector}} is C{None}. 

1229 

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

1231 ''' 

1232 return self.xyz if Vector is None else self._xnamed( 

1233 Vector(self.x, self.y, self.z, **Vector_kwds)) # PYCHOK Ecef9Tuple 

1234 

1235# def _T_x_M(self, T): 

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

1237# ''' 

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

1239 

1240 @Property_RO 

1241 def xyz(self): 

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

1243 ''' 

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

1245 

1246 @Property_RO 

1247 def xyzh(self): 

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

1249 ''' 

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

1251 

1252 

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

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

1255 ''' 

1256 if Ecef is None: 

1257 Ecef = EcefKarney 

1258 else: 

1259 _xinstanceof(*_Ecefs, Ecef=Ecef) 

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

1261 

1262 

1263def _xEcef(Ecef): # PYCHOK .latlonBase.py 

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

1265 ''' 

1266 if issubclassof(Ecef, _EcefBase): 

1267 return Ecef 

1268 raise _TypesError(_Ecef_, Ecef, *_Ecefs) 

1269 

1270 

1271_Ecefs = (EcefKarney, EcefSudano, EcefVeness, EcefYou, 

1272 EcefFarrell21, EcefFarrell22) 

1273 

1274__all__ += _ALL_DOCS(_EcefBase) 

1275 

1276# **) MIT License 

1277# 

1278# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

1279# 

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

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

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

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

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

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

1286# 

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

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

1289# 

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

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

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

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

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

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

1296# OTHER DEALINGS IN THE SOFTWARE.