Coverage for pyaxqg/axqgs.py: 97%

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

4u'''Classes L{Ax2QG} and L{Ax3QG} implement a WGS84 biaxial respectively EGM2008 triaxial

5reference ellipsoid and bilinear interpolation of quasi-geoid height C{N} and cartesian

6C{X}, C{Y}, C{Z} from self-contained, U{1-degree, whole Earth grids

7<https://link.Springer.com/article/10.1007/s00190-023-01717-1#Sec21>}.

9Class L{Ax5QG} --a sub-class of L{Ax3QG}-- is based on 5 bivariate spline interpolations

10of rectangular grids from C{SciPy} and requires both C{scipy} and C{numpy} to be installed.

12Each class provides a C{forward} method to transform geodetic lat-, longitude and ellipsoidal

13height to cartesian X, Y, Z and orthometric height and a C{reverse} method for converting

14cartesian to geodetic coordinates and orthometric to ellipsoidal height.

15'''

16# make sure int/int division yields float quotient in Py2-

17from __future__ import division as _; del _ # noqa: E702 ;

19from pyaxqg.__pygeodesy import (AxQGError, AxQG8Tuple, AxyzNgeoid4Tuple,

20 _0_0, _1_0, _90_0, _180_0, _isNAN,

21 _ALL_DOCS, _ALL_OTHER, _FOR_DOCS,

22 _xkwds, _xkwds_get,

23 _isinside, _name_, _NamedBase, Vector3Tuple)

24from pygeodesy import (Ang, NAN, INT0, typename, fdot_, # angles, "consterns", fmath

25 Bounds4Tuple, LatLonNgeoid3Tuple, # namedTuples

26 Property_RO, property_RO, property_ROver, # props

27 Triaxial3, # LLK as _LLK, # triaxials

28 Degrees, Height, Lat, Lon, Meter) # units

30from array import array as _array

31from math import ceil, floor

33__all__ = ()

34__version__ = '26.09.01'

36_1_Degree = Degrees(_1_0)

37_forward_ = 'forward'

38_nS2N = 181

39_nW2E = 361 # PYCHOK used!

40_outside__ = 'outside '

41_region4ax = Bounds4Tuple(-_90_0, -_180_0,

42 _90_0, _180_0, name='AxQG region ')

43_reverse_ = 'reverse'

46class _AxG(tuple):

47 '''(INTERNAL) [_nW2E] * [_nS2N] grid of float or double C{array}s.

48 '''

49 @property_RO

50 def dtype(self):

51 '''Return NumPy C{'f8'} for C{double-}, C{'f4'} for C{single-}precision floating point, otherwise C{None}.

52 '''

53 # dtype=float == 'f8' == numpy.float64, dtype='f4' == numpy.float32, dtype='f2' == numpy.float16

54 return {'d': 'f8', 'f': 'f4'}.get(self[0].typecode, None)

56# @property_RO

57# def iterate(self):

58# # iterate lon [_nW2E] lat [_nS2N]

59# return (k for m in self for k in m)

61 def _Nterpolate(self, c_latI, f_latI, latN_f,

62 c_lonI, f_lonI, lonN_f):

63 # bilinear interpolation at C{(lat, lon), normalized}

64 # in this C{Ax[_nW2E][_nS2N]} grid, col-major ordered

65 Ax = self # GeoidKarney.ev2k and geoids._Dotf

66 if c_latI != f_latI or c_lonI != f_lonI:

67 Me, Mw = Ax[c_lonI], Ax[f_lonI]

68 ne, nw = Me[c_latI], Mw[c_latI]

69 se, sw = Me[f_latI], Mw[f_latI]

70 lonN_f1 = _1_0 - lonN_f # == 1 - (lonN - f_lonN)

71 x = (ne * lonN_f + nw * lonN_f1) * latN_f + \

72 (se * lonN_f + sw * lonN_f1) * (_1_0 - latN_f)

73 else:

74 x = Ax[c_lonI][c_latI]

75 return x

77 @property_RO

78 def transpose(self):

79 # transpose from [_nW2E][_nS2N] to [_nS2N][_nW2E]

80 return (k for r in zip(*self) for k in r)

83class _AxQGbase(_NamedBase):

84 '''(INTERNAL) C{Ax*QG} base class.

85 '''

