Coverage for pygeodesy/ecef.py: 95%

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

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

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

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

8into pure Python class L{EcefKarney}, class L{EcefFukushima} from I{Toshio Fukushima}'s U{Fortran

9<https://www.ResearchGate.net/publication/277721539>} version, class L{EcefSudano} based on

10I{John Sudano}'s U{paper<https://www.ResearchGate.net/publication/3709199>}, class L{EcefUPC}

11using the I{Universitat Politècnica de Catalunya}'s U{method, page 186

12<https://GSSC.ESA.int/navipedia/GNSS_Book/ESA_GNSS-Book_TM-23_Vol_I.pdf>}, class L{EcefVeness}

13transcoded from I{Chris Veness}' JavaScript classes U{LatLonEllipsoidal, Cartesian

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

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

16and classes L{EcefFarrell22} and L{EcefFarrell22} from I{Jay A. Farrell}'s U{Table 2.1 and 2.2

17<https://Books.Google.com/books?id=fW4foWASY6wC>}, page 29-30.

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

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

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

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

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

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

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

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

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

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

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

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

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

37(2002) 76, page 451-454.

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

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

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

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

43in U{An analytical method to transform geocentric into geodetic coordinates

44<https://DOI.org/10.1007/s00190-010-0419-x>}, J. Geodesy (2011) 85, page 105-117. See U{Geocentric

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

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

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

49nm (7 nanometers) for the WGS84 ellipsoid. See U{Geocentric coordinates

50<https://GeographicLib.SourceForge.io/C++/doc/geocentric.html>} for further information on the errors.

52@note: The C{reverse} methods of all C{Ecef...} classes return by default C{INT0} as the (geodetic)

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

54C{lon00} to configure that value.

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

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

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

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

60'''

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

63 _xinstanceof, _xsubclassof, typename # _args_kwds_names

64from pygeodesy.constants import EPS, EPS0, EPS02, EPS1, INT0, PI, PI_2, _0_0, _0_5, \

65 _1_0, _1_0_1T, _1_5, _2_0, _3_0, _4_0, _6_0, _90_0, \

66 _copysign_1_0, _isNAN, _isNAN0, _over, isnon0 # PYCHOK used!

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

68from pygeodesy.ecefLocals import _EcefLocal

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

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

71 _xattr, _xdatum, _xkwds, _xkwds_get

72from pygeodesy.fmath import cbrt, _fdotf, hypot, hypot1, hypot2_

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

74# from pygeodesy.internals import typename # from .basics

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

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

77 _x_, _xyz_, _y_, _z_

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

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

80from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \

81 PhiLam2Tuple, Vector3Tuple, Vector4Tuple

82from pygeodesy.props import deprecated_method, deprecated_property, Property_RO, \

83 property_RO, property_ROver

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

85from pygeodesy.units import _isRadius, Degrees, Degrees_, Height, Int, Lam, Lat, \

86 Lon, Meter, Phi, Scalar, Scalar_

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

88 sincos2, sincos2_, sincos2d_

89# from pygeodesy.vector3d import Vector3d # _MODS

91from math import cos, degrees, fabs, radians, sqrt

93__all__ = _ALL_LAZY.ecef

94__version__ = '26.09.06'

96_Ecef_ = 'Ecef'

97_prolate_ = 'prolate'

98_TOL = 1.e-12 # degrees > 1.e-14

99_TRIPS = 33 # 8..9 sufficient

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

101

102

103def _Degrees2Radians(tol): # for EcefUPC

104 return Degrees_(tol=tol, low=EPS, Error=EcefError).toRadians()

105

106

107class _EcefBase(_NamedBase):

108 '''(INTERNAL) Base class for C{Ecef*} convertor classes.

109 '''

110 _datum = _WGS84

111 _e_e2 = None

112 _E = _EWGS84

113 _isYou = False

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

115

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

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

118

119 @arg a_ellipsoid: An ellipsoid (L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple})

120 or a datum (L{Datum}) or the ellipsoid's equatorial

121 radius (C{meter}).

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

123 C{B{a_ellipsoid} is scalar}.

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

125 C{reverse} method.

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

127

128 @raise EcefError: If B{C{a_ellipsoid}} is not an L{Ellipsoid}, L{Ellipsoid2},

129 L{a_f2Tuple} or L{Datum} instance or a positive C{scalar}

130 or if B{C{f}} is not C{scalar} and less than C{1.0}.

131 '''

132 try:

133 E = a_ellipsoid

134 if f is None:

135 pass

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

137 E = a_f2Tuple(E, f)

138 else:

139 raise ValueError() # _invalid_

140

141 if not _isin(E, _EWGS84, _WGS84):

142 d = _ellipsoidal_datum(E, **name)

143 E = d.ellipsoid

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

145 raise ValueError() # _invalid_

146 self._datum = d

147 self._E = E

148

149 if self._isYou:

150 E = self.ellipsoid

151 e2 = E.a2 - E.b2

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

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

154 self._e_e2 = sqrt(e2), e2

155

156 except (TypeError, ValueError) as x:

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

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

159

160 if name:

161 self.name = name

162 if lon00 is not INT0:

163 self.lon00 = lon00

164

165 def __eq__(self, other):

166 '''Compare this and an other Ecef.

167

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

169

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

171 '''

172 return other is self or (isinstance(other, type(self)) and

173 other.ellipsoid == self.ellipsoid)

174

175 @Property_RO

176 def datum(self):

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

178 '''

179 return self._datum

180

181 @Property_RO

182 def ellipsoid(self):

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

184 '''

185 return self._E

186

187 @Property_RO

188 def equatoradius(self):

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

190 '''

191 return self.ellipsoid.a

192

193 a = equatorialRadius = equatoradius # Karney property

194

195 @Property_RO

196 def flattening(self): # Karney property

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

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

199 '''

200 return self.ellipsoid.f

201

202 f = flattening

203

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

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

206

207 @note: From C{Karney}'s, let C{v} be a unit vector located at C{(lat,

208 lon, h)}. We can express C{v} as column vectors in one of two

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

