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@note: The C{reverse} methods of all C{Ecef...} classes return by default C{INT0} as the (geodetic)

51 longitude for I{polar} ECEF location C{x == y == 0}. Use keyword argument C{lon00} or property

52 C{lon00} to configure that value.

53

54@see: Module L{ltp} and class L{LocalCartesian}, a transcription of I{Charles Karney}'s C++ class

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

56 for conversion between geodetic and I{local cartesian} coordinates in a I{local tangent

57 plane} as opposed to I{geocentric} (ECEF) ones.

58'''

59

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

61 _xinstanceof, _xsubclassof # _args_kwds_names

62 from pygeodesy.constants import EPS, EPS0, EPS02, EPS1, EPS2, EPS_2, INT0, PI, PI_2, \

63 _0_0, _0_0001, _0_01, _0_5, _1_0, _1_0_1T, _N_1_0, \

64 _2_0, _N_2_0, _3_0, _4_0, _6_0, _60_0, _90_0, _N_90_0, \

65 _100_0, _copysign_1_0, isnon0 # PYCHOK used!

66 from pygeodesy.datums import a_f2Tuple, _ellipsoidal_datum, _WGS84, _EWGS84

67# from pygeodesy.ellipsoids import a_f2Tuple, _EWGS84 # from .datums

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

69 _xattr, _xdatum, _xkwds, _xkwds_get

70 from pygeodesy.fmath import cbrt, fdot, Fpowers, hypot, hypot1, hypot2_, sqrt0

71 from pygeodesy.fsums import Fsum, fsumf_, Fmt, unstr

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

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

74 _x_, _xyz_, _y_, _z_

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

76 from pygeodesy.named import _name__, _name1__, _NamedBase, _NamedTuple, _Pass, _xnamed

77 from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \

78 PhiLam2Tuple, Vector3Tuple, Vector4Tuple

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

80# from pygeodesy.streprs import Fmt, unstr # from .fsums

81 from pygeodesy.units import _isRadius, Degrees, Height, Int, Lam, Lat, Lon, Meter, \

82 Phi, Scalar, Scalar_

83 from pygeodesy.utily import atan1, atan1d, atan2d, degrees90, degrees180, \

84 sincos2, sincos2_, sincos2d, sincos2d_

85

86 from math import atan2, cos, degrees, fabs, radians, sqrt

87

88__all__ = _ALL_LAZY.ecef

89__version__ = '24.06.11'

90

91_Ecef_ = 'Ecef'

92_prolate_ = 'prolate'

93_TRIPS = 33 # 8..9 sufficient, EcefSudano.reverse

94_xyz_y_z = _xyz_, _y_, _z_ # _args_kwds_names(_xyzn4)[:3]

95

96

97 class EcefError(_ValueError):

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

99 '''

100 pass

101

102

103 class _EcefBase(_NamedBase):

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

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

106 '''

107 _datum = _WGS84

108 _E = _EWGS84

109 _lon00 = INT0 # arbitrary, "polar" lon for LocalCartesian, Ltp

110

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

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

113

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

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

116 equatorial radius (C{meter}).

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

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

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

120 @kwarg lon00: An arbitrary, I{"polar"} longitude (C{degrees}), see the

121 C{reverse} method.

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

123

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

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

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

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

128 '''

129 try:

130 E = a_ellipsoid

131 if f is None:

132 pass

133 elif _isRadius(E) and isscalar(f):

134 E = a_f2Tuple(E, f)

135 else:

136 raise ValueError() # _invalid_

137

138 if E not in (_EWGS84, _WGS84):

139 d = _ellipsoidal_datum(E, **name)

140 E = d.ellipsoid

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

142 raise ValueError() # _invalid_

143 self._datum = d

144 self._E = E

145

146 except (TypeError, ValueError) as x:

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

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

149

150 if name:

151 self.name = name

152 if lon00 is not INT0:

153 self.lon00 = lon00

154

155 def __eq__(self, other):

156 '''Compare this and an other Ecef.

157

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

159

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

161 '''

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

163 other.ellipsoid == self.ellipsoid)

164

165 @Property_RO

166 def datum(self):

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

168 '''

169 return self._datum

170

171 @Property_RO

172 def ellipsoid(self):

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

174 '''

175 return self._E

176

177 @Property_RO

178 def equatoradius(self):

179 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}).

180 '''

181 return self.ellipsoid.a

182

183 a = equatorialRadius = equatoradius # Karney property

184

185 @Property_RO

186 def flattening(self): # Karney property

187 '''Get the C{ellipsoid}'s flattening (C{scalar}), positive for

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