86 _Ax_grid = None # overloaded with _Ax2 or _Ax3

87 _ellipsoid = None # overloaded with WGS84 bi- or EGM2008 triaxial

88 _kind = 1 # bilinear

89 _latD = \

90 _lonD = _1_Degree

91 _onEPS = 2.515e-11 # min -2.507081e-11 max 2.514483e-11, EPS4 = 4.440892098500626e-12

92 _raiser = False

93 _smooth = None # n/a

94 _triaxial = None # overloaded with TriAxial

96 def __init__(self, raiser=False, **name):

97 '''New C{Ax*QG} transformer instance.

99 @kwarg raiser: If C{True} raise an L{AxQGError} for

100 points outside L{region4} (C{bool}).

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

102 '''

103 if raiser:

104 self.raiser = True

105 if name:

106 self.name = name # or typename(self)

107

108 def _Ax_G(self, ax):

109 # get the C{ax} tuple as an C{_AxG} instance

110 nlon, nlat = self._Ax_G2

111 assert isinstance(ax, tuple), type(ax)

112 assert nlon == len(ax) == _nW2E

113 assert nlat == len(ax[0]) == _nS2N # in _darray, _farray

114 return _AxG(ax)

115

116 @Property_RO

117 def _Ax_G2(self):

118 # get the required C{_Ax_G} shape

119 S, W, N, E = _region4ax

120 return (int(_degN(E, W, self._lonD) + _1_0),

121 int(_degN(N, S, self._latD) + _1_0))

122

123 def axN(self, lat, lon):

124 '''Interpolate the quasi-geoid height C{N} for a geodetic point.

125

126 @arg lat: Latitude (C{degrees}, geodetic).

127 @arg lon: Longitude (C{degrees}, geodetic).

128

129 @return: Normal or quasi-geoid height C{N} (C{meter}) or

130 C{NAN} if C{lat} or C{lon} is outside L{region4}.

131 '''

132 lat, lon, _NAN, _, _ = self._LatLon5(lat, lon, False)

133 return NAN if _NAN else self._axN(lat, lon)

134

135 def _axN(self, lat, lon, raiser=False):

136 # interpolate C{N} at C{(lat, lon)} or C{NAN} if

137 # outside or ... if _isNAN(lat) or _isNAN(lon)

138 if _isinside(lat, lon, 0, _region4ax):

139 c_f_N_f6_ll = self._c_f_N_f6_ll(lat, lon)

140 N = self._Ax_grid(*c_f_N_f6_ll)

141 N = Height(N=N)

142 elif raiser or (raiser is None and self._raiser):

143 raise self._outsidError(lat, lon)

144 else:

145 N = NAN

146 return N

147

148 def axN3(self, x, y, z):

149 '''Interpolate the quasi-geoid height C{N} for a cartesian point.

150

151 @arg x: X coordinate (C{meter}, cartesian).

152 @arg y: Y coordinate (C{meter}, cartesian).

153 @arg z: Z coordinate (C{meter}, cartesian).

154

155 @return: L{LatLonNgeoid3Tuple}C{(lat, lon, N)} with the

156 quasi-geoid height C{N} in C{meter} or C{NAN} if

157 the point is not on the L{triaxial}'s surface.

158 '''

159 return self.reverse(x, y, z, H=0, raiser=False).latlonNgeoid

160

161 def _c_f_N_f6_ll(self, lat, lon):

162 # return (int(ceil), int(floor), Normalized less floor) of C{lat}) + \

163 # (int(ceil), int(floor), Normalized less floor) of C{lon})

164 S, W, _, _ = _region4ax

165 return _c_f_N_f3(lat, S, self._latD) + \

166 _c_f_N_f3(lon, W, self._lonD)

167

168 @property_RO

169 def ellipsoid(self):

170 '''Get the C{WGS84} biaxial or C{EGM2008} triaxial ellipsoid.

171 '''

172 return self._ellipsoid

173

174 def forward(self, lat, lon, height=0, **raiser_name):