210 components are relative to a local coordinate system at C{C(lat0,

211 lon0, h0)}) or as C{v0} in geocentric C{x, y, z} coordinates.

212 Then, M{v0 = M ⋅ v1} where C{M} is the rotation matrix.

213 '''

214 if _philam: # lat, lon in radians

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

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

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

218 else:

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

220

221 E = self.ellipsoid

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

223 H = _isNAN0(h)

224 c = (H + n) * ca

225 x = cb * c

226 y = sb * c

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

228

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

230 n = self._name__(name)

231 return Ecef9Tuple(x, y, z, lat, lon, h, 0, # C=0, forward

232 m, self.datum, name=n)

233

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

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

236

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

238 latitude (C{degrees}).

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

240 (C{degrees}).

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

242 the surface of the ellipsoid.

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

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

245

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

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

248 C{(lat, lon, height)}, case C{C} (0, forward), rotation matrix

249 C{M} (L{EcefMatrix} or C{None}) and C{datum}.

250

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

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

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

254

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

256 and avoid double angle conversions.

257 '''

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

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

260

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

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

263 in I{radians}.

264

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

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

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

268 the surface of the ellipsoid.

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

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

271

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

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

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

275

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

277 '''

278 try: # like function C{_llhn4} below

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

280 except (TypeError, ValueError) as x:

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

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

283

284 @property_ROver

285 def _Geocentrics(self):

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

287 '''

