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 (geocentric) cartesian X, Y, Z and orthometric height and a C{reverse} method for

11converting cartesian 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, 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 _farray

28from math import ceil, floor

30__all__ = ()

31__version__ = '26.08.22'

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 _forward4(self, raiser, lat, lon, height): # in .__main__

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

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

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

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

162 return x, y, z, H

163

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

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

166

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

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

169

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

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

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

173

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

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

176

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

178 '''

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

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

181 _region4ax)

182

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

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

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

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

187 return lat, lon, _NAN, raiser, name

188

189 @property_RO

190 def lon0(self):

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

192 '''

193 return self.triaxial.lon0

194

195 @property

196 def onEPS(self):

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

198 '''

199 return self._onEPS

200

201 @onEPS.setter # PYCHOK setter!

202 def onEPS(self, eps):

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

204 '''

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

206

207 def _outsidError(self, llxyz, region):

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

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

210

211 @property

212 def raiser(self):

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

214 '''

215 return self._raiser

216

217 @raiser.setter # PYCHOK setter!

218 def raiser(self, raiser):

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

220 '''

221 self._raiser = bool(raiser)

222

223 def region4(self, **unused):

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

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

226 '''

227 return _region4ax

228

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

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

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

232

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

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

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

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

237 C{None} to ignore height interpolation.

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

239 for points outside L{region4}, overriding property

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

241 (C{str}).

242

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

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

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

246

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

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

249 C{B{raiser}=True}.

250

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

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

253 '''

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

255 if _NAN:

256 lat = lon = height = NAN

257 else:

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

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

260

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

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

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

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

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

266 return lat, lon, height

267

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

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

270

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

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

273

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

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

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

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

278 C{None} to ignore height interpolation.

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

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

281 this triaxial (C{bool}).

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

283

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

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

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

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

288 '''

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

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

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

292

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

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

295

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

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

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

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

300 overriding default L{onEPS}.

301

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

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

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

305 this triaxial.

306 '''

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

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

309

310 def _sideOfError(self, xyz, s2):

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

312 n = typename(self.sideOf)

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

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

315 return AxQGError(xyz, txt=t)

316

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

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

319

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

321

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

323 '''

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

325

326 @property_RO

327 def triaxial(self):

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

329 '''

330 return self._triaxial

331

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

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

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

335

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

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

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

339 @kwarg lon: Optional rotation (C{degrees}), overriding this triaxial's

340 prime-meridian rotation L{lon0<_AxQGbase.lon0>}.

341

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

343

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

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

346 '''

347 A = self.triaxial.Lon0 if lon is None else Ang(Lon(lon), unit=Degrees)

348 if A.degrees0:

349 s, c = A.sc2

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

351 fdot_(x, s, y, c)

352 return Vector3Tuple(x, y, z)

353

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

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

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

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

358 _NAN = True

359 else:

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

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

362 raise self._sideOfError(t, s2)

363 _NAN = bool(s2)

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

365

366

367class Ax2QG(_AxQGbase):

368 '''Biaxial C{WGS84} transformer.

369 '''

370 @property_ROver

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

372 try: # also if ax_grids.zip is unzipped

373 from pyaxqg.ax_grids import ax2grid

374 except ImportError:

375 _import_pyaxqg_ax_grids()

376 from pyaxqg.ax_grids import ax2grid

377 return self._Ax_assert(ax2grid._Ax2)

378

379 @property_ROver

380 def _ecef(self):

381 from pygeodesy import EcefKarney

382 return EcefKarney() # WGS84

383

384 @property_RO

385 def _ellipsoid(self):

386 return self._ecef.ellipsoid

387

388 def _forward3(self, lat, lon):

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

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

391

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

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

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

395

396 @property_ROver

397 def _triaxial(self):

398 E = self.ellipsoid

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

400

401

402class Ax3QG(_AxQGbase):

403 '''Triaxial C{EGM2008} transformer.