175 '''Convert a geodetic C{B{lat}}, C{B{lon}} point and ellipsoidal

176 B{C{height}} to cartesian C{x}, C{y}, C{z} on this triaxial and

177 orthometric height C{H}.

178

179 @arg lat: Latitude (C{degrees}, geodetic).

180 @arg lon: Longitude (C{degrees}, geodetic).

181 @kwarg height: The (ellipsoidal) height (C{meter}, conventionally)

182 or C{None} to ignore height interpolation.

183 @kwarg raiser_name: Use C{B{raiser}=True} to raise an L{AxQGError}

184 if B{C{lat}} or B{C{lon}} is outside L{region4},

185 overriding property C{raiser} (C{bool}) and optional

186 C{B{name}='forward'} (C{str}).

187

188 @return: An L{AxQG8Tuple}C{(x, y, z, H, lat, lon, height, axQG)} with

189 cartesian C{x}, C{y}, C{z} and (orthometric) height C{H} all

190 in C{meter} or C{NAN} and C{axQG} this C{Ax*BG} instance.

191

192 @raise AxQGError: If the geodetic point is outside L{region4} and property

193 C{raiser is True} or keyword argument C{B{raiser}=True}.

194

195 @note: Orthometric height C{(H = h - N)} equals ellipsoidal height C{h}

196 less (quasi-)geoid height C{N}.

197 '''

198 lat, lon, _NAN, raiser, name = self._LatLon5(lat, lon, **raiser_name)

199 if _NAN:

200 x = y = z = H = NAN

201 else:

202 x, y, z, H = self._forward4(raiser, lat, lon, height)

203 return AxQG8Tuple(x, y, z, H, lat, lon, height, self, name=name)

204

205 def _forward3(self, lat, lon): # must be overloaded!

206 self._notOverloaded(lat, lon) # PYCHOK no cover

207

208 def _forward4(self, raiser, lat, lon, height): # in .__main__

209 # C{forward} core, returning C{(easting, northing, H)}

210 H = NAN if height is None or _isNAN(height) else (

211 Height(height) - self._axN(lat, lon, raiser))

212 x, y, z = self._forward3(lat, lon)

213 return x, y, z, H

214

215 def forwardOn(self, lat, lon, **raiser_name):

216 '''Interpolate the (reference) triaxial's geoid grid C{x}, C{y} and C{z}

217 and geoid height C{N} at a geodetic C{lat}- and C{lon}gitude.

218

219 @arg lat: Latitude (C{degrees}, geodetic).

220 @arg lon: Longitude (C{degrees}, geodetic).

221 @kwarg raiser_name: Use C{B{raiser}=True} to raise an L{AxQGError}

222 if B{C{lat}} or B{C{lon}} is outside L{region4},

223 overriding property C{raiser} (C{bool}) and optional

224 C{B{name}='forwardOn'} (C{str}).

225

226 @return: An L{AxyzNgeoid4Tuple}C{(x, y, z, N)} with cartesian C{x}, C{y},

227 C{z} and geoid height C{N}, all I{interpolated} and in C{meter}

228 or C{NAN}.

229

230 @raise AxQGError: If the geodetic point is outside L{region4} and property

231 C{raiser is True} or keyword argument C{B{raiser}=True}.

232

233 @see: B{Geoid grid file format} in U{Supplementary File 3

234 <https://link.Springer.com/article/10.1007/s00190-023-01717-1#Sec21>}.

235 '''

236 kwds = _xkwds(raiser_name, name=typename(self.forwardOn))

237 lat, lon, _NAN, raiser, name = self._LatLon5(lat, lon, **kwds)

238 if _NAN:

239 x = y = z = N = NAN

240 elif _isinside(lat, lon, 0, _region4ax):

241 c_f_N_f6_ll = self._c_f_N_f6_ll(lat, lon)

242 N = self._Ax_grid(*c_f_N_f6_ll) # imports ax_grids

243 x, y, z = self._forwardOn3(c_f_N_f6_ll, lat, lon)

244 elif raiser or (raiser is None and self._raiser):

245 raise self._outsidError(lat, lon)

246 else:

247 x = y = z = N = NAN

248 return AxyzNgeoid4Tuple(x, y, z, N, name=name)

249

250 def _forwardOn3(self, unused, lat, lon): # must be overloaded!

251 self._notOverloaded(lat, lon) # PYCHOK no cover

252

253 def isinside(self, lat, lon, eps=0):