189 '''

190 return self.ellipsoid.f

191

192 f = flattening

193

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

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

196 '''

197 if _philam: # lat, lon in radians

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

199 lat = Lat(degrees90( lat), Error=EcefError)

200 lon = Lon(degrees180(lon), Error=EcefError)

201 else:

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

203

204 E = self.ellipsoid

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

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

207 x = (h + n) * ca

208

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

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

211 0, m, self.datum,

212 name=self._name__(name))

213

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

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

216

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

218 latitude (C{degrees}).

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

220 (C{degrees}).

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

222 the surface of the ellipsoid.

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

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

225

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

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

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

229 and C{datum} if available.

230

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

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

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

234

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

236 and avoid double angle conversions.

237 '''

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

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

240

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

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

243 in I{radians}.

244

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

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

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

248 the surface of the ellipsoid.

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

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

251

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

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

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

255

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

257 '''

258 try: # like function C{_llhn4} below

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

260 except (TypeError, ValueError) as x:

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

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

263

264 @property_RO

265 def _Geocentrics(self):

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

267 '''

268 _EcefBase._Geocentrics = t = (Ecef9Tuple, # overwrite property_RO

269 _MODS.cartesianBase.CartesianBase)

270 return t

271

272 @Property_RO

273 def _isYou(self):

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

275 '''

276 return isinstance(self, EcefYou)

277

278 @property

279 def lon00(self):

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

281 '''

282 return self._lon00

283

284 @lon00.setter # PYCHOK setter!

285 def lon00(self, lon00):

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

287 '''

288 self._lon00 = Degrees(lon00=lon00)

289

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

291 '''Creation a rotation matrix.

292

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

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

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

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

297

298 @return: An L{EcefMatrix}.

299 '''

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

301

302 def _polon(self, y, x, R, **lon00_name):

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

304 '''

305 return atan2d(y, x) if R else _xkwds_get(lon00_name, lon00=self.lon00)

306

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

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

309 self._notOverloaded(xyz, y=y, z=z, M=M, **lon00_name)

310

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

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

313

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

315

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

317 '''

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

319

320

321 class EcefFarrell21(_EcefBase):

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

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

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

325 '''

326

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

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

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

330 page 29.

331

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

333 coordinate (C{meter}).

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

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

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

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

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

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

340

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

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

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

344 if available.

345

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

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

348 zero division error.

349

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

351 '''

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

353

354 E = self.ellipsoid

355 a = E.a

356 a2 = E.a2

357 b2 = E.b2

358 e2 = E.e2

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

360 e4 = E.e4

361

362 try: # names as page 29

363 z2 = z**2

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

365

366 p = hypot(x, y)

367 p2 = p**2

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

369 F = b2 * z2 * 54

370 t = e4 * p2 * F / G**3

371 t = cbrt(sqrt(t * (t + _2_0)) + t + _1_0)

372 G *= fsumf_(t , _1_0, _1_0 / t)

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

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

375 Q1 = Q + _1_0

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

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

378 -P * p2 * _0_5))

379 r = p + e2 * r0

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

381

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

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

384 lon = self._polon(y, x, p, **lon00_name)

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

386

387 except (ValueError, ZeroDivisionError) as e:

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

389

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

391 1, None, self.datum,

392 name=self._name__(name))

393

394

395 class EcefFarrell22(_EcefBase):

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

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

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

399 '''

400

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

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

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

404 page 30.

405

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

407 coordinate (C{meter}).

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

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

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

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

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

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

414

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

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

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

418 if available.

419

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

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

422 zero division error.

423

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

425 '''

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

427

428 E = self.ellipsoid

429 a = E.a

430 b = E.b

431

432 try: # see EcefVeness.reverse

433 p = hypot(x, y)

434 lon = self._polon(y, x, p, **lon00_name)

435

436 s, c = sincos2(atan2(z * a, p * b)) # == _norm3

437 lat = atan1d(z + s**3 * b * E.e22,

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

439

440 s, c = sincos2d(lat)

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

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

443 else: # polar

444 h = fabs(z) - b

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

446

447 except (ValueError, ZeroDivisionError) as e:

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

449

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

451 1, None, self.datum,

452 name=self._name__(name))

453

454

455 class EcefKarney(_EcefBase):

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

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

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

459

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

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

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

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

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

465 the rotation matrix.

466 '''

467

468 @Property_RO

469 def hmax(self):

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

471 '''

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

