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.08.28'

36_1_Degree = Degrees(_1_0)

37_forward_ = 'forward'

38_outside__ = 'outside '

39_region4ax = Bounds4Tuple(-_90_0, -_180_0,

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

41_reverse_ = 'reverse'

42_S2N = 181

43_W2E = 361 # PYCHOK in .ax*grid

46class _AxG(tuple):

47 '''(INTERNAL) [_W2E] * [_S2N] 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 [_W2E] lat [_S2N]

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[_W2E][_S2N]} grid, col-major ordered

65 Ax = self # see GeoidKarney, _Dotf and _Hornerf

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 [_W2E][_S2N] to [_S2N][_W2E]

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_assert(self, ax):

109 S_N, W_E = self._Ax_assert2

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

111 assert W_E == len(ax) == _W2E # in _AxG.__init__

112 assert S_N == len(ax[0]) == _S2N # in _darray, _farray

113 # assert all(len(m) == _S2N for m in ax) # in _darray, _farray

114 return _AxG(ax)

115

116 @Property_RO

117 def _Ax_assert2(self):

118 S, W, N, E = _region4ax

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

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

121

122 def axN(self, lat, lon):

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

124

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

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

127

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

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

130 '''

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

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

133

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

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

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

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

138 c_f_N_f6_ll = self._c_f_N_f6_ll(lat, lon)

139 N = self._Ax_grid(*c_f_N_f6_ll)

140 N = Height(N=N)

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

142 raise self._outsidError(lat, lon)

143 else:

144 N = NAN

145 return N

146

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

148 '''Interpolate the quasi-geoid height C{H} for a cartesian point.

149

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

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

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

153

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

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

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

157 '''

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

159

160 def _c_f_N_f6_ll(self, lat, lon):

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

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

163 S, W, _, _ = _region4ax

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

165 _c_f_N_f3(lon, W, self._lonD)

166

167 @property_RO

168 def ellipsoid(self):

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

170 '''

171 return self._ellipsoid

172

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

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

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

176 orthometric height C{H}.

177

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

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

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

181 or C{None} to ignore height interpolation.

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

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

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

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

186

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

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

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

190

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

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

193

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

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

196 '''

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

198 if _NAN:

199 x = y = z = H = NAN

200 else:

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

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

203

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

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

206

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

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

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

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

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

212 return x, y, z, H

213

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

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

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

217

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

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

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

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

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

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

224

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

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

227 or C{NAN}.

228

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

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

231

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

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

234 '''

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

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

237 if _NAN:

238 x = y = z = N = NAN

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

240 c_f_N_f6_ll = self._c_f_N_f6_ll(lat, lon)

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

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

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

244 raise self._outsidError(lat, lon)

245 else:

246 x = y = z = N = NAN

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

248

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

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

251

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

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

254

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

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

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

258

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

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

261

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

263 '''

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

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

266 _region4ax)

267

268 @property_RO

269 def kind(self):

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

271 '''

272 return self._kind

273

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

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

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

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

278 return lat, lon, _NAN, raiser, name

279

280 @property_RO

281 def lon0(self):

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

283 '''

284 return self.triaxial.lon0

285

286 @property

287 def onEPS(self):

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

289 '''

290 return self._onEPS

291

292 @onEPS.setter # PYCHOK setter!

293 def onEPS(self, eps):

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

295 '''

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

297

298 def _outsidError(self, llxyz, region):

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

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

301

302 @property

303 def raiser(self):

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

305 '''

306 return self._raiser

307

308 @raiser.setter # PYCHOK setter!

309 def raiser(self, raiser):

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

311 '''

312 self._raiser = bool(raiser)

313

314 def region4(self, **unused):

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

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

317 '''

318 return _region4ax

319

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

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

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

323

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

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

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

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

328 C{None} to ignore height interpolation.

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

330 for points outside L{region4}, overriding property

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

332 (C{str}).

333

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

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

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

337

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

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

340 C{B{raiser}=True}.

341

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

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

344 '''

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

346 if _NAN:

347 lat = lon = height = NAN

348 else:

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

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

351

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

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

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

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

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

357 return lat, lon, height

358

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

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

361

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

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

364

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

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

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

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

369 C{None} to ignore height interpolation.

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

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

372 this triaxial (C{bool}).

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

374

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

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

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

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

379 '''

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

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

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

383

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

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

386

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

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

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

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

391 overriding default L{onEPS}.

392

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

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

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

396 this triaxial.

397 '''

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

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

400

401 def _sideOfError(self, xyz, s2):

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

403 n = typename(self.sideOf)

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

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

406 return AxQGError(xyz, txt=t)

407

408 @property_RO

409 def smooth(self):

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

411 '''

412 return self._smooth

413

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

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

416

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

418

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

420 '''

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

422

423 @property_RO

424 def triaxial(self):

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

426 '''

427 return self._triaxial

428

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

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

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

432

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

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

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

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

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

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

439 positive counter-clockwise.

440

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

442

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

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

445 '''

446

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

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

449 if A.degrees0:

450 s, c = A.sc2

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

452 fdot_(x, s, y, c)

453 return Vector3Tuple(x, y, z)

454

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

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

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

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

459 _NAN = True

460 else:

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

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

463 raise self._sideOfError(t, s2)

464 _NAN = bool(s2)

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

466

467

468class Ax2QG(_AxQGbase):

469 '''Biaxial C{WGS84} transformer.

470 '''

471 if _FOR_DOCS:

472 __init__ = _AxQGbase.__init__

473

474 @property_ROver

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

476 try: # also if ax_grids.zip is unzipped

477 from pyaxqg.ax_grids import ax2grid

478 except ImportError:

479 _import_pyaxqg_ax_grids()

480 from pyaxqg.ax_grids import ax2grid

481 return self._Ax_assert(ax2grid._Ax2)._Nterpolate

482

483 @property_ROver

484 def _ecef(self):

485 from pygeodesy import EcefKarney

486 return EcefKarney() # WGS84

487

488 @property_RO

489 def _ellipsoid(self):

490 return self._ecef.ellipsoid

491

492 def _forward3(self, lat, lon):

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

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

495

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

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

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

499 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

500 return t.xyz

501

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

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

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

505

506 @property_ROver

507 def _triaxial(self):

508 E = self.ellipsoid

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

510

511

512class Ax3QG(_AxQGbase):

513 '''Triaxial C{EGM2008} transformer.

514 '''

515 @property_ROver

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

517 try: # also if ax_grids.zip is unzipped

518 from pyaxqg.ax_grids import ax3grid

519 except ImportError:

520 _import_pyaxqg_ax_grids()

521 from pyaxqg.ax_grids import ax3grid

522 return self._Ax_assert(ax3grid._Ax3)._Nterpolate

523

524 @property_ROver

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

526 from pyaxqg.ax_grids import axXgrid

527 return self._Ax_assert(axXgrid._AxX)._Nterpolate

528

529 @property_ROver

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

531 from pyaxqg.ax_grids import axYgrid

532 return self._Ax_assert(axYgrid._AxY)._Nterpolate

533

534 @property_ROver

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

536 from pyaxqg.ax_grids import axZgrid

537 return self._Ax_assert(axZgrid._AxZ)._Nterpolate

538

539 @property_RO

540 def _ellipsoid(self):

541 return self._triaxial

542

543 def _forward3(self, lat, lon):

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

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

546 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

547 return t.xyz

548

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

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

551 return (self._AxXgrid(*c_f_N_f6_ll),

552 self._AxYgrid(*c_f_N_f6_ll),

553 self._AxZgrid(*c_f_N_f6_ll))

554

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

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

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

558 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

559 return t.lat, t.lon

560

561 @property_ROver

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

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

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

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

566 6378102.104632902,

567 6356752.334340346)

568

569

570class Ax5QG(Ax3QG):

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

572 '''

573 _kind = \

574 _smooth = None # bilinear

575

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

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

578

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

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

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

582 the bivariate spline kinds and further information.

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

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

585 '''