254 '''Is geodetic C{B{lat}} and C{B{lon}} inside L{region4}?

255

256 @arg lat: Latitude (C{degrees}, geodetic).

257 @arg lon: Longitude (C{degrees}, geodetic).

258 @kwarg eps: Over-/undersize L{region4} (C{degrees}).

259

260 @return: C{None} if B{C{lat}} or B{C{lon}} is NAN, C{False}

261 if outside L{region4}, C{True} otherwise.

262

263 @see: Method C{Bounds4Tuple.isinside}.

264 '''

265 lat, lon, _NAN, _, _ = self._LatLon5(lat, lon, False)

266 return None if _NAN else _isinside(lat, lon, Degrees(eps=eps),

267 _region4ax)

268

269 @property_RO

270 def kind(self):

271 '''Get the interpolation kind (C{int} or C{None}).

272 '''

273 return self._kind

274

275 def _LatLon5(self, lat, lon, raiser=None, name=_forward_):

276 # return lat, lon, ... if non-NAN

277 lat, lon = Lat(lat, clip=0), Lon(lon, clip=0)

278 _NAN = _isNAN(lat) or _isNAN(lon)

279 return lat, lon, _NAN, raiser, name

280

281 @property_RO

282 def lon0(self):

283 '''Get the triaxial's prime-meridian rotation (C{degrees}).

284 '''

285 return self.triaxial.lon0

286

287 @property

288 def onEPS(self):

289 '''Get the default L{sideOf} tolerance (C{meter}, I{squared}).

290 '''

291 return self._onEPS

292

293 @onEPS.setter # PYCHOK setter!

294 def onEPS(self, eps):

295 '''Set the default L{sideOf} tolerance (C{meter}, I{squared}).

296 '''

297 self._onEPS = max(float(eps), _0_0)

298

299 def _outsidError(self, llxyz, region):

300 # format an AxQGError for C{llxyz} outside C{region*}

301 return AxQGError(llxyz, txt=_outside__ + region.toRepr())

302

303 @property

304 def raiser(self):

305 '''Do points outside L{region4} cause an C{AxQGError}?

306 '''

307 return self._raiser

308

309 @raiser.setter # PYCHOK setter!

310 def raiser(self, raiser):

311 '''Use C{True} to throw an C{AxQGError} for points outside L{region4}.

312 '''

313 self._raiser = bool(raiser)

314

315 def region4(self, **unused):

316 '''Get the South, West, North and East bounds of the C{axQG} region as

317 L{Bounds4Tuple}C{(latS, lonW, latN, lonE)}.

318 '''

319 return _region4ax

320

321 def reverse(self, x, y, z, H=0, **raiser_name):

322 '''Convert a cartesian C{x}, C{y}, C{z} and orthometric height B{C{H}}

323 point to geodetic C{lat-}, C{longitude} and ellipsoidal C{height}.

324

325 @arg x: X coordinate (C{meter}, cartesian).

326 @arg y: Y coordinate (C{meter}, cartesian).

327 @arg z: Z coordinate (C{meter}, cartesian).

328 @kwarg H: The (orthometric) height (C{meter}, conventionally) or

329 C{None} to ignore height interpolation.

330 @kwarg raiser_name: Use C{B{raiser}=True} to raise an L{AxQGError}

331 for points outside L{region4}, overriding property

332 C{raiser} (C{bool}) and optional C{B{name}='reverse'}

333 (C{str}).

334

335 @return: An L{AxQG8Tuple}C{(x, y, z, H, lat, lon, height, axQG)} with

336 geodetic C{lat} and C{lon} and (ellipsoidal) C{height} in

337 C{meter} or C{NAN} and C{axQG} is this C{Ax*QG} instance.

338

339 @raise AxQGError: If the point is not on the L{triaxial}'s surface

340 and property C{raiser is True} or keyword argument

341 C{B{raiser}=True}.

342

343 @note: Ellipsoidal height C{(h = H + N)} equals orthometric height C{H}

344 plus (hybrid quasi-) geoid height C{N}.

345 '''

346 x, y, z, _NAN, raiser, name = self._xyz6(x, y, z, **raiser_name)

347 if _NAN:

348 lat = lon = height = NAN

349 else:

350 lat, lon, height = self._reverse3(raiser, x, y, z, H)