473

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

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

476

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

478 coordinate (C{meter}).

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

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

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

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

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

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

485

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

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

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

489 C{datum} if available.

490

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

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

493

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

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

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

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

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

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

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

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

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

503 '''

504 def _norm3(y, x):

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

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

507

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

509

510 E = self.ellipsoid

511 f = E.f

512

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

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

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

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

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

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

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

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

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

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

523 C = 1

524

525 elif E.e4: # E.isEllipsoidal

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

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

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

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

530 if f < 0:

531 p, q = q, p

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

533 e = E.e4 * q

534 if e or r > 0:

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

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

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

538 u = r = r / _6_0

539 r2 = r**2

540 r3 = r2 * r

541 t3 = r3 + s

542 d *= t3 + r3

543 if d < 0:

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

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

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

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

548 else:

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

550 # minimizes loss of precision due to cancellation. The

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

552 # in definition of u.

553 if d > 0:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

568 if f < 0:

569 k1 -= E.e2

570 else:

571 k2 += E.e2

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

573 h *= k1 - E.e21

574 C = 2

575

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

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

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

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

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

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

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

583 q = E.e4 - p

584 if f < 0:

585 p, q = q, p

586 e = E.a

587 else:

588 e = E.b2_a

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

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

591 sa = neg(sa)

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

593 C = 3

594

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

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

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

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

599 h -= E.a

600 C = 4

601

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

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

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

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

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

607

608

609 class EcefSudano(_EcefBase):

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

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

612 '''

613 _tol = EPS2

614

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

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

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

618

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

620 coordinate (C{meter}).

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

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

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

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

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

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

627

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

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

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

631

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

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

634 '''

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

636

637 E = self.ellipsoid

638 e = E.e2 * E.a

639 R = hypot(x, y) # Rh

640 d = e - R

641

642 _a = fabs

643 lat = atan1d(z, R * E.e21)

644 sa, ca = sincos2d(_a(lat))

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

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

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

648 # = ca**2 * (z * ca + E.e2 * E.a * sa - R * sa)

649 # = ca**2 * (z * ca + (E.e2 * E.a - R) * sa)

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

651 # = ca**2 * (E.e2 * E.a / E.e2s2(sa) - R / ca**2)

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

653 # (E.e2 * E.a / E.e2s2(sa) - R / ca**2)

654 _E = EPS_2

655 tol = self.tolerance

656 _S2 = Fsum(sa).fsum2f_

657 _rt = sqrt

658 for i in range(1, _TRIPS):

659 ca2 = _1_0 - sa**2

660 if ca2 < _E: # PYCHOK no cover

661 ca = _0_0

662 break

663 ca = _rt(ca2)

664 r = e / E.e2s2(sa) - R / ca2

665 if _a(r) < _E:

666 break

667 lat = None

668 sa, r = _S2(-z * ca / r, -d * sa / r)

669 if _a(r) < tol:

670 break

671 else:

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

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

674

675 if lat is None:

676 lat = copysign0(atan1d(_a(sa), ca), z)

677 lon = self._polon(y, x, R, **lon00_name)

678

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

680 # because Sudano's Eq (7) doesn't produce the correct height

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

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

683 i, None, self.datum, # C=i, M=None

684 iteration=i, name=self._name__(name))

685

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

687 def tolerance(self):

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

689 '''

690 return self._tol

691

692 @tolerance.setter # PYCHOK setter!

693 def tolerance(self, tol):

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

695

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

697 '''

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

699

700

701 class EcefVeness(_EcefBase):

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

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

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

705

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

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

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

709 '''

710

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

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

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

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

715

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

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

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

719

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

721 coordinate (C{meter}).

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

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

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

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

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

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

728

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

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

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

732

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

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

735

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

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

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

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

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

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

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

743 publication/3709199>}.

