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>}.

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

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

11cartesian to geodetic coordinates and orthometric to ellipsoidal height.

12'''

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

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

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

17 _0_0, _1_0, _90_0, _180_0, _isNAN,

18 _ALL_DOCS, _ALL_OTHER, _FOR_DOCS,

19 _xkwds, _xkwds_get,

20 _isinside, _name_, _NamedBase, Vector3Tuple)

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

22 Bounds4Tuple, LatLonNgeoid3Tuple, # namedTuples

23 Property_RO, property_RO, property_ROver, # props

24 Triaxial3, # LLK as _LLK, # triaxials

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

27from array import array as _array

28from math import ceil, floor

30__all__ = ()

31__version__ = '26.08.23'

33_1_Degree = Degrees(_1_0)

34_forward_ = 'forward'

35_outside__ = 'outside '

36_region4ax = Bounds4Tuple(-_90_0, -_180_0,

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

38_reverse_ = 'reverse'

39_S2N = 181

40_W2E = 361 # PYCHOK in .ax*grid

43class _AxQGbase(_NamedBase):

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

45 '''

46 _Ax_grid = None # overloaded with _Ax2 or _Ax3

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

48 _latD = \

49 _lonD = _1_Degree

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

51 _raiser = False

52 _triaxial = None # overloaded with TriAxial

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

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

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

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

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

60 '''

61 if raiser:

62 self.raiser = True

63 if name:

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

66 def _Ax_assert(self, _Ax):

67 # assert len(_Ax) == _W2E # in ax*grid.py

68 # assert all(len(m) == _S2N for m in _Ax) # in _meridian

70 S, W, N, E = _region4ax

71 assert int(_degN(E, W, self._lonD) + _1_0) == len(_Ax) == _W2E

72 assert int(_degN(N, S, self._latD) + _1_0) == len(_Ax[0]) == _S2N

73 return _Ax

75 def axN(self, lat, lon):

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

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

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

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

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

83 '''

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

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

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

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

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

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

91 c_f_N_f6 = self._c_f_N_f6(lat, lon)

92 N = _bilinear(self._Ax_grid, *c_f_N_f6)

93 N = Height(N=N)

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

95 raise self._outsidError(lat, lon)

96 else:

97 N = NAN

98 return N

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

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

102

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

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

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

106

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

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

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

110 '''

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

112

113 def _c_f_N_f6(self, lat, lon):

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

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

116 S, W, _, _ = _region4ax

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

118 _c_f_N_f3(lon, W, self._lonD)

119

120 @property_RO

121 def ellipsoid(self):

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

123 '''

124 return self._ellipsoid

125

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

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

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

129 orthometric height C{H}.

130

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

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

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

134 or C{None} to ignore height interpolation.

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

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

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

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

139

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

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

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

143

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

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

146

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

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

149 '''

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

151 if _NAN:

152 x = y = z = H = NAN

153 else:

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

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

156

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

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

159

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

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

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

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

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

165 return x, y, z, H

166

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

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

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

170

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

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

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

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

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

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

177

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

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

180 or C{NAN}.

181

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

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

184

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

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

187 '''

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

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

190 if _NAN:

191 x = y = z = N = NAN

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

193 c_f_N_f6 = self._c_f_N_f6(lat, lon)

194 N = _bilinear(self._Ax_grid, *c_f_N_f6) # imports ax_grids

195 x, y, z = self._forwardOn3(c_f_N_f6, lat, lon)

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

197 raise self._outsidError(lat, lon)

198 else:

199 x = y = z = N = NAN

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

201

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

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

204

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

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

207

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

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

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

211

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

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

214

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

216 '''

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

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

219 _region4ax)

220

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

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

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

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

225 return lat, lon, _NAN, raiser, name

226

227 @property_RO

228 def lon0(self):

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

230 '''

231 return self.triaxial.lon0

232

233 @property

234 def onEPS(self):

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

236 '''

237 return self._onEPS

238

239 @onEPS.setter # PYCHOK setter!

240 def onEPS(self, eps):

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

242 '''

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

244

245 def _outsidError(self, llxyz, region):

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

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

248

249 @property

250 def raiser(self):

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

252 '''

253 return self._raiser

254

255 @raiser.setter # PYCHOK setter!

256 def raiser(self, raiser):

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

258 '''

259 self._raiser = bool(raiser)

260

261 def region4(self, **unused):

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

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

264 '''

265 return _region4ax

266

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

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

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

270

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

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

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

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

275 C{None} to ignore height interpolation.

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

277 for points outside L{region4}, overriding property

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

279 (C{str}).

280

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

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

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

284

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

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

287 C{B{raiser}=True}.

288

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

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

291 '''

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

293 if _NAN:

294 lat = lon = height = NAN

295 else:

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

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

298

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

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

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

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

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

304 return lat, lon, height

305

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

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

308

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

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

311

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

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

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

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

316 C{None} to ignore height interpolation.

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

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

319 this triaxial (C{bool}).

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

321

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

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

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

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

326 '''

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

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

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

330

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

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

333

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

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

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

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

338 overriding default L{onEPS}.

339

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

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

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

343 this triaxial.

344 '''

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

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

347

348 def _sideOfError(self, xyz, s2):

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

350 n = typename(self.sideOf)

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

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

353 return AxQGError(xyz, txt=t)

354

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

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

357

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

359

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

361 '''

362 return self.attrs(_name_, 'ellipsoid', 'raiser', prec=prec)

363

364 @property_RO

365 def triaxial(self):

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

367 '''

368 return self._triaxial

369

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

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

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