351 return AxQG8Tuple(x, y, z, H, lat, lon, height, self, name=name)

352

353 def _reverse3(self, raiser, x, y, z, H): # in .__main__

354 # C{reverse} core, returning C{(lat, lon, height)}

355 lat, lon = self._reverse2(x, y, z)

356 height = NAN if H is None or _isNAN(H) else (

357 Height(H=H) + self._axN(lat, lon, raiser))

358 return lat, lon, height

359

360 def _reverse2(self, x, y, z): # must be overloaded!

361 self._notOverloaded(x, y, z) # PYCHOK no cover

362

363 def reverseOn(self, x, y, z, H=0, normal=True, **name):

364 '''Project cartesian C{x}, C{y}, C{z} onto this triaxial's surface.

365

366 @arg x: X coordinate (C{meter}, cartesian).

367 @arg y: Y coordinate (C{meter}, cartesian).

368 @arg z: Z coordinate (C{meter}, cartesian).

369 @kwarg H: The (orthometric) height (C{meter}, conventionally) or

370 C{None} to ignore height interpolation.

371 @kwarg normal: If C{True}, the projection is C{perpendicular} to

372 the surface, otherwise C{radial} to the center of

373 this triaxial (C{bool}).

374 @kwarg name: Optional C{B{name}='reverseOn'} (C{str}).

375

376 @return: An L{AxQG8Tuple}C{(x, y, z, H, lat, lon, height, axQG)} with

377 cartesian C{x}, C{y}, C{z} I{on this triaxial's surface},

378 geodetic C{lat} and C{lon} and (ellipsoidal) C{height} in

379 C{meter} and C{axQG} is this C{Ax*QG} instance.

380 '''

381 t = self.triaxial.forwardCartesian(x, y, z, normal=normal)

382 n = _xkwds_get(name, name=typename(self.reverseOn))

383 return self.reverse(t.x, t.y, t.z, H=H, raiser=False, name=n)

384

385 def sideOf(self, x, y, z, **eps):

386 '''Is a cartesian on, above or below this triaxial's surface?

387

388 @arg x: X coordinate (C{meter}, cartesian).

389 @arg y: Y coordinate (C{meter}, cartesian).

390 @arg z: Z coordinate (C{meter}, cartesian).

391 @kwarg eps: Optional on-surface tolerance (C{meter}, squared),

392 overriding default L{onEPS}.

393

394 @return: Signed, radial distance to this triaxial's surface

395 (C{meter} I{squared}), C{INT0} if within tolerance

396 B{C{eps}}, positive if outside or negative if inside

397 this triaxial.

398 '''

399 eps = _xkwds_get(eps, eps=self.onEPS)

400 return self.triaxial.sideOf(x, y, z, eps=eps)

401

402 def _sideOfError(self, xyz, s2):

403 # format an AxQGError for C{xyz} not on this triaxial

404 n = typename(self.sideOf)

405 s = 'in' if s2 < 0 else 'out'

406 t = '%s (%.3f) %sside, not on %r' % (n, s2, s, self.triaxial)

407 return AxQGError(xyz, txt=t)

408

409 @property_RO

410 def smooth(self):

411 '''Get the smoothing factor (C{int} or C{None}).

412 '''

413 return self._smooth

414

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

416 '''Return this C{Ax*QG} instance as a string.

417

418 @kwarg prec: Precision, number of decimal digits (C{int}, 0..9).

419

420 @return: This C{Ax*QG} (C{str}).

421 '''

422 return self.attrs(_name_, 'ellipsoid', 'kind', 'smooth', 'raiser', Nones=False, prec=prec)

423

424 @property_RO

425 def triaxial(self):

426 '''Get the C{WGS84} or C{EGM2008} reference triaxial (L{TriAxial}).

427 '''

428 return self._triaxial

429

430 def unrotate(self, x, y, z=INT0, lon=None):