288 return (Ecef9Tuple, # overwrite property_ROver

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

290

291 @property

292 def lon00(self):

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

294 '''

295 return self._lon00

296

297 @lon00.setter # PYCHOK setter!

298 def lon00(self, lon00):

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

300 '''

301 self._lon00 = Degrees(lon00=lon00)

302

303 def _Matrix(self, *sa_ca_sb_cb):

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

305 '''

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

307

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

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

310 '''

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

312

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

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

315

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

317 coordinate (C{meter}).

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

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

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

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

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

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

324

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

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

327 ones C{(x, y, z)}, case indicator C{C} (C{int} 1..5), rotation matrix

328 C{M} (L{EcefMatrix} or C{None}) and C{datum}.

329

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

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

332 '''

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

334

335 E = self.ellipsoid

336 i = None

337

338 sa, ca, sb, cb, h, p, C = _norm7(y, x, z, E)

339 if C: # PYCHOK no cover

340 pass # too high, too far

341

342 elif p < EPS: # near polar

343 p = 0 # force lon00

344 sa = _copysign_1_0(z)

345 ca = _0_0

346 h = fabs(z) - E.b

347 C = 2 # polar

348

349 elif E.e4: # E.isEllipsoidal

350 s, t = _equatorial2(E, p, z)

351 r = fsumf_(s, t, -E.e4)

352 if t or r > 0:

353 try:

354 sa, ca, h, i, C = self._reverse5(1, h, p, z, y, x, r, s, t)

355 except (TypeError, ValueError) as X:

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

357 raise EcefError(t, cause=X)

358

359 else: # near equatorial plane: e = E.e4 * q == 0 and r <= 0

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

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

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

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

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

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

366 sa, ca, h = _equatorial3(E, s, z)

367 C = 3 # equatorial

368

369 else: # E.isSpherical: E.e4 == 0

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

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

372 sa, ca, _ = _norm3((z if h else _1_0), p)

373 h -= E.a

374 C = 4 # spherical

375

376 lat = atan1d(sa, ca)

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

378 lon = self._polon(sb, cb, p, **lon00_name)

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

380 return Ecef9Tuple(x, y, z, lat, lon, h, C, m, self.datum,

381 iteration=i, name=self._name__(name)) # PYCHOK return

382

383 def _reverse5(self, *C_h_p_z_y_x_r_s_t): # PYCHOK no cover

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

385 self._notOverloaded(*C_h_p_z_y_x_r_s_t)

386

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

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

389

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

391

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

393 '''

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

395

396

397class EcefError(_ValueError):

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

399 '''

400 pass

401

402

403class EcefFarrell21(_EcefBase):

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

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

406 books?id=fW4foWASY6wC>}, page 29, aka the I{Heikkinen application} of U{Ferrari's

407 solution<https://WikiPedia.org/wiki/Geographic_coordinate_conversion>}.

408

409 @see: Classes L{EcefFarrell22} and L{EcefVeness}.

410 '''

411 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature

412 E = self.ellipsoid

413 a = E.a

414 a2 = E.a2

415 b2 = E.b2

416 e2 = E.e2

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

418 e4 = E.e4

419

420 z2 = z**2 # names as page 29

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

422

423 p2 = p**2

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

425 F = b2 * z2 * 54

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

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

428 c = cbrt(s + c + _1_0)

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

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

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

432 Q1 = Q + _1_0

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

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

435 -P * p2 * _0_5)

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

437 r = p + r * e2

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

439

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

441 z += e2_ * v * z # lat = atan1d(z, p)

442 return z, p, h, None, C

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

444

445

446class EcefFarrell22(_EcefBase):

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

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

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

450

451 @see: Classes L{EcefFarrell21} and L{EcefVeness}.

452 '''

453 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature

454 E = self.ellipsoid

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

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

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

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

459 if c:

460 h = p / fabs(c)

461 if s:

462 h -= E.roc1_(s)

463 else:

464 h -= a

465# C = 3 # XXX 1?

466 else:

467 h = fabs(z) - b

468 C = 2

469 # lat = atan1d(s, c)

470 return s, c, h, None, C

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

472

473

474class EcefFukushima(_EcefBase):

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

476 transcoded from I{Toshio Fukushima}'s U{Fortran<https://www.ResearchGate.net/publication/277721539>}

477 implementation.

478

479 @see: Fukushima, T. U{Transformation from Cartesian to Geodetic Coordinates Accelerated by

480 Halley’s Method<https://www.researchgate.net/publication/227215135>} and Eleiche, M.

481 U{A comparison between Fukushima-Halley algorithm and Trilateration algorithm for

482 geodetic conversion<https://link.Springer.com/article/10.1007/s12145-022-00779-7>}.

483 '''

484 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature

485 E = self.ellipsoid

486 a = E.a

487 e2 = E.e2

488 e4 = _1_5 * E.e4 # e4T

489 ec = _1_0 - E.f # sqrt(_1_0 - e2)

490# assert (a * ec) == E.b

491

492 za = fabs(z)

493 s0 = za / a

494 zc = s0 * ec

495 pn = p / a

496 # Newton Correction Factors

497 c0 = pn * ec

498 c2 = c0**2

499 c3 = c2 * c0

500 s2 = s0**2

501 s3 = s2 * s0

502# a2 = s2 + c2

503 a0 = hypot(s0, c0) # sqrt(a2)

504 a3 = a0**3 # a0 * a2

505 d0 = a3 * zc + e2 * s3

506 f0 = a3 * pn - e2 * c3