586 Ax3QG.__init__(self, **raiser_name)

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

588 self._smooth = smooth # float or 0

589

590 @property_ROver

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

592 try: # also if ax_grids.zip is unzipped

593 from pyaxqg.ax_grids import ax3grid

594 except ImportError:

595 _import_pyaxqg_ax_grids()

596 from pyaxqg.ax_grids import ax3grid

597 return self._toNterpolate(ax3grid._Ax3)

598

599 @property_ROver

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

601 from pyaxqg.ax_grids import axXgrid

602 return self._toNterpolate(axXgrid._AxX)

603

604 @property_ROver

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

606 from pyaxqg.ax_grids import axYgrid

607 return self._toNterpolate(axYgrid._AxY)

608

609 @property_ROver

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

611 from pyaxqg.ax_grids import axZgrid

612 return self._toNterpolate(axZgrid._AxZ)

613

614 def _c_f_N_f6_ll(self, lat, lon):

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

616 return lat, lon

617

618 @property_ROver

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

620 from pygeodesy import GeoidQuasi

621 return GeoidQuasi

622

623 def _toNterpolate(self, ax):

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

625 '''

626 ax = self._Ax_assert(ax)

627 q = self._GeoidQuasi(ax.transpose, dtype=ax.dtype, kind=self.kind, smooth=self.smooth)

628 # assert q.dtype == ax.dtype

629 # assert q.shape == (_S2N, _W2E)

630 # assert q.lowerleft()[:2] == _region4ax[:2]

631 # assert q.upperright()[:2] == _region4ax[2:]

632 return q._Nterpolate # (lat, lon)

633

634

635class TriAxial(Triaxial3):

636 '''Ordered C{pygeodesy.Triaxial3} for L{Ax3QG} and L{Ax2QG}.

637 '''

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

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

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

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

642 '''

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

644 self.Lon0 = lon0 # Lon is Ang

645

646 def __repr__(self):

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

648 '''

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

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

651

652 @Property_RO

653 def lon0(self):

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

655 '''

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

657

658

659def _c_f_N_f3(*deg_SWD):

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

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

662 N = _degN(*deg_SWD)

663 # assert N >= 0, N

664 f = floor(N)

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

666

667

668def _darray(mx):

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

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

671 assert len(m) == _S2N

672 return m

673

674

675def _degN(deg, degSW, degD):

676 # return C{deg} Normalized

677 d = float(deg - degSW)

678 if degD is not _1_Degree:

679 d = d / degD # /= chokes PyChecker

680 return d

681

682

683def _farray(mx):

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

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

686 assert len(m) == _S2N

687 return m

688

689

690def _import_pyaxqg_ax_grids():

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

692 from pyaxqg import _sys_modules_pyaxqg

693 ax_grids = _sys_modules_pyaxqg('ax_grids')

694 if not ax_grids:

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

696

697

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

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

700 A.axN = _AxQGbase.axN

701 A.axN3 = _AxQGbase.axN3

702 A.ellipsoid = _AxQGbase.ellipsoid

703 A.forward = _AxQGbase.forward

704 A.forwardOn = _AxQGbase.forwardOn

705 A.isinside = _AxQGbase.isinside

706 A.onEPS = _AxQGbase.onEPS

707 A.region4 = _AxQGbase.region4

708 A.reverse = _AxQGbase.reverse

709 A.reverseOn = _AxQGbase.reverseOn

710 A.sideOf = _AxQGbase.sideOf

711 A.triaxial = _AxQGbase.triaxial

712 A.unrotate = _AxQGbase.unrotate

713

714__all__ += _ALL_DOCS(_AxQGbase)

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

716 Bounds4Tuple, LatLonNgeoid3Tuple)

717del _ALL_DOCS, _ALL_OTHER

718

719# **) MIT License

720#

722#

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

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

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

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

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

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

729#

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

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

732#

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

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

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

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

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

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

739# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyaxqg / axqgs.py: 97%

296 statements