431 '''Reverse a cartesian to C{Earth-Centered, Earth-Fixed (ECEF)} by this

432 triaxial's prime-meridian rotation C{lon0}.

433

434 @arg x: Rotated X coordinate (C{meter}, cartesian).

435 @arg y: Rotated Y coordinate (C{meter}, cartesian).

436 @arg z: Rotated Z coordinate (C{meter}, cartesian).

437 @kwarg lon: Optional rotation (C{Ang}, C{Degrees}, C{degrees} or C{str}),

438 overriding this triaxial's prime-meridian rotation L{lon0

439 <_AxQGbase.lon0>}. Negative B{C{lon}} rotates clockwise,

440 positive counter-clockwise.

441

442 @return: L{Vector3Tuple}C{(x, y, z)} with C{x} and C{y} reversed to C{ECEF}.

443

444 @see: B{Geoid grid file format} in U{Supplementary File 3

445 <https://link.Springer.com/article/10.1007/s00190-023-01717-1#Sec21>}.

446 '''

447

448 A = self.triaxial.Lon0 if lon is None else (lon if isinstance(lon, Ang) else

449 Ang(Lon(lon), unit=Degrees))

450 if A.degrees0:

451 s, c = A.sc2

452 x, y = fdot_(x, c, -y, s), \

453 fdot_(x, s, y, c)

454 return Vector3Tuple(x, y, z)

455

456 def _xyz6(self, x, y, z, raiser=None, name=_reverse_):

457 # return x, y, z, ... if non-NAN and on triaxial's surface

458 x, y, z = t = Meter(x=x), Meter(y=y), Meter(z=z)

459 if _isNAN(x) or _isNAN(y) or _isNAN(z):

460 _NAN = True

461 else:

462 s2 = self.triaxial.sideOf(*t, eps=self.onEPS)

463 if s2 and (raiser or (raiser is None and self.raiser)):

464 raise self._sideOfError(t, s2)

465 _NAN = bool(s2)

466 return x, y, z, _NAN, raiser, name

467

468

469class Ax2QG(_AxQGbase):

470 '''Biaxial C{WGS84} transformer.

471 '''

472 if _FOR_DOCS:

473 __init__ = _AxQGbase.__init__

474

475 @property_ROver

476 def _Ax_grid(self): # load the _Ax2 grid, I{once}

477 try: # also if ax_grids.zip is unzipped

478 from pyaxqg.ax_grids import ax2grid

479 except ImportError:

480 _import_pyaxqg_ax_grids()

481 from pyaxqg.ax_grids import ax2grid

482 return self._Ax_G(ax2grid._Ax2)._Nterpolate

483

484 @property_ROver

485 def _ecef(self):

486 from pygeodesy import EcefKarney

487 return EcefKarney() # WGS84

488

489 @property_RO

490 def _ellipsoid(self):

491 return self._ecef.ellipsoid

492

493 def _forward3(self, lat, lon):

494 # geodetic C{{lat, lon)} to cartesian C{(x, y, z)}

495 return self._ecef.forward(lat, lon, 0).xyz

496

497 def _forwardOn3(self, unused, lat, lon):

498 # geodetic C{{lat, lon)} to C{(x, y, z)} on triaxial

499 t = self._triaxial.forward(lat, lon) # unit=Degrees

500 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

501 return t.xyz

502

503 def _reverse2(self, x, y, z):

504 # cartesian C{(x, y, z)} to geodetic C{{lat, lon)}

505 return self._ecef.reverse(x, y, z).latlon

506

507 @property_ROver

508 def _triaxial(self):

509 E = self.ellipsoid

510 return TriAxial(E.name, _0_0, E.a, E.a, E.b)

511

512

513class Ax3QG(_AxQGbase):

514 '''Triaxial C{EGM2008} transformer.

515 '''

516 @property_ROver

517 def _Ax_grid(self): # load the _Ax3 grid, I{once}

518 try: # also if ax_grids.zip is unzipped

519 from pyaxqg.ax_grids import ax3grid

520 except ImportError:

521 _import_pyaxqg_ax_grids()

522 from pyaxqg.ax_grids import ax3grid

523 return self._Ax_G(ax3grid._Ax3)._Nterpolate

524

525 @property_ROver

526 def _AxXgrid(self): # load the _AxX grid, I{once}

527 from pyaxqg.ax_grids import axXgrid

528 return self._Ax_G(axXgrid._AxX)._Nterpolate

529

530 @property_ROver

531 def _AxYgrid(self): # load the _AxY grid, I{once}

532 from pyaxqg.ax_grids import axYgrid

533 return self._Ax_G(axYgrid._AxY)._Nterpolate

534