507 # Halley Correction Factor

508 b0 = e4 * s2 * c2 * pn * (a0 - ec)

509 sa = d0 * f0 - b0 * s0

510 ca = (f0**2 - b0 * c0) * ec

511

512 # lat = atan1d(sa, ca)

513 h = hypot(ec * sa, ca)

514 h = fsumf_(p * ca, za * sa, -h * a)

515 h = _over(h, hypot(sa, ca))

516 return sa, ca, h, None, C

517

518

519class EcefKarney(_EcefBase):

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

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

522 classGeographicLib_1_1Geocentric.html>} methods.

523

524 @note: In general, there are multiple solutions and the result which minimizes C{height} is

525 returned, i.e., the C{(lat, lon)} corresponding to the closest point on the ellipsoid.

526 If there are still multiple solutions with different latitudes (applies only if C{z}

527 = 0), then the solution with C{lat} > 0 is returned. If there are still multiple

528 solutions with different longitudes (applies only if C{x} = C{y} = 0), then C{lon00}

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

530 below M{−E.a * (1 − E.e2) / sqrt(1 − E.e2 * sin(lat)**2)}. Like C{forward} above,

531 M{v1 = Transpose(M) ⋅ v0}.

532 '''

533 def _reverse5(self, C, h, p, z, y, x, r, s, q): # PYCHOK unused y, x

534 E = self.ellipsoid

535 e = E.e4 * q # p renamed to s

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

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

538 d = s = e * s / _4_0 # s = r^3 * s

539 u = r = r / _6_0

540 r2 = r**2

541 r3 = r2 * r

542 t3 = r3 + s

543 d *= t3 + r3

544 if d < 0:

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

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

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

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

549 else:

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

551 # minimizes loss of precision due to cancellation. The

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

553 # in definition of u.

554 if d > 0:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

569 if E.f < 0:

570 k1 -= E.e2

571 else:

572 k2 += E.e2

573 sa, ca, h = _norm3(z / k1, p / k2)

574 h *= k1 - E.e21

575 return sa, ca, h, None, C

576

577

578class EcefSudano(_EcefBase):

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

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

581 '''

582 _TOL = \

583 _tol = EPS

584

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

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

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

588

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

590

591 @see: L{Parent method<_EcefBase.reverse>} for all other information.

592

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

594 C{scalar} for C{scalar} B{C{xyz}} or no convergence for C{B{tol}}.

595 '''

596 if tol != self._TOL:

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

598 return _EcefBase.reverse(self, xyz, y=y, z=z, M=M, **lon00_name)

599

600 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature

601 E = self.ellipsoid

602 e = E.e2 * E.a

603 d = e - p

604

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

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

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

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

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

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

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

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

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

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

615 tol = self._tol

616 _S2 = Fsum(sa).fsum2f_

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

618 ca2 = _1_0 - sa**2

619 if ca2 < EPS02:

620 break

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

622 if fabs(D) < EPS0:

623 break

624 ca = sqrt(ca2)

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

626 if fabs(D) < tol:

627 break

628 else: # PYCHOK no cover

629 raise ValueError(Fmt.no_convergence(fabs(D), tol))

630

631 sa = copysign0(sa, z)

632 # lat = atan1d(sa, ca)

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

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

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

636 return sa, ca, h, i, C

637

638 @deprecated_property

639 def tolerance(self):

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

641 return self._tol

642

643 @tolerance.setter # PYCHOK setter!

644 def tolerance(self, tol):

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

646

647

648class EcefUPC(_EcefBase):

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

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

651 '''

652 _TOL = _TOL

653 _tol = _Degrees2Radians(_TOL)

654

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

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

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

658

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

660

661 @see: L{Parent method<_EcefBase.reverse>} for all other information.

662

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

664 C{scalar} for C{scalar} B{C{xyz}} or no convergence for C{B{tol}}.