404 '''

405 @property_ROver

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

407 try: # also if ax_grids.zip is unzipped

408 from pyaxqg.ax_grids import ax3grid

409 except ImportError:

410 _import_pyaxqg_ax_grids()

411 from pyaxqg.ax_grids import ax3grid

412 return self._Ax_assert(ax3grid._Ax3)

413

414 @property_ROver

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

416 try: # also if ax_grids.zip is unzipped

417 from pyaxqg.ax_grids import axXgrid

418 except ImportError:

419 _import_pyaxqg_ax_grids()

420 from pyaxqg.ax_grids import axXgrid

421 return self._Ax_assert(axXgrid._AxX)

422

423 @property_ROver

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

425 try: # also if ax_grids.zip is unzipped

426 from pyaxqg.ax_grids import axYgrid

427 except ImportError:

428 _import_pyaxqg_ax_grids()

429 from pyaxqg.ax_grids import axYgrid

430 return self._Ax_assert(axYgrid._AxY)

431

432 @property_ROver

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

434 try: # also if ax_grids.zip is unzipped

435 from pyaxqg.ax_grids import axZgrid

436 except ImportError:

437 _import_pyaxqg_ax_grids()

438 from pyaxqg.ax_grids import axZgrid

439 return self._Ax_assert(axZgrid._AxZ)

440

441 @property_RO

442 def _ellipsoid(self):

443 return self._triaxial

444

445 def _forward3(self, lat, lon):

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

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

448 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

449 return t.xyz

450

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

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

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

454

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

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

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

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

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

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

461

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

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

464 or C{NAN}.

465

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

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

468

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

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

471 '''

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

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

474 if _NAN:

475 x = y = z = N = NAN

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

477 c_f_N_f6 = self._c_f_N_f6(lat, lon)

478 x = _bilinear(self._AxXgrid, *c_f_N_f6)

479 y = _bilinear(self._AxYgrid, *c_f_N_f6)

480 z = _bilinear(self._AxZgrid, *c_f_N_f6)

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

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

483 raise self._outsidError(lat, lon)

484 else:

485 x = y = z = N = NAN

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

487

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

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

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

491 # assert t.llk == _LLK_GEODETIC or _GEODETIC_LON0

492 return t.lat, t.lon

493

494 @property_ROver

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

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

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

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

499 6378102.104632902,

500 6356752.334340346)

501

502

503class TriAxial(Triaxial3):

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

505 '''

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

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

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

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

510 '''

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

512 self.Lon0 = lon0 # Lon is Ang

513

514 def __repr__(self):

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

516 '''

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

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

519

520 @Property_RO

521 def lon0(self):

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

523 '''

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

525

526

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

528 c_lonI, f_lonI, lonN_f):

529 # interpolate quasi-geoid height from Ax2/3

530 # grid, col-ordered [_W2E][_S2N] meridians

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

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

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

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

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

536 Nh = (ne * lonN_f + nw * lonN_f1) * latN_f + \

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

538 else:

539 Nh = Ax[c_lonI][c_latI]

540 return Nh

541

542

543def _c_f_N_f3(*deg_SWD):

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

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

546 N = _degN(*deg_SWD)

547 # assert N >= 0, N

548 f = floor(N)

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

550

551

552def _degN(deg, degSW, degD):

553 # return C{deg} Normalized

554 d = float(deg - degSW)

555 if degD is not _1_Degree:

556 d = d / degD # /= chokes PyChecker

557 return d

558

559

560def _import_pyaxqg_ax_grids():

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

562 from pyaxqg import _sys_modules_pyaxqg

563 ax_grids = _sys_modules_pyaxqg('ax_grids')

564 if not ax_grids:

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

566

567

568def _meridian(Mx):

569 # meridian C{Mx}, col-ordered Ax2/3 grid

570 m = _farray('f', map(float, Mx.split()))

571 assert len(m) == _S2N

572 return m

573

574

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

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

577 A.__init__ = _AxQGbase.__init__

578 A.axN = _AxQGbase.axN

579 A.axN3 = _AxQGbase.axN3

580 A.ellipsoid = _AxQGbase.ellipsoid

581 A.forward = _AxQGbase.forward

582 A.isinside = _AxQGbase.isinside

583 A.onEPS = _AxQGbase.onEPS

584 A.region4 = _AxQGbase.region4

585 A.reverse = _AxQGbase.reverse

586 A.reverseOn = _AxQGbase.reverseOn

587 A.sideOf = _AxQGbase.sideOf

588 A.triaxial = _AxQGbase.triaxial

589 A.unrotate = _AxQGbase.unrotate

590

591__all__ += _ALL_DOCS(_AxQGbase)

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

593del _ALL_DOCS, _ALL_OTHER

594

595# **) MIT License

596#

598#

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

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

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

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

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

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

605#

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

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

608#

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

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

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

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

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

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

615# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pyaxqg / axqgs.py: 97%

237 statements