535 @property_ROver

536 def _AxZgrid(self): # load the _AxZ grid, I{once}

537 from pyaxqg.ax_grids import axZgrid

538 return self._Ax_G(axZgrid._AxZ)._Nterpolate

539

540 @property_RO

541 def _ellipsoid(self):

542 return self._triaxial

543

544 def _forward3(self, lat, lon):

545 # geodetic C{{lat, lon)} to cartesian C{(x, y, z)}

546 t = self._triaxial.forward(lat, lon) # unit=Degrees

547 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

548 return t.xyz

549

550 def _forwardOn3(self, c_f_N_f6_ll, *unused): # PYCHOK signature

551 # interpolate C{x}, C{y}, C{z} B{after} C{N} imports ax_grids!

552 return (self._AxXgrid(*c_f_N_f6_ll),

553 self._AxYgrid(*c_f_N_f6_ll),

554 self._AxZgrid(*c_f_N_f6_ll))

555

556 def _reverse2(self, x, y, z):

557 # cartesian C{(x, y, z)} to geodetic C{{lat, lon)}

558 t = self._triaxial.reverse(x, y, z)

559 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

560 return t.lat, t.lon

561

562 @property_ROver

563 def _triaxial(self): # "A reference triaxial ellipsoid of the Earth"

564 # <https://link.Springer.com/article/10.1007/s00190-023-01717-1> ...

565 lon0 = -14.92850851 # TriAxial._lon0WGS84_3 # .triaxials.bases._Triaxial3Base

566 return TriAxial('EGM2008', lon0, 6378171.860779762, # ... Table 5

567 6378102.104632902,

568 6356752.334340346)

569

570

571class Ax5QG(Ax3QG):

572 '''Triaxial C{EGM2008} transformer using C{SciPy} bivariate spline interpolations.

573 '''

574 def __init__(self, kind=3, smooth=0, **raiser_name):

575 '''New C{Ax5QG} transformer instance.

576

577 @kwarg kind: C{scipy.interpolate} order (C{int}, -1, -3, -5 or 1..5),

578 see class U{pygeodesy.GeoidQuasi<https://mrJean1.GitHub.io/

579 PyGeodesy/docs/pygeodesy.geoids.GeoidQuasi-class.html>} for

580 the bivariate spline kinds and further information.

581 @kwarg smooth: Spline smoothing factor for C{B{kind}=1..5} only (C{float}).

582 @kwarg raiser_name: See L{Ax3QG<_AxQGbase.__init__>}.

583 '''

584 self._kind = kind # int -5, -3, -2, -1, 1..5

585 self._smooth = smooth # float or 0 (or None}

586 Ax3QG.__init__(self, **raiser_name)

587

588 @property_ROver

589 def _Ax_grid(self): # load the _Ax3 grid, I{once}

590 try: # also if ax_grids.zip is unzipped

591 from pyaxqg.ax_grids import ax3grid

592 except ImportError:

593 _import_pyaxqg_ax_grids()

594 from pyaxqg.ax_grids import ax3grid

595 return self._toNterpolate(ax3grid._Ax3)

596

597 @property_ROver

598 def _AxXgrid(self): # load the _AxX grid, I{once}

599 from pyaxqg.ax_grids import axXgrid

600 return self._toNterpolate(axXgrid._AxX)

601

602 @property_ROver

603 def _AxYgrid(self): # load the _AxY grid, I{once}

604 from pyaxqg.ax_grids import axYgrid

605 return self._toNterpolate(axYgrid._AxY)

606

607 @property_ROver

608 def _AxZgrid(self): # load the _AxZ grid, I{once}

609 from pyaxqg.ax_grids import axZgrid

610 return self._toNterpolate(axZgrid._AxZ)

611

612 def _c_f_N_f6_ll(self, lat, lon):

613 # pass C{(lat, lon)} to ._Nterpolate

614 return lat, lon

615

616 @property_ROver

617 def _GeoidQuasi(self): # lazily import GeoidQuasi, I{once}

618 from pygeodesy import GeoidQuasi

619 return GeoidQuasi

620

621 def _toNterpolate(self, ax):

622 '''(INTERNAL) Set up C{GeoidQuasi._Nterpolate(lat, lon)} for C{ax} grid.

623 '''