665 '''

666 if tol != _TOL:

667 self._tol = _Degrees2Radians(tol)

668 return _EcefBase.reverse(self, xyz, y=y, z=z, M=M, **lon00_name)

669

670 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature

671 E = self.ellipsoid

672 a = E.a

673 e2 = E.e2 # signed

674

675 za = fabs(z)

676 ph_ = atan1(za, E.e21 * p)

677 tol = self._tol

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

679 s, c = sincos2(ph_)

680 N = a / sqrt(_1_0 - s**2 * e2) # N + h == N + p / c - N == p / c

681 ca = p - N * c * e2 # == p * (1 - N * e2 / (N + h)) == p * (1 - N * e2 * c / p)

682 ph = atan1(za, ca) # atan1(z / p, 1 - N * e2 / (N + h)) == atan1(z, ca)

683 r = fabs(ph - ph_)

684 if r < tol:

685 # lat = copysign0(degrees(ph), z)

686 # == atan1d(z, ca)

687 h = p / c - N

688 break

689 ph_ = ph

690 else: # PYCHOK no cover

691 r, tol = map1(degrees, r, tol)

692 raise ValueError(Fmt.no_convergence(r, tol))

693 return z, ca, h, i, C

694

695

696class EcefVeness(_EcefBase):

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

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

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

700

701 @note: Uses B. R. Bowring’s formulation for μm precision in concise form U{The accuracy of

702 geodetic latitude and height equations<https://www.ResearchGate.net/publication/233668213>},

703 Survey Review, Vol 28, 218, Oct 1985.

704

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

706 resources/guide-coordinate-systems-great-britain.pdf>}, section I{B) Converting between 3D

707 Cartesian and ellipsoidal latitude, longitude and height coordinates}.

708

709 @see: Toms, Ralph M. U{An Efficient Algorithm for Geocentric to Geodetic Coordinate Conversion

710 <https://www.OSTI.gov/scitech/biblio/110235>}, Sept 1995 and U{An Improved Algorithm for

711 Geocentric to Geodetic Coordinate Conversion<https://www.OSTI.gov/scitech/servlets/purl/231228>},

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

713 '''

714 def _reverse5(self, C, h, p, z, *unused): # PYCHOK signature

715 # assert h >= p > 0 # h = hypot(z, p)

716 E = self.ellipsoid

717 a = E.a

718 B = E.b * E.e22

719 # parametric latitude (Bowring eqn 17, replaced)

720 t = (E.b * z) / (a * p) * (B / h + _1_0) # theta

721 c = _1_0 / hypot1(t) # t == atan2(z * a, p * E.b)

722 s = c * t # s, c == sincos2(t)

723 # geodetic latitude (Bowring eqn 18)

724 sa, ca, _ = _norm3(z + s**3 * B,

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

726 h = E._heightB(sa, ca, z, p) # height (Bowring eqn 7)

727 # lat = atan1d(sa, ca)

728 return sa, ca, h, None, C

729

730

731class EcefYou(_EcefBase):

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

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

734 for I{non-prolate} ellipsoids.

735

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

737 ellipsoidal coordinates<https://Espace.Curtin.edu.AU/bitstream/handle/20.500.11937/11589/

738 115114_9021_geod2ellip_final.pdf>} Studia Geophysica et Geodaetica, 2008, 52, pages 1-18

739 and U{PyMap3D<https://PyPI.org/project/pymap3d>}.

740 '''

741 _isYou = True

742

743 def _reverse5(self, C, h, p, z, y, x, *unused): # PYCHOK signature

744 E = self.ellipsoid

745 a, b = E.a, E.b

746 e, e2 = self._e_e2

747

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

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

750 u *= _0_5

751 if u > EPS02:

752 u = sqrt(u)

753 q = hypot(u, e)

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

755 sB, cB = sincos2(B)

756 if cB and sB:

757 q *= a

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

759 if isnon0(d):

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

761 sB, cB = sincos2(B)

762 elif u < (-EPS02):

763 raise EcefError(u=u, txt=_singular_)

764 else: # near polar # PYCHOK no cover

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

766

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

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

769 h = neg(h) # inside ellipsoid

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

771 return (a * sB), (b * cB), h, None, C

772

773

774class EcefMatrix(_NamedTuple):

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

776

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

778 Geographic_coordinate_conversion#From_ECEF_to_ENU>} and

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

780 '''

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

782 '_1_0_', '_1_1_', '_1_2_',

783 '_2_0_', '_2_1_', '_2_2_')

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

785

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

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

788 '''

789 _NamedTuple._validate(self, underOK=True)

790

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

792 '''New L{EcefMatrix} matrix.

793

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

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

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

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

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

799

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

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