744 '''

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

746

747 E = self.ellipsoid

748 a = E.a

749

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) / (a * p) * (_1_0 + b / r)

756 c = _1_0 / hypot1(t)

757 s = c * t

758 # geodetic latitude (Bowring eqn 18)

759 lat = atan1d(z + s**3 * b,

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

761

762 # height above ellipsoid (Bowring eqn 7)

763 sa, ca = sincos2d(lat)

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

765# h = p * ca + z * sa - (a * a / r)

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

767 C = 1

768

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

770 elif p > EPS: # lat arbitrarily zero, equatorial lon

771 C, lat, h = 2, _0_0, (p - a)

772

773 else: # polar lat, lon arbitrarily lon00

774 C, lat, h = 3, (_N_90_0 if z < 0 else _90_0), (fabs(z) - E.b)

775

776 lon = self._polon(y, x, p, **lon00_name)

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

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

779 name=self._name__(name))

780

781

782 class EcefYou(_EcefBase):

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

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

785 for I{non-prolate} ellipsoids.

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=_EWGS84, f=None, **lon00_name): # PYCHOK signature

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

795 self._ee2 = EcefYou._ee2(self.ellipsoid)

796

797 @staticmethod

798 def _ee2(E):

799 e2 = E.a2 - E.b2

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

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

802 return sqrt0(e2), e2

803

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

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

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

807

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

809 coordinate (C{meter}).

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

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

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

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

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

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

816

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

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

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

820 available.

821

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

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

824 ellipsoid is I{prolate}.

825 '''

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

827

828 E = self.ellipsoid

829 e, e2 = self._ee2

830

831 q = hypot(x, y) # R

832 u = Fpowers(2, x, y, z) - e2

833 u = u.fadd_(hypot(u, e * z * _2_0)).fover(_2_0)

834 if u > EPS02:

835 u = sqrt(u)

836 p = hypot(u, e)

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

838 sB, cB = sincos2(B)

839 if cB and sB:

840 p *= E.a

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

842 if isnon0(d):

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

844 sB, cB = sincos2(B)

845 elif u < (-EPS2):

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

847 else:

848 sB, cB = _copysign_1_0(z), _0_0

849

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

851 lon = self._polon(y, x, q, **lon00_name)

852

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

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

855 h = neg(h) # inside ellipsoid

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

857 1, None, self.datum, # C=1, M=None

858 name=self._name__(name))

859

860

861 class EcefMatrix(_NamedTuple):

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

863

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

865 Geographic_coordinate_conversion#From_ECEF_to_ENU>} and

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

867 '''

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

869 '_1_0_', '_1_1_', '_1_2_',

870 '_2_0_', '_2_1_', '_2_2_')

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

872

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

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

875 '''

876 _NamedTuple._validate(self, underOK=True)

877

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

879 '''New L{EcefMatrix} matrix.

880

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

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

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

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

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

886

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

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

889 '''

890 t = sa, ca, sb, cb

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

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

893

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

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

896

897 else: # build matrix from the following quaternion operations

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

899 # or

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

901 # where

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

903

904 # Local X axis (East) in geocentric coords

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

906 # Local Y axis (North) in geocentric coords

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

908 # Local Z axis (Up) in geocentric coords

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

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

911 cb, -sb * sa, sb * ca,

912 _0_0, ca, sa)

913

914 return _NamedTuple.__new__(cls, *t)

915

916 def column(self, column):

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

918 '''

919 if 0 <= column < 3:

920 return self[column::3]

921 raise _IndexError(column=column)

922

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

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

925

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

927 '''

928 return self.classof(*self)

929

930 __copy__ = __deepcopy__ = copy

931

932 @Property_RO

933 def matrix3(self):

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

935 '''

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

937

938 @Property_RO

939 def matrixTransposed3(self):

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

941 '''

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

943

944 def multiply(self, other):

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

946 with an other matrix.

947

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

949

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

951

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

953 '''

954 _xinstanceof(EcefMatrix, other=other)

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

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

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

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

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

960

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

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

963

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

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

966

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

968

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

970 '''

971 if xyz0:

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

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

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

975

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

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

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