624 ax = self._Ax_G(ax)

625 gq = self._GeoidQuasi(ax.transpose, dtype=ax.dtype, kind=self.kind, smooth=self.smooth)

626 # assert gq.dtype == ax.dtype

627 # assert gq.shape == (_S2N, _W2E) # reversed(self._Ax_G2)

628 # assert gq.lowerleft()[:2] == _region4ax[:2]

629 # assert gq.upperright()[:2] == _region4ax[2:]

630 return gq._Nterpolate # (lat, lon)

631

632

633class TriAxial(Triaxial3):

634 '''Ordered C{pygeodesy.Triaxial3} for C{Ax*QG}.

635 '''

636 def __init__(self, name, lon0, *abc): # PYCHOK signature

637 '''New L{TriAxial} named B{C{name}} (C{str}), prime-meridian

638 rotated to B{C{lon0}} (C{degrees}) and I{ordered} axes

639 B{C{a}}, B{C{b}} and B{C{c}} (C{meter}).

640 '''

641 Triaxial3.__init__(self, *abc, name=name)

642 self.Lon0 = lon0 # type(Lon) is Ang

643

644 def __repr__(self):

645 '''Default C{repr(self)}.

646 '''

647 t = self.toRepr(terse=-5)

648 return t.replace(' Lon0=', ' lon0=')

649

650 @Property_RO

651 def lon0(self):

652 '''Get the prime-meridian rotation (C{degrees}).

653 '''

654 return Degrees(lon0=self.Lon0.degrees)

655

656

657def _c_f_N_f3(*deg_SWD):

658 # return int(ceil) and int(floor) of Normalized

659 # and (Normalized less floor) of C{deg} degrees

660 N = _degN(*deg_SWD)

661 # assert N >= 0, N

662 f = floor(N)

663 return int(ceil(N)), int(f), (N - f)

664

665

666def _darray(mx):

667 # meridian C{mx}, col-ordered _AxX/Y/Z grid

668 m = _array('d', map(float, mx.split()))

669 assert len(m) == _nS2N

670 return m

671

672

673def _degN(deg, degSW, degD):

674 # return C{deg} Normalized

675 d = float(deg - degSW)

676 if degD is not _1_Degree:

677 d = d / degD # /= chokes PyChecker

678 return d

679

680

681def _farray(mx):

682 # meridian C{mx}, col-ordered _Ax2/3 grid

683 m = _array('f', map(float, mx.split()))

684 assert len(m) == _nS2N

685 return m

686

687

688def _import_pyaxqg_ax_grids():

689 # set sys.modules['pyaxqg.ax_grids'] to ax_grids

690 from pyaxqg import _sys_modules_pyaxqg

691 ax_grids = _sys_modules_pyaxqg('ax_grids')

692 if not ax_grids:

693 raise AxQGError(_sys_modules_pyaxqg='ax_grids', txt=str(ax_grids))

694

695

696if _FOR_DOCS: # force epydoc to document all ...

697 for A in (Ax2QG, Ax3QG, Ax5QG): # ... public methods

698 A.axN = _AxQGbase.axN

699 A.axN3 = _AxQGbase.axN3

700 A.ellipsoid = _AxQGbase.ellipsoid

701 A.forward = _AxQGbase.forward

702 A.forwardOn = _AxQGbase.forwardOn

703 A.isinside = _AxQGbase.isinside

704 A.onEPS = _AxQGbase.onEPS

705 A.region4 = _AxQGbase.region4

706 A.reverse = _AxQGbase.reverse

707 A.reverseOn = _AxQGbase.reverseOn

708 A.sideOf = _AxQGbase.sideOf

709 A.triaxial = _AxQGbase.triaxial

710 A.unrotate = _AxQGbase.unrotate

711

712__all__ += _ALL_DOCS(_AxQGbase)

713__all__ += _ALL_OTHER(Ax2QG, Ax3QG, Ax5QG, TriAxial,

714 Bounds4Tuple, LatLonNgeoid3Tuple)

715del _ALL_DOCS, _ALL_OTHER, _FOR_DOCS

716

717# **) MIT License

718#

720#

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

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

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

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

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

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

727#

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

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

730#

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

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

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

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

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

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

737# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyaxqg / axqgs.py: 97%

295 statements