802 '''

803 t = sa, ca, sb, cb

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

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

806

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

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

809

810 else: # build matrix from the following quaternion operations

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

812 # or

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

814 # where

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

816

817 # Local X axis (East) in geocentric coords

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

819 # Local Y axis (North) in geocentric coords

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

821 # Local Z axis (Up) in geocentric coords

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

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

824 cb, -sb * sa, sb * ca,

825 _0_0, ca, sa)

826

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

828

829 def column(self, column):

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

831 '''

832 if 0 <= column < 3:

833 return self[column::3]

834 raise _IndexError(column=column)

835

836 @property_RO

837 def _columns(self):

838 for c in range(3):

839 yield self[c::3]

840

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

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

843

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

845 '''

846 return self.classof(*self)

847

848 __copy__ = __deepcopy__ = copy

849

850 @Property_RO

851 def matrix3(self):

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

853 '''

854 return tuple(self._rows)

855

856 @Property_RO

857 def matrixTransposed3(self):

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

859 '''

860 return tuple(self._columns)

861

862 def multiply(self, other):

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

864

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

866

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

868

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

870 '''

871 _xinstanceof(EcefMatrix, other=other)

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

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

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

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

876

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

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

879

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

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

882

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

884

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

886 '''

887 if xyz0:

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

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

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

891

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

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

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

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

896

897 def row(self, row):

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

899 '''

900 if 0 <= row < 3:

901 r = row * 3

902 return self[r:r+3]

903 raise _IndexError(row=row)

904

905 @property_RO

906 def _rows(self):

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

908 yield self[r:r+3]

909

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

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

912

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

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

915

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

917

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

919 '''

920 if xyz0:

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

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

923 _xyz = _1_0_1T + xyz

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

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

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

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

928 else:

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

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

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

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

933 return tuple(xyz_)

934

935

936class Ecef9Tuple(_NamedTuple, _EcefLocal):

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

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

939 optionally, rotation matrix C{M} (L{EcefMatrix} or C{None}) and C{datum},

940 with C{lat} and C{lon} in C{degrees} and C{x}, C{y}, C{z} and C{height} in

941 C{meter}, conventionally. Case C{C=0} means C{x, y,z} from foward, C{C=1}

942 C{lat, lon, height} from reverse, C{C=2} near-polar C{lat, lon}, C{C=3}

943 near-equatorial C{lat, lon}, C{C=4} spherical C{lat, lon} and C{C=5} means

944 the C{height} exceeds C{datum}'s C{ellipsoid.heightMax}.

945 '''

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

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

948

949 @property_ROver

950 def _CartesianBase(self):

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

952 '''

953 return _MODS.cartesianBase.CartesianBase # overwrite property_ROver

954

955 @deprecated_method

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

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

958 return self.toDatum(datum2)

959

960 @property_RO

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

962 return self

963

964 @property_RO

965 def ellipsoid(self):

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

967 '''

968 return (self.datum or _WGS84).ellipsoid

969

970 @Property_RO

971 def lam(self):

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

973 '''

974 return self.philam.lam

975

976 @Property_RO

977 def lamVermeille(self):

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

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

980

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

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

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

984 '''

985 x, y = self.x, self.y

986 a = fabs(y)

987 if a > EPS0:

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

989 if y < 0:

990 r = -r

991 else: # y == 0

992 r = PI if x < 0 else _0_0

993 return Lam(Vermeille=r)

994

995 @Property_RO

996 def latlon(self):

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

998 '''

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

1000

1001 @Property_RO

1002 def latlonheight(self):

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

1004 '''

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

1006

1007 @Property_RO

1008 def latlonheightdatum(self):

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

1010 '''

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

1012

1013 @Property_RO

1014 def latlonVermeille(self):

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

1016

1017 @see: Property C{lonVermeille}.

1018 '''

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

1020

1021 @Property_RO

1022 def lonVermeille(self):

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

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

1025

1026 @see: Property C{lamVermeille}.

1027 '''

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

1029

1030 @Property_RO

1031 def Mx(self):

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

1033 '''

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

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

1036

1037 @Property_RO

1038 def phi(self):

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

1040 '''

1041 return self.philam.phi

1042

1043 @Property_RO

1044 def philam(self):

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

1046 '''

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

1048

1049 @Property_RO

1050 def philamheight(self):

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

1052 '''

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

1054

1055 @Property_RO

1056 def philamheightdatum(self):

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

1058 '''

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

1060

1061 @Property_RO