373

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

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

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

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

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

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

380 positive counter-clockwise.

381

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

383

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

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

386 '''

387

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

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

390 if A.degrees0:

391 s, c = A.sc2

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

393 fdot_(x, s, y, c)

394 return Vector3Tuple(x, y, z)

395

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

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

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

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

400 _NAN = True

401 else:

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

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

404 raise self._sideOfError(t, s2)

405 _NAN = bool(s2)

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

407

408

409class Ax2QG(_AxQGbase):

410 '''Biaxial C{WGS84} transformer.

411 '''

412 @property_ROver

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

414 try: # also if ax_grids.zip is unzipped

415 from pyaxqg.ax_grids import ax2grid

416 except ImportError:

417 _import_pyaxqg_ax_grids()

418 from pyaxqg.ax_grids import ax2grid

419 return self._Ax_assert(ax2grid._Ax2)

420

421 @property_ROver

422 def _ecef(self):

423 from pygeodesy import EcefKarney

424 return EcefKarney() # WGS84

425

426 @property_RO

427 def _ellipsoid(self):

428 return self._ecef.ellipsoid

429

430 def _forward3(self, lat, lon):

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

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

433

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

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

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

437 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

438 return t.xyz

439

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

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

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

443

444 @property_ROver

445 def _triaxial(self):

446 E = self.ellipsoid

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

448

449

450class Ax3QG(_AxQGbase):

451 '''Triaxial C{EGM2008} transformer.

452 '''

453 @property_ROver

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

455 try: # also if ax_grids.zip is unzipped

456 from pyaxqg.ax_grids import ax3grid

457 except ImportError:

458 _import_pyaxqg_ax_grids()

459 from pyaxqg.ax_grids import ax3grid

460 return self._Ax_assert(ax3grid._Ax3)

461

462 @property_ROver

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

464 from pyaxqg.ax_grids import axXgrid

465 return self._Ax_assert(axXgrid._AxX)

466

467 @property_ROver

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

469 from pyaxqg.ax_grids import axYgrid

470 return self._Ax_assert(axYgrid._AxY)

471

472 @property_ROver

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

474 from pyaxqg.ax_grids import axZgrid

475 return self._Ax_assert(axZgrid._AxZ)

476

477 @property_RO

478 def _ellipsoid(self):

479 return self._triaxial

480

481 def _forward3(self, lat, lon):

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

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

484 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

485 return t.xyz

486

487 def _forwardOn3(self, c_f_N_f6, *unused): # PYCHOK signature

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

489 return (_bilinear(self._AxXgrid, *c_f_N_f6),

490 _bilinear(self._AxYgrid, *c_f_N_f6),

491 _bilinear(self._AxZgrid, *c_f_N_f6))

492

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

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

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

496 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

497 return t.lat, t.lon

498

499 @property_ROver

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

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

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

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

504 6378102.104632902,

505 6356752.334340346)

506

507

508class TriAxial(Triaxial3):

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

510 '''

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

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

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

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

515 '''

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

517 self.Lon0 = lon0 # Lon is Ang

518

519 def __repr__(self):

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

521 '''

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

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

524

525 @Property_RO

526 def lon0(self):

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

528 '''

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

530

531

532def _bilinear(Ax, c_latI, f_latI, latN_f,

533 c_lonI, f_lonI, lonN_f):

534 # interpolate an C{Ax[_W2E][_S2N]} grid, col-ordered

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

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

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

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

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

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

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

542 else:

543 x = Ax[c_lonI][c_latI]

544 return x

545

546

547def _c_f_N_f3(*deg_SWD):

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

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

550 N = _degN(*deg_SWD)

551 # assert N >= 0, N

552 f = floor(N)

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

554

555

556def _darray(mx):

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

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

559 assert len(m) == _S2N

560 return m

561

562

563def _degN(deg, degSW, degD):

564 # return C{deg} Normalized

565 d = float(deg - degSW)

566 if degD is not _1_Degree:

567 d = d / degD # /= chokes PyChecker

568 return d

569

570

571def _import_pyaxqg_ax_grids():

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

573 from pyaxqg import _sys_modules_pyaxqg

574 ax_grids = _sys_modules_pyaxqg('ax_grids')

575 if not ax_grids:

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

577

578

579def _farray(mx):

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

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

582 assert len(m) == _S2N

583 return m

584

585

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

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

588 A.__init__ = _AxQGbase.__init__

589 A.axN = _AxQGbase.axN

590 A.axN3 = _AxQGbase.axN3

591 A.ellipsoid = _AxQGbase.ellipsoid

592 A.forward = _AxQGbase.forward

593 A.forwardOn = _AxQGbase.forwardOn

594 A.isinside = _AxQGbase.isinside

595 A.onEPS = _AxQGbase.onEPS

596 A.region4 = _AxQGbase.region4

597 A.reverse = _AxQGbase.reverse

598 A.reverseOn = _AxQGbase.reverseOn

599 A.sideOf = _AxQGbase.sideOf

600 A.triaxial = _AxQGbase.triaxial

601 A.unrotate = _AxQGbase.unrotate

602

603__all__ += _ALL_DOCS(_AxQGbase)

604__all__ += _ALL_OTHER(Ax2QG, Ax3QG, TriAxial, Bounds4Tuple, LatLonNgeoid3Tuple)

605del _ALL_DOCS, _ALL_OTHER

606

607# **) MIT License

608#

610#

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

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

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

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

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

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

617#

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

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

620#

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

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

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

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

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

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

627# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyaxqg / axqgs.py: 97%

242 statements