1

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

3

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

5

6 Geocentric conversions transcoded from I{Charles Karney}'s C++ class U{Geocentric

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

8 into pure Python class L{EcefKarney}, class L{EcefSudano} based on I{John Sudano}'s

9 U{paper<https://www.ResearchGate.net/publication/

10 3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>},

11 class L{EcefVeness} transcoded from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal,

12 Cartesian<https://www.Movable-Type.co.UK/scripts/geodesy/docs/latlon-ellipsoidal.js.html>}, class L{EcefYou}

13 implementing I{Rey-Jer You}'s U{transformations <https://www.ResearchGate.net/publication/240359424>} and

14 classes 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

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

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

19

20 Convert between geodetic coordinates C{lat}-, C{lon}gitude and height C{h} (measured vertically

21 from the surface of the ellipsoid) to geocentric C{x}, C{y} and C{z} coordinates, also known as

22 I{Earth-Centered, Earth-Fixed} (U{ECEF<https://WikiPedia.org/wiki/ECEF>}).

23

24 The origin of geocentric coordinates is at the center of the earth. The C{z} axis goes thru

25 the North pole, C{lat} = 90°. The C{x} axis goes thru C{lat} = 0°, C{lon} = 0°.

26

27 The I{local (cartesian) origin} is at (C{lat0}, C{lon0}, C{height0}). The I{local} C{x} axis points

28 East, the I{local} C{y} axis points North and the I{local} C{z} axis is normal to the ellipsoid. The

29 plane C{z = -height0} is tangent to the ellipsoid, hence the alternate name I{local tangent plane}.

30

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

32

33 For the reverse transformation we use Hugues Vermeille's U{I{Direct transformation from geocentric

34 coordinates to geodetic coordinates}<https://DOI.org/10.1007/s00190-002-0273-6>}, J. Geodesy

35(2002) 76, page 451-454.

36

37 Several changes have been made to ensure that the method returns accurate results for all finite

38 inputs (even if h is infinite). The changes are described in Appendix B of C. F. F. Karney

39 U{I{Geodesics on an ellipsoid of revolution}<https://ArXiv.org/abs/1102.1215v1>}, Feb. 2011, 85,

40 105-117 (U{preprint<https://ArXiv.org/abs/1102.1215v1>}). Vermeille similarly updated his method

41 in 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

43 coordinates<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>} for more information.

44

45 The errors in these routines are close to round-off. Specifically, for points within 5,000 Km of

46 the surface of the ellipsoid (either inside or outside the ellipsoid), the error is bounded by 7

47 nm (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

51 U{LocalCartesian <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1LocalCartesian.html>},

52 providing conversion to and from I{local} cartesian cordinates in a I{local tangent plane} as

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

54'''

55

56 from pygeodesy.basics import copysign0, isscalar, issubclassof, neg, map1, \

57 _xinstanceof, _xsubclassof

58 from 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!

62 from pygeodesy.datums import a_f2Tuple, _ellipsoidal_datum

63# from pygeodesy.ellipsoids import a_f2Tuple # from .datums

64 from pygeodesy.errors import _IndexError, LenError, _ValueError, _TypesError, \

65 _xattr, _xdatum, _xkwds

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

67 from pygeodesy.fsums import Fsum, fsumf_

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

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

70 _x_, _xyz_, _y_, _z_

71 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS

72 from pygeodesy.named import _NamedBase, _NamedTuple, notOverloaded, _Pass, _xnamed

73 from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \

74 PhiLam2Tuple, Vector3Tuple, Vector4Tuple

75 from pygeodesy.props import deprecated_method, Property_RO, property_RO, property_doc_

76 from pygeodesy.streprs import Fmt, unstr

77 from pygeodesy.units import Height, Int, Lam, Lat, Lon, Meter, Phi, Scalar, Scalar_

78 from pygeodesy.utily import atan2d, degrees90, degrees180, sincos2, sincos2_, \

79 sincos2d_

80

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

82

83__all__ = _ALL_LAZY.ecef

84__version__ = '23.05.15'

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

92 class EcefError(_ValueError):

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

94 '''

95 pass

96

97

98 def _llhn4(latlonh, lon, height, suffix=NN, Error=EcefError, name=NN): # in .ltp.LocalCartesian.forward and -.reset

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

100 '''

101 try:

102 lat, lon = latlonh.lat, latlonh.lon

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

104 n = _xattr(latlonh, name=NN)

105 except AttributeError:

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

107

108 try:

109 llhn = Lat(lat), Lon(lon), Height(h), (name or n)

110 except (TypeError, ValueError) as x:

111 t = _lat_, _lon_, _height_

112 if suffix:

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

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

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

116 return llhn

117

118

119 def _xyzn4(xyz, y, z, Types, Error=EcefError, name=NN, # in .ltp

120 _xyz_y_z_names=_xyz_y_z):

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

122 '''

123 try:

124 try:

125 t = xyz.x, xyz.y, xyz.z, _xattr(xyz, name=name)

126 if not isinstance(xyz, Types):

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

128 except AttributeError:

129 t = map1(float, xyz, y, z) + (name,)

130

131 except (TypeError, ValueError) as x:

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

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

134 return t

135

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

137

138

139 class _EcefBase(_NamedBase):

140 '''(INTERNAL) Base class for L{EcefFarrell21}, L{EcefFarrell22}, L{EcefKarney},

141 L{EcefSudano}, L{EcefVeness} and L{EcefYou}.

142 '''

143 _datum = None

144 _E = None

145

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

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

148

149 @arg a_ellipsoid: A (non-prolate) ellipsoid (L{Ellipsoid}, L{Ellipsoid2},

150 L{Datum} or L{a_f2Tuple}) or C{scalar} ellipsoid's

151 equatorial radius (C{meter}).

152 @kwarg f: C{None} or the ellipsoid flattening (C{scalar}), required

153 for C{scalar} B{C{a_ellipsoid}}, C{B{f}=0} represents a

154 sphere, negative B{C{f}} a prolate ellipsoid.

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

156

157 @raise EcefError: If B{C{a_ellipsoid}} not L{Ellipsoid}, L{Ellipsoid2},

158 L{Datum} or L{a_f2Tuple} or C{scalar} or B{C{f}} not

159 C{scalar} or if C{scalar} B{C{a_ellipsoid}} not positive

160 or B{C{f}} not less than 1.0.

161 '''

162 if name:

163 self.name = name

164 try:

165 E = a_ellipsoid

166 if f is None:

167 pass

168 elif isscalar(E) and isscalar(f):

169 E = a_f2Tuple(E, f)

170 else:

171 raise ValueError # _invalid_

172

173 d = _ellipsoidal_datum(E, name=name)

174 E = d.ellipsoid

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

176 raise ValueError # _invalid_

177

178 except (TypeError, ValueError) as x:

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

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

181

182 self._datum = d

183 self._E = E

184

185 def __eq__(self, other):

186 '''Compare this and an other Ecef.

187

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

189

190 @return: C{True} if equal, C{False} otherwise.

191 '''

192 return other is self or (isinstance(other, self.__class__) and

193 other.ellipsoid == self.ellipsoid)

194

195 @Property_RO

196 def equatoradius(self):

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

198 '''

199 return self.ellipsoid.a

200

201 a = equatorialRadius = equatoradius # Karney property

202

203 @Property_RO

204 def datum(self):

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

206 '''

207 return self._datum

208

209 @Property_RO

210 def ellipsoid(self):

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

212 '''

213 return self._E

214

215 @Property_RO

216 def flattening(self): # Karney property

217 '''Get the I{flattening} (C{float}), M{(a - b) / a}, positive for

218 I{oblate}, negative for I{prolate} or C{0} for I{near-spherical}.

219 '''

220 return self.ellipsoid.f

221

222 f = flattening

223

224 def _forward(self, lat, lon, h, name, M=False, _philam=False): # in .ltp.LocalCartesian.forward and -.reset

225 '''(INTERNAL) Common for all C{Ecef*}.

226 '''

227 E = self.ellipsoid

228

229 if _philam:

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

231 lat = Lat(degrees90( lat))

232 lon = Lon(degrees180(lon))

233 else:

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

235

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 _EcefBase._forward(self, *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

336 class 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

406 class 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

462 class 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 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 disc

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 = 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 disc = s * (r3 + t3)

547 if disc < 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, disc < 0 implies r < 0.

551 u += _2_0 * r * cos(atan2(sqrt(-disc), -t3) / _3_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 disc > 0:

558 t3 += copysign0(sqrt(disc), 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 # t can be zero; but then r2 / t -> 0.

562 if t:

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

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

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

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

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

568 # Need to guard against w going negative due to roundoff in uv - q.

569 w = max(_0_0, E.e2abs * (uv - q) / (_2_0 * v))

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

571 # subtraction. Division by 0 not possible because uv > 0, w >= 0.

572 k1 = k2 = uv / (sqrt(uv + w**2) + w)

573 if E.isProlate:

574 k1 -= E.e2

575 else:

576 k2 += E.e2

577 sa, ca, h = norm3(z / k1, R / k2)

578 h *= k1 - E.e21

579 C = 2

580

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

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

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

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

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

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

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

588 q = E.e4 - p

589 if E.isProlate:

590 p, q = q, p

591 e = E.a

592 else:

593 e = E.b2_a

594 sa, ca, h = norm3(sqrt(q * E._1_e21), sqrt(p))

595 if z < 0:

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

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

598 C = 3

599

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

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

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

603 sa, ca, _ = norm3((z if h else _1_0), R)

604 h -= E.a

605 C = 4

606

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

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

609 atan2d(sb, cb), h,

610 C, m, self.datum,

611 name=name or self.name)

612

613

614 class EcefSudano(_EcefBase):

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

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

617 3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}.

618 '''

619 _tol = EPS2

620

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

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

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

624 3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}.

625

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

627 coordinate (C{meter}).

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

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

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

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

632

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

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

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

636

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

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

639 '''

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

641

642 E = self.ellipsoid

643 e = E.e2 * E.a

644 h = hypot(x, y) # Rh

645 d = e - h

646

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

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

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

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

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

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

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

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

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

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

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

658 tol = self.tolerance

659 _S2_ = Fsum(sa).fsum2_

660 for C in range(1, _TRIPS):

661 ca2 = _1_0 - sa**2

662 if ca2 < EPS_2: # PYCHOK no cover

663 ca = _0_0

664 break

665 ca = sqrt(ca2)

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

667 if fabs(r) < EPS_2:

668 break

669 a = None

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

671 if fabs(r) < tol:

672 break

673 else:

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

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

676

677 if a is None:

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

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

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

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

682

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

684 C, None, self.datum,

685 name=name or self.name)

686 r._iteration = C

687 return r

688

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

690 def tolerance(self):

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

692 '''

693 return self._tol

694

695 @tolerance.setter # PYCHOK setter!

696 def tolerance(self, tol):

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

698

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

700 '''

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

702

703

704 class EcefVeness(_EcefBase):

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

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

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

708

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

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

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

712 latitude, longitude and height coordinates}.

713 '''

714

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

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

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

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

719

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

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

722 233668213_The_Accuracy_of_Geodetic_Latitude_and_Height_Equations>}, Survey Review,

723 Vol 28, 218, Oct 1985.

724

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

726 coordinate (C{meter}).

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

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

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

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

731

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

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

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

735

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

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

738

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

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

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

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

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

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

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

746 publication/3709199_An_exact_conversion_from_an_Earth-centered_coordinate_system_to_latitude_longitude_and_altitude>}.

747 '''

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

749

750 E = self.ellipsoid

751

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

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

754 if min(p, r) > EPS0:

755 b = E.b * E.e22

756 # parametric latitude (Bowring eqn 17, replaced)

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

758 c = _1_0 / hypot1(t)

759 s = t * c

760

761 # geodetic latitude (Bowring eqn 18)

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

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

764

765 # height above ellipsoid (Bowring eqn 7)

766 sa, ca = sincos2(a)

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

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

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

770

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

772

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

774 elif p > EPS: # lat arbitrarily zero

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

776

777 else: # polar lat, lon arbitrarily zero

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

779

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

781 C, None, # M=None

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

783

784

785 class EcefYou(_EcefBase):

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

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

788

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

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

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

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

793 '''

794

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

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

797 E = self.ellipsoid

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

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

800

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

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

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

804

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

806 coordinate (C{meter}).

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

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

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

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

811

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

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

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

815 available.

816

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

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

819 '''

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

821

822 r2 = hypot2_(x, y, z)

823

824 E = self.ellipsoid

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

826 if e2 < 0:

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

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

829

830 q = hypot(x, y)

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

832 if u > EPS02:

833 u = sqrt(u)

834 p = hypot(u, e)

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

836 sB, cB = sincos2(B)

837 if cB and sB:

838 p *= E.a

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

840 if isnon0(d):

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

842 sB, cB = sincos2(B)

843 elif u < 0:

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

845 else:

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

847

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

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

850 h = neg(h) # inside ellipsoid

851

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

853 atan2d(y, x), h,

854 1, None, # C=1, M=None

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

856

857

858 class EcefMatrix(_NamedTuple):

859 '''A rotation matrix.

860 '''

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

862 '_1_0_', '_1_1_', '_1_2_',

863 '_2_0_', '_2_1_', '_2_2_')

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

865

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

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

868 '''

869 _NamedTuple._validate(self, _OK=True)

870

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

872 '''New L{EcefMatrix} matrix.

873

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

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

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

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

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

879

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

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

882 '''

883 t = sa, ca, sb, cb

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

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

886

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

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

889

890 else: # build matrix from the following quaternion operations

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

892 # or

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

894 # where

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

896

897 # Local X axis (east) in geocentric coords

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

899 # Local Y axis (north) in geocentric coords

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

901 # Local Z axis (up) in geocentric coords

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

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

904 cb, -sb * sa, sb * ca,

905 _0_0, ca, sa)

906

907 return _NamedTuple.__new__(cls, *t)

908

909 def column(self, column):

910 '''Get matrix B{C{column}} as 3-tuple.

911 '''

912 if 0 <= column < 3:

913 return self[column::3]

914 raise _IndexError(column=column)

915

916 @Property_RO

917 def _column_0(self):

918 return self.column(0)

919

920 @Property_RO

921 def _column_1(self):

922 return self.column(1)

923

924 @Property_RO

925 def _column_2(self):

926 return self.column(2)

927

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

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

930

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

932 '''

933 return self.classof(*self)

934

935 __copy__ = __deepcopy__ = copy

936

937 def multiply(self, other):

938 '''Matrix multiply M{M0' ⋅ M} this matrix transposed with

939 an other matrix.

940

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

942

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

944

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

946 '''

947 _xinstanceof(EcefMatrix, other=other)

948

949 # like LocalCartesian.MatrixMultiply, transposed(self) X other

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

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

952 return _xnamed(EcefMatrix(*M), EcefMatrix.multiply.__name__)

953

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

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

956

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

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

959

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

961

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

963 '''

964 if xyz0:

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

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

967

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

969

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

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

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

973 return (fdot(xyz, *self._column_0),

974 fdot(xyz, *self._column_1),

975 fdot(xyz, *self._column_2))

976

977 def row(self, row):

978 '''Get matrix B{C{row}} as 3-tuple.

979 '''

980 if 0 <= row < 3:

981 r = row * 3

982 return self[r:r+3]

983 raise _IndexError(row=row)

984

985 @Property_RO

986 def _row_0(self):

987 return self.row(0)

988

989 @Property_RO

990 def _row_1(self):

991 return self.row(1)

992

993 @Property_RO

994 def _row_2(self):

995 return self.row(2)

996

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

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

999

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

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

1002

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

1004

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

1006 '''

1007 if xyz0:

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

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

1010

1011 _xyz = _1_0_1T + xyz

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

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

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

1015 xyz_ = (fdot(_xyz, xyz0[0], *self._row_0),

1016 fdot(_xyz, xyz0[1], *self._row_1),

1017 fdot(_xyz, xyz0[2], *self._row_2))

1018 else:

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

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

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

1022 xyz_ = (fdot(xyz, *self._row_0),

1023 fdot(xyz, *self._row_1),

1024 fdot(xyz, *self._row_2))

1025 return xyz_

1026

1027

1028 class Ecef9Tuple(_NamedTuple):

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

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

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

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

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

1034 '''

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

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

1037

1038 _IndexM = _Names_.index(_M_) # for ._M_x_M

1039

1040 @property_RO

1041 def _CartesianBase(self):

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

1043 '''

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

1045 return C

1046

1047 @deprecated_method

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

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

1050 return self.toDatum(datum2)

1051

1052 @Property_RO

1053 def lam(self):

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

1055 '''

1056 return self.philam.lam

1057

1058 @Property_RO

1059 def lamVermeille(self):

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

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

1062

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

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

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

1066 '''

1067 x, y = self.x, self.y

1068 if y > EPS0:

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

1070 elif y < -EPS0:

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

1072 else: # y == 0

1073 r = PI if x < 0 else _0_0

1074 return Lam(Vermeille=r)

1075

1076 @Property_RO

1077 def latlon(self):

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

1079 '''

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

1081

1082 @Property_RO

1083 def latlonheight(self):

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

1085 '''

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

1087

1088 @Property_RO

1089 def latlonheightdatum(self):

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

1091 '''

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

1093

1094 @Property_RO

1095 def latlonVermeille(self):

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

1097

1098 @see: Property C{lonVermeille}.

1099 '''

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

1101

1102 @Property_RO

1103 def lonVermeille(self):

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

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

1106

1107 @see: Property C{lamVermeille}.

1108 '''

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

1110

1111 def _T_x_M(self, T):

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

1113 '''

1114 t = list(self)

1115 M = self._IndexM

1116 t[M] = T.multiply(t[M])

1117 return self.classof(*t)

1118

1119 @Property_RO

1120 def phi(self):

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

1122 '''

1123 return self.philam.phi

1124

1125 @Property_RO

1126 def philam(self):

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

1128 '''

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

1130

1131 @Property_RO

1132 def philamheight(self):

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

1134 '''

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

1136

1137 @Property_RO

1138 def philamheightdatum(self):

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

1140 '''

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

1142

1143 @Property_RO

1144 def philamVermeille(self):

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

1146

1147 @see: Property C{lamVermeille}.

1148 '''

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

1150

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

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

1153 C{Cartesian}.

1154

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

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

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

1158 or C{None}.

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

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

1161

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

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

1164

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

1166 '''

1167 if Cartesian in (None, Vector4Tuple):

1168 r = self.xyzh

1169 elif Cartesian is Vector3Tuple:

1170 r = self.xyz

1171 else:

1172 _xsubclassof(self._CartesianBase, Cartesian=Cartesian)

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

1174 return r

1175

1176 def toDatum(self, datum2):

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

1178

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

1180

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

1182

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

1184 '''

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

1186 r = self.copy()

1187 else:

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

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

1190 # and returns another Ecef9Tuple iff LatLon is None

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

1192 return r

1193

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

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

1196

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

1198 or C{None}.

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

1200 B{C{LatLon}} keyword arguments.

1201

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

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

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

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

1206

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

1208 '''

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

1210 d = kwds[_datum_]

1211 if LatLon is None:

1212 r = LatLon3Tuple(self.lat, self.lon, kwds[_height_], name=kwds[_name_]) # PYCHOK Ecef9Tuple

1213 if d:

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

1215 else:

1216 if d is None: # remove None datum

1217 _ = kwds.pop[_datum_]

1218 r = LatLon(self.lat, self.lon, **kwds) # PYCHOK Ecef9Tuple

1219 _xdatum(_xattr(r, datum=self.datum), self.datum) # PYCHOK Ecef9Tuple

1220 return r

1221

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

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

1224

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

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

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

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

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

1230

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

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

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

1234

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

1236 '''

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

1238

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

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

1241

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

1243 C{None}.

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

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

1246

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

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

1249

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

1251 '''

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

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

1254

1255 @Property_RO

1256 def xyz(self):

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

1258 '''

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

1260

1261 @Property_RO

1262 def xyzh(self):

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

1264 '''

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

1266

1267

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

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

1270 '''

1271 if Ecef is None:

1272 Ecef = EcefKarney

1273 else:

1274 _xinstanceof(*_Ecefs, Ecef=Ecef)

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

1276

1277

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

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

1280 '''

1281 if issubclassof(Ecef, _EcefBase):

1282 return Ecef

1283 raise _TypesError(_Ecef_, Ecef, *_Ecefs)

1284

1285

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

1287 EcefFarrell21, EcefFarrell22)

1288

1289__all__ += _ALL_DOCS(_EcefBase)

1290

1291# **) MIT License

1292#

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

1294#

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

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

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

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

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

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

1301#

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

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

1304#

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

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

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

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

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

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

1311# OTHER DEALINGS IN THE SOFTWARE.