1062 def philamVermeille(self):

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

1064

1065 @see: Property C{lamVermeille}.

1066 '''

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

1068

1069 phiVermeille = phi

1070

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

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

1073 C{Cartesian}.

1074

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

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

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

1078 or C{None}.

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

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

1081

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

1083 is None}.

1084

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

1086 '''

1087 if _isin(Cartesian, None, Vector4Tuple):

1088 r = self.xyzh

1089 elif Cartesian is Vector3Tuple:

1090 r = self.xyz

1091 else:

1092 _xsubclassof(self._CartesianBase, Cartesian=Cartesian)

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

1094 return r

1095

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

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

1098

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

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

1101

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

1103

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

1105 '''

1106 n = _name__(name, _or_nameof=self)

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

1108 r = self.copy(name=n)

1109 else:

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

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

1112 # and returns another Ecef9Tuple iff LatLon is None

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

1114 return r

1115

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

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

1118

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

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

1121 keyword arguments.

1122

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

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

1125 specified or not.

1126

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

1128 '''

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

1130 kwds = _name1__(LatLon_kwds, _or_nameof=self)

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

1132 d = kwds.get(_datum_, LatLon)

1133 if LatLon is None:

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

1135 if d is not None:

1136 # assert d is not LatLon

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

1138 else:

1139 if d is None:

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

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

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

1143 return r

1144

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

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

1147

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

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

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

1151

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

1153 C{B{Vector} is None}.

1154

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

1156

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

1158 '''

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

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

1161

1162# def _T_x_M(self, T):

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

1164# '''

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

1166

1167 @Property_RO

1168 def xyz(self):

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

1170 '''

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

1172

1173 @Property_RO

1174 def xyzh(self):

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

1176 '''

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

1178

1179

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

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

1182 '''

1183 if Ecef is None:

1184 Ecef = EcefKarney

1185 else:

1186 _xinstanceof(*_Ecefs, Ecef=Ecef)

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

1188

1189

1190def _equatorial2(E, p, z):

1191 '''(INTERNAL) Equatorial plane from C{EcefKarney}.

1192 '''

1193 # Treat prolate spheroids by swapping p and z here and by

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

1195 # of method C{EcefKarney._reverse5}

1196 p = (p / E.a)**2

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

1198 return (q, p) if E.f < 0 else (p, q)

1199

1200

1201def _equatorial3(E, s, z):

1202 '''(INTERNAL) Equatorial plane from C{EcefKarney}.

1203 '''

1204 t = E.e4 - s

1205 if E.f < 0:

1206 s, t = t, s

1207 e = E.a

1208 else:

1209 e = E.b2_a

1210 sa, ca, h = _norm3(*map1(sqrt, E._1_e21 * t, s))

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

1212 sa = neg(sa)

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

1214 return sa, ca, h

1215

1216

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

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

1219 '''

1220 try:

1221 lat, lon = latlonh.lat, latlonh.lon

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

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

1224 except AttributeError:

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

1226 try:

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

1228 except (TypeError, ValueError) as x:

1229 t = _lat_, _lon_, _height_

1230 if suffix:

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

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

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

1234

1235

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

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

1238 '''

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

1240 return (y / h, x / h, h) if h > eps else (_0_0, _1_0, h) # copysign_1_0(x)

1241

1242

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

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

1245 '''

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

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

1248 if h > E.heightMax:

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

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

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

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

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

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

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

1256 C = 5

1257 else:

1258 C = 0

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

1260

1261

1262def _xEcef(Ecef): # PYCHOK .latlonBase

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

1264 '''

1265 if issubclassof(Ecef, _EcefBase):

1266 return Ecef

1267 raise _TypesError(_Ecef_, Ecef, *_Ecefs)

1268

1269

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

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

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

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

1274 '''

1275 try:

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

1277 try:

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

1279 if not isinstance(xyz, Types):

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

1281 except AttributeError:

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

1283 except (TypeError, ValueError) as x:

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

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

1286 return t

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

1288

1289

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

1291 if issubclassof(_, _EcefBase) and

1292 _ is not _EcefBase)

1293__all__ += _ALL_DOCS(_EcefBase)

1294

1295# **) MIT License

1296#

1298#

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

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

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

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

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

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

1305#

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

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

1308#

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

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

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

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

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

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

1315# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pygeodesy / ecef.py: 95%

527 statements