979 return (fdot(xyz, *self[0::3]), # .column(0)

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

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

982

983 def row(self, row):

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

985 '''

986 if 0 <= row < 3:

987 r = row * 3

988 return self[r:r+3]

989 raise _IndexError(row=row)

990

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

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

993

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

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

996

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

998

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

1000 '''

1001 if xyz0:

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

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

1004 _xyz = _1_0_1T + xyz

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

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

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

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

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

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

1011 else:

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

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

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

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

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

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

1018 return xyz_

1019

1020

1021 class Ecef9Tuple(_NamedTuple):

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

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

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

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

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

1027 '''

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

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

1030

1031 @property_RO

1032 def _CartesianBase(self):

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

1034 '''

1035 Ecef9Tuple._CartesianBase = C = _MODS.cartesianBase.CartesianBase # overwrite property_RO

1036 return C

1037

1038 @deprecated_method

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

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

1041 return self.toDatum(datum2)

1042

1043 @Property_RO

1044 def lam(self):

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

1046 '''

1047 return self.philam.lam

1048

1049 @Property_RO

1050 def lamVermeille(self):

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

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

1053

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

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

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

1057 '''

1058 x, y = self.x, self.y

1059 if y > EPS0:

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

1061 elif y < -EPS0:

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

1063 else: # y == 0

1064 r = PI if x < 0 else _0_0

1065 return Lam(Vermeille=r)

1066

1067 @Property_RO

1068 def latlon(self):

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

1070 '''

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

1072

1073 @Property_RO

1074 def latlonheight(self):

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

1076 '''

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

1078

1079 @Property_RO

1080 def latlonheightdatum(self):

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

1082 '''

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

1084

1085 @Property_RO

1086 def latlonVermeille(self):

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

1088

1089 @see: Property C{lonVermeille}.

1090 '''

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

1092

1093 @Property_RO

1094 def lonVermeille(self):

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

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

1097

1098 @see: Property C{lamVermeille}.

1099 '''

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

1101

1102 @Property_RO

1103 def phi(self):

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

1105 '''

1106 return self.philam.phi

1107

1108 @Property_RO

1109 def philam(self):

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

1111 '''

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

1113

1114 @Property_RO

1115 def philamheight(self):

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

1117 '''

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

1119

1120 @Property_RO

1121 def philamheightdatum(self):

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

1123 '''

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

1125

1126 @Property_RO

1127 def philamVermeille(self):

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

1129

1130 @see: Property C{lamVermeille}.

1131 '''

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

1133

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

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

1136 C{Cartesian}.

1137

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

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

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

1141 or C{None}.

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

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

1144

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

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

1147

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

1149 '''

1150 if Cartesian in (None, Vector4Tuple):

1151 r = self.xyzh

1152 elif Cartesian is Vector3Tuple:

1153 r = self.xyz

1154 else:

1155 _xsubclassof(self._CartesianBase, Cartesian=Cartesian)

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

1157 return r

1158

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

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

1161

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

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

1164

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

1166

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

1168 '''

1169 n = _name__(name, _or_nameof=self)

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

1171 r = self.copy(name=n)

1172 else:

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

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

1175 # and returns another Ecef9Tuple iff LatLon is None

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

1177 return r

1178

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

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

1181

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

1183 or C{None}.

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

1185 B{C{LatLon}} keyword arguments.

1186

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

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

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

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

1191

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

1193 '''

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

1195 kwds = _name1__(LatLon_kwds, _or_nameof=self)

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

1197 d = kwds.get(_datum_, LatLon)

1198 if LatLon is None:

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

1200 if d is not None:

1201 # assert d is not LatLon

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

1203 else:

1204 if d is None:

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

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

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

1208 return r

1209

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

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

1212

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

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

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

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

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

1218

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

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

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

1222

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

1224 '''

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

1226

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

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

1229

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

1231 C{None}.

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

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

1234

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

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

1237

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

1239 '''

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

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

1242

1243# def _T_x_M(self, T):

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

1245# '''

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

1247

1248 @Property_RO

1249 def xyz(self):

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

1251 '''

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

1253

1254 @Property_RO

1255 def xyzh(self):

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

1257 '''

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

1259

1260

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

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

1263 '''

1264 if Ecef is None:

1265 Ecef = EcefKarney

1266 else:

1267 _xinstanceof(*_Ecefs, Ecef=Ecef)

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

1269

1270

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

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

1273 '''

1274 try:

1275 lat, lon = latlonh.lat, latlonh.lon

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

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

1278 except AttributeError:

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

1280 try:

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

1282 except (TypeError, ValueError) as x:

1283 t = _lat_, _lon_, _height_

1284 if suffix:

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

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

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

1288

1289

1290 def _xEcef(Ecef): # PYCHOK .latlonBase

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

1292 '''

1293 if issubclassof(Ecef, _EcefBase):

1294 return Ecef

1295 raise _TypesError(_Ecef_, Ecef, *_Ecefs)

1296

1297

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

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

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

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

1302 '''

1303 try:

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

1305 try:

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

1307 if not isinstance(xyz, Types):

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

1309 except AttributeError:

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

1311 except (TypeError, ValueError) as x:

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

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

1314 return t

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

1316

1317

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

1319 EcefFarrell21, EcefFarrell22)

1320__all__ += _ALL_DOCS(_EcefBase)

1321

1322# **) MIT License

1323#

1324# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1325#

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

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

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

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

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

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

1332#

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

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

1335#

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

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

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

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

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

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

1342# OTHER DEALINGS IN THE SOFTWARE.