Coverage for pygeodesy/geoids.py: 96%

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

4u'''Geoid models and geoid height interpolations.

6Classes L{GeoidEGM96}, L{GeoidG2012B}, L{GeoidKarney} and L{GeoidPGM} to

7interpolate the height of various U{geoid<https://WikiPedia.org/wiki/Geoid>}s

8at C{LatLon} locations or separate lat-/longitudes using various interpolation

9methods and C{geoid} model files.

11L{GeoidKarney} is a transcoding of I{Charles Karney}'s C++ class U{Geoid

12<https://GeographicLib.SourceForge.io/C++/doc/geoid.html>} to pure Python.

14The L{GeoidEGM96}, L{GeoidG2012B} and L{GeoidPGM} interpolators both depend on

15U{scipy<https://SciPy.org>} and U{numpy<https://PyPI.org/project/numpy>} and

16require those packages to be installed.

18In addition, each geoid interpolator needs C{grid knots} (down)loaded from

19a C{geoid} model file, I{specific to the interpolator}. More details below

20and in the documentation of the interpolator class. For each interpolator,

21there are several interpolation choices, like I{linear}, I{cubic}, etc.

23Typical usage

24=============

261. Choose an interpolator class L{GeoidEGM96}, L{GeoidG2012B}, L{GeoidKarney}

27or L{GeoidPGM} and download a C{geoid} model file, containing locations with

28known heights also referred to as the C{grid knots}. See the documentation of

29the interpolator class for references to available C{grid} models.

31C{>>> from pygeodesy import GeoidEGM96 as GeoidXyz # or GeoidG2012B, -Karney or -PGM}

332. Instantiate an interpolator with the C{geoid} model file and use keyword

34arguments to select different interpolation options

36C{>>> ginterpolator = GeoidXyz(geoid_model_file, **options)}

383. Get the interpolated geoid height of C{LatLon} location(s) with

40C{>>> ll = LatLon(1, 2, ...)}

42C{>>> h = ginterpolator(ll)}

44or

46C{>>> h1, h2, h3, ... = ginterpolator(ll1, ll2, ll3, ...)}

48or a list, tuple, generator, etc. of C{LatLon}s

50C{>>> hs = ginterpolator(lls)}

524. For separate lat- and longitudes invoke the C{height} method as

54C{>>> h = ginterpolator.height(lat, lon)}

56or as 2 lists, 2 tuples, etc.

58C{>>> hs = ginterpolator.height(lats, lons)}

60or for several positionals use the C{height_} method

62C{>>> h1, h2, ... = ginterpolator.height_(lat1, lon1, lat2, lon2, ...)}

645. An example is in U{issue #64<https://GitHub.com/mrJean1/PyGeodesy/issues/64>},

65courtesy of SBFRF.

67@note: Classes L{GeoidEGM96}, L{GeoidG2012B} and L{GeoidPGM} require both U{numpy

68 <https://PyPI.org/project/numpy>} and U{scipy<https://PyPI.org/project/scipy>}

69 to be installed.

71@note: Errors from C{scipy} are raised as L{SciPyError}s. Warnings issued by C{scipy} can

72 be thrown as L{SciPyWarning} exceptions, provided Python C{warnings} are filtered

73 accordingly, see L{SciPyWarning}.

75@see: I{Karney}'s U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/index.html>},

76 U{Geoid height<https://GeographicLib.SourceForge.io/C++/doc/geoid.html>} and U{Installing

77 the Geoid datasets<https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>},

78 World Geodetic System 1984 (WG84) and U{Earth Gravitational Model 96 (EGM96) Data and

79 Apps<https://earth-info.NGA.mil/index.php?dir=wgs84&action=wgs84>},

80 U{SciPy<https://docs.SciPy.org/doc/scipy/reference/interpolate.html>} interpolation

81 U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.

82 RectBivariateSpline.html>}, U{bisplrep/-ev<https://docs.scipy.org/doc/scipy/reference/generated/

83 scipy.interpolate.bisplrep.html>} and U{interp2d<https://docs.SciPy.org/doc/scipy/reference/

84 generated/scipy.interpolate.interp2d.html>}, functions L{elevations.elevation2} and

85 L{elevations.geoidHeight2}, U{I{Ellispoid vs Orthometric Elevations}<https://www.YouTube.com/

86 watch?v=dX6a6kCk3Po>} and U{6.22.1 Avoiding Pitfalls Related to Ellipsoid Height and Height

87 Above Mean Sea Level<https://Wiki.ROS.org/mavros>}.

88'''

89# make sure int/int division yields float quotient, see .basics

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

92from pygeodesy.basics import _float0d, _isin, len2, min2, isodd, _splituple

93from pygeodesy.constants import EPS, _float as _F, _1_0, _N_90_0, _180_0, \

94 _N_180_0, _360_0

95from pygeodesy.datums import Datums, _ellipsoidal_datum, _WGS84

96# from pygeodesy.dms import parseDMS2 # _MODS

97from pygeodesy.errors import _incompatible, LenError, RangeError, _SciPyIssue, \

98 _xkwds_pop2

99from pygeodesy.fmath import favg, frange, Fsum

100# from pygoedesy.formy import heightOrthometric # _MODS

101# from pygeodesy.fsums import Fsum # from .fmath

102from pygeodesy.heights import _as_llis2, _ascalar, _HeightBase, HeightError, _Wrap

103# from pygeodesy.internals import typename, _version2 # _MODS

104from pygeodesy.interns import NN, _COLONSPACE_, _COMMASPACE_, _DMAIN_, _E_, \

105 _height_, _in_, _kind_, _lat_, _lon_, _mean_, _N_, \

106 _n_a_, _numpy_, _on_, _outside_, _S_, _s_, _scipy_, \

107 _SPACE_, _stdev_, _tbd_, _W_, _width_, _4_

108from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, _FOR_DOCS

109from pygeodesy.named import _name__, _Named, _NamedTuple

110# from pygeodesy.namedTuples import LatLon3Tuple # _MODS

111from pygeodesy.props import deprecated_property_RO, Property_RO, property_RO, \

112 property_ROver

113from pygeodesy.streprs import attrs, Fmt, fstr, pairs

114from pygeodesy.units import Float_, Height, Int_, Lat, Lon

115# from pygeodesy.utily import _Wrap # from .heights

116

117from math import floor as _floor

118# import os as _os # _MODS

119# import os.path # _os.path

120from struct import calcsize as _calcsize, unpack as _unpack

121try:

122 from StringIO import StringIO as _BytesIO # reads bytes

123 _ub2str = str # PYCHOK convert bytes to str for egm*.pgm text

124

125except ImportError: # Python 3+

126 from io import BytesIO as _BytesIO # PYCHOK expected

127 from pygeodesy.basics import ub2str as _ub2str

128

129__all__ = _ALL_LAZY.geoids

130__version__ = '26.08.26'

131

132_assert_ = 'assert'

133_bHASH_ = b'#'

134_endian_ = 'endian'

135_format_ = '%s %r'

136_header_ = 'header'

137_intCs = {} # cache int value, del below

138_lli_ = 'lli'

139_os = _MODS.os

140_rb_ = 'rb'

141_supported_ = 'supported'

142

143if __debug__:

144 from pygeodesy.fmath import Fdot as _Dotf, Fhorner as _Hornerf

145

146else: # -OO ... runs GeoidKarney 8+X faster (w/o Fwelford)

147 from pygeodesy.fsums import _fsum

148 from operator import mul as _mul

149

150 class _GKsum(Fsum): # for GeoidKarney only

151

152 def __iadd__(self, other):

153 xs = other._ps if isinstance(other, _GKsum) else (other,)

154 self._ps_acc(self._ps, xs, up=False)

155 return self

156

157 def __imul__(self, x):

158 self._ps[:] = (x * p for p in self._ps)

159 return self

160

161 fmul = __imul__ # like .fsums.Fsum

162

163 def fsum(self): # PYCHOK signature

164 return _fsum(self._ps)

165

166 class _Dotf(_GKsum): # PYCHOK redef

167 '''Fast precision dot product M{sum(a[i] * b[i] for i=0..len(a))}

168 but for C{float} C{a} and C{b} only.

169 '''

170 def __init__(self, a, *b):

171 self._ps = self._ps_acc([], map(_mul, a, b), up=False)

172

173 class _Hornerf(_GKsum): # PYCHOK redef

174 '''Fast polynomial evaluation M{sum(cs[i] * x**i for i=0..len(cs))}

175 but for C{float} C{x} and C{cs} and C{incx=True} only.

176 '''

177 def __init__(self, x, *cs): # incx=True

178 self._ps = []

179 for c in reversed(cs): # multiply-accumulate

180 self *= x # ps[:] = (x * p for p in ps)

181 self += c # self._ps_acc(ps, (c,), up=False)

182

183

184class GeoidError(HeightError):

185 '''Geoid interpolator C{Geoid...} or interpolation issue.

186 '''

187 pass

188

189

190class _GeoidBase(_HeightBase):

191 '''(INTERNAL) Base class for C{Geoid...}s.

192 '''

193# _center = None # (lat, lon, height)

194 _cropped = None

195# _datum = _WGS84 # from _HeightBase

196 _egm = None # open C{egm*.pgm} geoid file

197 _endian = _tbd_

198 _Error = GeoidError # in ._HeightBase._as_lls, ...

199 _geoid = _n_a_

200 _fudge = 0 # for _iscipy._lon_hi

201 _hs_y_x = None # numpy 2darray, row-major order

202 _iscipy = True # scipy or Karney's interpolation

203 _kind = 3 # order for interp2d, RectBivariateSpline

204# _kmin = 2 # min number of knots

205 _mean = None # fixed in GeoidKarney

206 _nBytes = 0 # numpy size in bytes, float64

207 _nots = 0 # nlat * nlon

208 _pgm = None # PGM attributes, C{_PGM} or C{None}

209 _sizeB = 0 # geoid file size in bytes

210 _smooth = 0 # used only for RectBivariateSpline

211 _stdev = None # fixed in GeoidKarney

212 _u2B = 0 # np.itemsize or undefined

213 _yx_hits = None # cache hits, ala Karney's

214

215# _lat_d = _0_0 # increment, +tive

216# _lat_lo = _0_0 # lower lat, south

217# _lat_hi = _0_0 # upper lat, noth

218# _lon_d = _0_0 # increment, +tive

219# _lon_lo = _0_0 # left lon, west

220# _lon_hi = _0_0 # right lon, east

221# _lon_of = _0_0 # forward lon offset

222# _lon_og = _0_0 # reverse lon offset

223

224 def __init__(self, hs, p):

225 '''(INTERNAL) Set up the grid axes, the C{SciPy} interpolator and

226 several internal geoid attributes.

227

228 @arg hs: Grid knots with known height (C{numpy 2darray}).

229 @arg p: The C{slat, wlon, nlat, nlon, dlat, dlon} and other

230 geoid parameters (C{INTERNAL}).

231 '''

232 spi = self.scipy_interpolate

233 # for 2d scipy.interpolate.interp2d(xs, ys, hs, ...) and

234 # scipy.interpolate.RectBivariateSpline(ys, xs, hs, ...)

235 # require the shape of hs to be (len(ys), len(xs)), note

236 # the different (xs, ys, ...) and (ys, xs, ...) orders

237 if (p.nlat, p.nlon) != hs.shape:

238 raise GeoidError(shape=hs.shape, txt=_incompatible((p.nlat, p.nlon)))

239

240 # both axes and bounding box

241 ys, self._lat_d = self._gaxis2(p.slat, p.dlat, p.nlat, _lat_ + _s_)

242 xs, self._lon_d = self._gaxis2(p.wlon, p.dlon, p.nlon, _lon_ + _s_)

243

244 bb = ys[0], ys[-1], xs[0], xs[-1] + p.dlon # fudge lon_hi

245 # geoid grids are typically stored in row-major order, some

246 # with rows (90..-90) reversed and columns (0..360) wrapped

247 # to Easten longitude, 0 <= east < 180 and 180 <= west < 360

248 k = self.kind

249 if k in self._k2interp2d: # see _HeightBase

250 self._interp2d(xs, ys, hs, kind=k)

251 else: # XXX order ys and xs, see HeightLSQBiSpline

252 k = self._kxky(k)

253 self._ev = spi.RectBivariateSpline(ys, xs, hs, bbox=bb, ky=k, kx=k,

254 s=self._smooth).ev

255 self._hs_y_x = hs # numpy 2darray, row-major

256 self._nBytes = hs.nbytes # numpy size in bytes

257 self._nots = p.nots # grid nots len(hs)

258 self._fudge = p.dlon # see bb above

259 self._lon_of = float(p.flon) # forward offset

260 self._lon_og = g = float(p.glon) # reverse offset

261 # shrink the bounding box by 1 unit on every side:

262 # +self._lat_d, -self._lat_d, +self._lon_d, -self._lon_d

263 self._lat_lo, \

264 self._lat_hi, \

265 self._lon_lo, \

266 self._lon_hi = map(float, bb)

267 self._lon_lo -= g

268 self._lon_hi -= g

269

270 def __call__(self, *llis, **wrap_H):

271 '''Interpolate the geoid (or orthometric) height for one or more locations.

272

273 @arg llis: One or several locations (each C{LatLon}), all positional.

274 @kwarg wrap_H: Keyword arguments C{B{wrap}=False} (C{bool}) and

275 C{B{H}=False} (C{bool}). Use C{B{wrap}=True} to wrap

276 or I{normalize} all B{C{llis}} locations. If C{B{H}

277 is True}, return the I{orthometric} height instead of

278 the I{geoid} height at each location.

279

280 @return: A single geoid (or orthometric) height (C{float}) or

281 a list or tuple of geoid (or orthometric) heights (each

282 C{float}).

283

284 @raise GeoidError: Insufficient number of B{C{llis}}, an invalid

285 B{C{lli}} or the C{egm*.pgm} geoid file is closed.

286

287 @raise RangeError: An B{C{lli}} is outside this geoid's lat- or

288 longitude range.

289

290 @raise SciPyError: A C{scipy} issue.

291

292 @raise SciPyWarning: A C{scipy} warning as exception.

293

294 @note: To obtain I{orthometric} heights, each B{C{llis}} location

295 must have an ellipsoid C{height} or C{h} attribute, otherwise

296 C{height=0} is used.

297

298 @see: Function L{pygeodesy.heightOrthometric}.

299 '''

300 return self._called(llis, **wrap_H)

301

302 def __enter__(self):

303 '''Open context.

304 '''

305 return self

306

307 def __exit__(self, *unused): # PYCHOK exc_type, exc_value, exc_traceback)

308 '''Close context.

309 '''

310 self.close()

311 # return None # XXX False

312

313 def __repr__(self):

314 return self.toStr()

315

316 def __str__(self):

317 return Fmt.PAREN(self.classname, repr(self.name))

318

319 def _called(self, llis, wrap=False, H=False):

320 # handle __call__

321 _H = self._heightOrthometric if H else None

322 _as, llis = _as_llis2(llis, Error=GeoidError)

323 _w, hs = _Wrap._latlonop(wrap), []

324 _h, _a = self._hGeoid, hs.append

325 try:

326 for i, lli in enumerate(llis):

327 N = _h(*_w(lli.lat, lli.lon))

328 # orthometric or geoid height

329 _a(_H(lli, N) if _H else N)

330 return _as(hs)

331 except (GeoidError, RangeError) as x:

332 # XXX avoid str(LatLon()) degree symbols

333 n = _lli_ if _as is _ascalar else Fmt.INDEX(llis=i)

334 t = fstr((lli.lat, lli.lon), strepr=repr)

335 E = type(x)

336 raise E(n, t, wrap=wrap, H=H, cause=x)

337 except Exception as x:

338 if self._iscipy and self.scipy:

339 raise _SciPyIssue(x, self._ev_name)

340 else:

341 raise

342

343 @Property_RO

344 def _center(self):

345 '''(INTERNAL) Cache for method L{center}.

346 '''

347 return self._llh3(favg(self._lat_lo, self._lat_hi),

348 favg(self._lon_lo, self._lon_hi))

349

350 def center(self, LatLon=None):

351 '''Return the center location and height of this geoid.

352

353 @kwarg LatLon: Optional class to return the location and height

354 (C{LatLon}) or C{None}.

355

356 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,

357 height)} otherwise a B{C{LatLon}} instance with the lat-,

358 longitude and geoid height of the center grid location.

359 '''

360 return self._llh3LL(self._center, LatLon)

361

362 def close(self):

363 '''Close the C{egm*.pgm} geoid file if open (and applicable).

364 '''

365 if not self.closed:

366 self._egm.close()

367 self._egm = None

368

369 @property_RO

370 def closed(self):

371 '''Get the C{egm*.pgm} geoid file status.

372 '''

373 return self._egm is None

374

375 @property_RO

376 def cropped(self):

377 '''Is geoid cropped (C{bool} or C{None} if crop not supported).

378 '''

379 return self._cropped

380

381 @property_RO

382 def dtype(self):

383 '''Get the grid C{scipy} U{dtype<https://docs.SciPy.org/doc/numpy/

384 reference/generated/numpy.ndarray.dtype.html>} (C{numpy.dtype}).

385 '''

386 return self._hs_y_x.dtype

387

388 @property_RO

389 def endian(self):

390 '''Get the geoid endianess and U{dtype<https://docs.SciPy.org/

391 doc/numpy/reference/generated/numpy.dtype.html>} (C{str}).

392 '''

393 return self._endian

394

395 def _ev(self, y, x): # PYCHOK overwritten with .RectBivariateSpline.ev

396 # see methods _HeightBase._ev and -._interp2d

397 return self._ev2d(x, y) # (y, x) flipped!

398

399 def _gaxis2(self, lo, d, n, name):

400 # build grid axis, hi = lo + (n - 1) * d

401 m, a = len2(frange(lo, n, d))

402 if m != n:

403 raise LenError(type(self), grid=m, **{name: n})

404 if d < 0:

405 d, a = -d, list(reversed(a))

406 a = self.numpy.array(a)

407 m, i = min2(*map(float, a[1:] - a[:-1]))

408 if m < EPS: # non-increasing axis

409 i = Fmt.INDEX(name, i + 1)

410 raise GeoidError(i, m, txt_not_='increasing')

411 return a, d

412

413 def _g2ll2(self, lat, lon): # PYCHOK no cover

414 '''(INTERNAL) I{Must be overloaded}.'''

415 self._notOverloaded(lat, lon)

416

417 def _gyx2g2(self, y, x):

418 # convert grid (y, x) indices to grid (lat, lon)

419 return ((self._lat_lo + self._lat_d * y),

420 (self._lon_lo + self._lon_of + self._lon_d * x))

421

422 def height(self, lats, lons, **wrap):

423 '''Interpolate the geoid height for one or several lat-/longitudes.

424

425 @arg lats: Latitude or latitudes (each C{degrees}).

426 @arg lons: Longitude or longitudes (each C{degrees}).

427 @kwarg wrap: Use C{B{wrap}=True} to wrap or I{normalize} all

428 B{C{lats}} and B{C{lons}}.

429

430 @return: A single geoid height (C{float}) or a list of geoid

431 heights (each C{float}).

432

433 @raise GeoidError: Insufficient or unequal number of B{C{lats}}

434 and B{C{lons}}.

435

436 @raise RangeError: A B{C{lat}} or B{C{lon}} is outside this geoid's

437 lat- or longitude range.

438

439 @raise SciPyError: A C{scipy} issue.

440

441 @raise SciPyWarning: A C{scipy} warning as exception.

442 '''

443 return _HeightBase.height(self, lats, lons, **wrap)

444

445 def height_(self, *latlons, **wrap):

446 '''Interpolate the geoid height for each M{(latlons[i], latlons[i+1])

447 pair for i in range(0, len(latlons), B{2})}.

448

449 @arg latlons: Alternating lat-/longitude pairs (each C{degrees}),

450 all positional.

451

452 @see: Method L{height} for further details.

453

454 @return: A tuple of geoid heights (each C{float}).

455 '''

456 lls = tuple(self._as_lls(latlons[0::2], latlons[1::2]))

457 return self._called(lls, **wrap)

458

459 @property_ROver

460 def _heightOrthometric(self):

461 return _MODS.formy.heightOrthometric # overwrite property_ROver

462

463 def _hGeoid(self, lat, lon): # like GeoidQuasi._Nterpolate

464 out = self.outside(lat, lon)

465 if out: # XXX avoid str(LatLon()) degree symbols

466 t = fstr((lat, lon), strepr=repr)

467 raise RangeError(lli=t, txt=_SPACE_(_outside_, _on_, out))

468 return _float0d(self._ev(*self._ll2g2(lat, lon))) # scipy 1.18.0

469

470 @Property_RO

471 def _highest(self):

472 '''(INTERNAL) Cache for C{.highest}.

473 '''

474 return self._LL3T(self._llh3minmax(True), name__=self.highest)

475

476 def highest(self, LatLon=None, **unused):

477 '''Return the location and largest height of this geoid.

478

479 @kwarg LatLon: Optional class to return the location and height

480 (C{LatLon}) or C{None}.

481

482 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,

483 height)} otherwise a B{C{LatLon}} instance with the lat-,

484 longitude and geoid height of the highest grid location.

485 '''

486 return self._llh3LL(self._highest, LatLon)

487

488 @property_RO

489 def hits(self):

490 '''Get the number of cache hits (C{int} or C{None}).

491 '''

492 return self._yx_hits

493

494 @property_RO

495 def kind(self):

496 '''Get the interpolator kind and order (C{int}).

497 '''

498 return self._kind

499

500 def _kind_smooth(self, kind, smooth, name): # **name

501 # set C{kind}, C{smooth} and C{name}

502 if kind != 3:

503 self._kind = Int_(kind=kind, Error=GeoidError, low=-5, high=5)

504 if smooth:

505 self._smooth = Float_(smooth=smooth, Error=GeoidError, low=0)

506 if name:

507 _HeightBase.name.fset(self, _name__(**name)) # rename

508

509 @deprecated_property_RO

510 def knots(self):

511 '''DEPRECATED on 2026.08.26, use property C{nots}.'''

512 return self.nots

513

514 @property_RO

515 def nots(self):

516 '''Get the number of grid knots (C{int}).

517 '''

518 return self._nots

519

520 def _ll2g2(self, lat, lon): # PYCHOK no cover

521 '''(INTERNAL) I{Must be overloaded}.'''

522 self._notOverloaded(lat, lon)

523

524 @property_ROver

525 def _LL3T(self):

526 '''(INTERNAL) Get L{LatLon3Tuple}, I{once}.

527 '''

528 return _MODS.namedTuples.LatLon3Tuple # overwrite property_ROver

529

530 def _llh3(self, lat, lon):

531 return self._LL3T(lat, lon, self._hGeoid(lat, lon), name=self.name)

532

533 def _llh3LL(self, llh, LatLon):

534 return llh if LatLon is None else self._xnamed(LatLon(*llh))

535

536 def _llh3minmax(self, highest, *unused):

537 hs, np = self._hs_y_x, self.numpy

538 # <https://docs.SciPy.org/doc/numpy/reference/generated/

539 # numpy.argmin.html#numpy.argmin>

540 arg = np.argmax if highest else np.argmin

541 y, x = np.unravel_index(arg(hs, axis=None), hs.shape)

542 return self._g2ll2(*self._gyx2g2(y, x)) + (float(hs[y, x]),)

543

544 def _load(self, g, dtype=float, n=-1, offset=0, **sep): # sep=NN

545 # numpy.fromfile, like .frombuffer

546 g.seek(offset, _os.SEEK_SET)

547 return self.numpy.fromfile(g, dtype, count=n, **sep)

548

549 @Property_RO

550 def _lowerleft(self):

551 '''(INTERNAL) Cache for C{.lowerleft}.

552 '''

553 return self._llh3(self._lat_lo, self._lon_lo)

554

555 def lowerleft(self, LatLon=None):

556 '''Return the lower-left location and height of this geoid.

557

558 @kwarg LatLon: Optional class to return the location

559 (C{LatLon}) and height or C{None}.

560

561 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon, height)}

562 otherwise a B{C{LatLon}} instance with the lat-, longitude and

563 geoid height of the lower-left, SW grid corner.

564 '''

565 return self._llh3LL(self._lowerleft, LatLon)

566

567 @Property_RO

568 def _loweright(self):

569 '''(INTERNAL) Cache for C{.loweright}.

570 '''

571 return self._llh3(self._lat_lo, self._lon_hi - self._fudge)

572

573 def loweright(self, LatLon=None):

574 '''Return the lower-right location and height of this geoid.

575

576 @kwarg LatLon: Optional class to return the location and height

577 (C{LatLon}) or C{None}.

578

579 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon, height)}

580 otherwise a B{C{LatLon}} instance with the lat-, longitude and

581 geoid height of the lower-right, SE grid corner.

582 '''

583 return self._llh3LL(self._loweright, LatLon)

584

585 lowerright = loweright # synonymous

586

587 @Property_RO

588 def _lowest(self):

589 '''(INTERNAL) Cache for C{.lowest}.

590 '''

591 return self._LL3T(self._llh3minmax(False), name__=self.lowest)

592

593 def lowest(self, LatLon=None, **unused):

594 '''Return the location and lowest height of this geoid.

595

596 @kwarg LatLon: Optional class to return the location and height

597 (C{LatLon}) or C{None}.

598

599 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,

600 height)} otherwise a B{C{LatLon}} instance with the lat-,

601 longitude and geoid height of the lowest grid location.

602 '''

603 return self._llh3LL(self._lowest, LatLon)

604

605 @Property_RO

606 def mean(self):

607 '''Get the mean of this geoid's heights (C{float}).

608 '''

609 if self._mean is None: # see GeoidKarney

610 self._mean = float(self.numpy.mean(self._hs_y_x))

611 return self._mean

612

613 @property_RO

614 def name(self):

615 '''Get the name of this geoid (C{str}).

616 '''

617 return _HeightBase.name.fget(self) or self._geoid # recursion

618

619 @property_RO

620 def nBytes(self):

621 '''Get the grid in-memory size in bytes (C{int}).

622 '''

623 return self._nBytes

624

625 def _open(self, geoid, datum, kind, smooth, name):

626 # open the geoid file

627 if not _isin(datum, None, self._datum):

628 self._datum = _ellipsoidal_datum(datum, **name)

629 self._kind_smooth(kind, smooth, name)

630 try:

631 self._geoid = _os.path.basename(geoid)

632 self._sizeB = _os.path.getsize(geoid)

633 g = open(geoid, _rb_)

634 except (IOError, OSError) as x:

635 raise GeoidError(geoid=geoid, cause=x)

636 return g

637

638 def outside(self, lat, lon):

639 '''Check whether a location is outside this geoid's lat-/longitude

640 or crop range.

641

642 @arg lat: The latitude (C{degrees}).

643 @arg lon: The longitude (C{degrees}).

644

645 @return: A 1- or 2-character C{str} if outside, an empty C{str} otherwise.

646 '''

647 lat = _S_ if lat < self._lat_lo else (_N_ if lat > self._lat_hi else NN)

648 lon = _W_ if lon < self._lon_lo else (_E_ if lon > self._lon_hi else NN)

649 return NN(lat, lon) if lat and lon else (lat or lon)

650

651 @property_RO

652 def pgm(self):

653 '''Get the PGM attributes (C{_PGM} or C{None} if not available/applicable).

654 '''

655 return self._pgm

656

657 @property_RO

658 def shape(self):

659 '''Get the grid C{scipy} U{shape<https://docs.SciPy.org/doc/numpy/

660 reference/generated/numpy.ndarray.shape.html>} (C{tuple}).

661 '''

662 return tuple(self._hs_y_x.shape)

663

664 @property_RO

665 def sizeB(self):

666 '''Get the geoid grid file size in bytes (C{int}).

667 '''

668 return self._sizeB

669

670 @property_RO

671 def smooth(self):

672 '''Get the C{RectBivariateSpline} smoothing (C{float}).

673 '''

674 return self._smooth

675

676 @Property_RO

677 def stdev(self):

678 '''Get the standard deviation of this geoid's heights (C{float}) or C{None}.

679 '''

680 if self._stdev is None: # see GeoidKarney

681 self._stdev = float(self.numpy.std(self._hs_y_x))

682 return self._stdev

683

684 def _swne(self, crop):

685 # crop box to 4-tuple (s, w, n, e)

686 try:

687 if len(crop) == 2:

688 try: # sw, ne LatLons

689 swne = (crop[0].lat, crop[0].lon,

690 crop[1].lat, crop[1].lon)

691 except AttributeError: # (s, w), (n, e)

692 swne = tuple(crop[0]) + tuple(crop[1])

693 else: # (s, w, n, e)

694 swne = crop

695 if len(swne) == 4:

696 s, w, n, e = map(float, swne)

697 if _N_90_0 <= s <= (n - _1_0) <= 89.0 and \

698 _N_180_0 <= w <= (e - _1_0) <= 179.0:

699 return s, w, n, e

700 except (IndexError, TypeError, ValueError):

701 pass

702 raise GeoidError(crop=crop)

703

704 def toStr(self, prec=3, sep=_COMMASPACE_): # PYCHOK signature

705 '''This geoid and all geoid attributes as a string.

706

707 @kwarg prec: Number of decimal digits (0..9 or C{None} for

708 default). Trailing zero decimals are stripped

709 for B{C{prec}} values of 1 and above, but kept

710 for negative B{C{prec}} values.

711 @kwarg sep: Separator to join the attributes (C{str}).

712

713 @return: Geoid name and attributes (C{str}).

714 '''

715 t = attrs(self, _kind_, Nones=False) if self.kind < 0 else \

716 attrs(self, _kind_, 'smooth', Nones=False)

717 t += attrs(self, 'cropped', 'dtype', _endian_, 'hits', _mean_, 'nBytes', 'nots',

718 'shape', 'sizeB', _stdev_, prec=prec, Nones=False)

719 if self._iscipy: # all except GeoidKarney

720 for _n in (_numpy_, _scipy_):

721 try:

722 t += Fmt.EQUAL(_n, getattr(self, _n).version.version),

723 except ImportError:

724 pass

725 _n = _MODS.internals.typename

726 t += tuple(Fmt.EQUAL(_n(m), m().toStr(prec=prec)) for m in (self.center,

727 self.highest, self.lowest, self.lowerleft,

728 self.lowerright, self.upperleft, self.upperright))

729 return _COLONSPACE_(self, sep.join(t))

730

731 @property_RO

732 def u2B(self):

733 '''Get the PGM itemsize in bytes (C{int}).

734 '''

735 return self._u2B

736

737 @Property_RO

738 def _upperleft(self):

739 '''(INTERNAL) Cache for C{.upperleft}.

740 '''

741 return self._llh3(self._lat_hi, self._lon_lo)

742

743 def upperleft(self, LatLon=None):

744 '''Return the upper-left location and height of this geoid.

745

746 @kwarg LatLon: Optional class to return the location and height

747 (C{LatLon}) or C{None}.

748

749 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon, height)}

750 otherwise a B{C{LatLon}} instance with the lat-, longitude and

751 geoid height of the upper-left, NW grid corner.

752 '''

753 return self._llh3LL(self._upperleft, LatLon)

754

755 @Property_RO

756 def _upperright(self):

757 '''(INTERNAL) Cache for C{.upperright}.

758 '''

759 return self._llh3(self._lat_hi, self._lon_hi - self._fudge)

760

761 def upperright(self, LatLon=None):

762 '''Return the upper-right location and height of this geoid.

763

764 @kwarg LatLon: Optional class to return the location and height

765 (C{LatLon}) or C{None}.

766

767 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon, height)}

768 otherwise a B{C{LatLon}} instance with the lat-, longitude and

769 geoid height of the upper-right, NE grid corner.

770 '''

771 return self._llh3LL(self._upperright, LatLon)

772

773

774class GeoidEGM96(_GeoidBase):

775 '''Geoid height interpolator for the EGM96 U{15 Minute Interpolation Grid<https://earth-info.NGA.mil>}

776 based on C{SciPy} interpolation U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/reference/

777 generated/scipy.interpolate.RectBivariateSpline.html>}, U{interp2d<https://docs.SciPy.org/doc/scipy/

778 reference/generated/scipy.interpolate.interp2d.html>} or U{bisplrep/-ev<https://docs.scipy.org/doc/

779 scipy/reference/generated/scipy.interpolate.bisplrep.html>}.

780

781 Use only the C{WW15MGH.GRD} file, unzipped from the EGM96 U{15 Minute Interpolation Grid

782 <https://earth-info.NGA.mil/index.php?dir=wgs84&action=wgs84>} download.

783 '''

784 def __init__(self, EGM96_grd, datum=_WGS84, kind=3, smooth=0, **name_crop):

785 '''New L{GeoidEGM96} interpolator.

786

787 @arg EGM96_grd: An C{EGM96_grd} grid file name (C{.GRD}).

788 @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}),

789 overriding C{WGS84}.

790 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline

791 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.

792 RectBivariateSpline.html>} or -1, -3 or -5 for U{bisplrep/-ev<https://

793 docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}

794 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.

795 interpolate.interp2d.html>} C{linear}, C{cubic} respectively C{quintic},

796 see note for more details.

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

798 @kwarg name_crop: Optional geoid C{B{name}=NN} (C{str}) and UNSUPPORTED keyword argument

799 C{B{crop}=None}.

800

801 @raise GeoidError: Invalid B{C{crop}}, B{C{kind}} or B{C{smooth}} or a ECM96 grid file

802 B{C{ECM96_grd}} issue.

803

804 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.

805

806 @raise LenError: Grid file B{C{EGM96_grd}} axis mismatch.

807

808 @raise SciPyError: A C{scipy} issue.

809

810 @raise SciPyWarning: A C{scipy} warning as exception.

811

812 @raise TypeError: Invalid B{C{datum}}.

813

814 @note: Specify C{B{kind}=-1, -3 or -5} to use C{scipy.interpolate.interp2d}

815 before or C{scipy.interpolate.bisplrep/-ev} since C{Scipy} version 1.14.

816 '''

817 crop, name = _xkwds_pop2(name_crop, crop=None)

818 if crop is not None:

819 raise GeoidError(crop=crop, txt_not_=_supported_)

820

821 g = self._open(EGM96_grd, datum, kind, smooth, name)

822 _ = self.numpy # import numpy for .fromfile, .reshape

823

824 try:

825 p, hs = _Gpars(), self._load(g, sep=_SPACE_) # text

826 p.slat, n, p.wlon, e, p.dlat, p.dlon = hs[:6] # n-s, 0-E

827 p.nlat = int((n - p.slat) / p.dlat) + 1 # include S

828 p.nlon = int((e - p.wlon) / p.dlon) + 1 # include W

829 p.nots = p.nlat * p.nlon # inverted lats N downto S

830 p.glon = _180_0 # Eastern lons 0-360

831 hs = hs[6:].reshape(p.nlat, p.nlon)

832 _GeoidBase.__init__(self, hs, p)

833

834 except Exception as x:

835 raise _SciPyIssue(x, _in_, repr(EGM96_grd))

836 finally:

837 g.close()

838

839 def _g2ll2(self, lat, lon):

840 # convert grid (lat, lon) to earth (lat, lon)

841 return -lat, _lonE2lon(lon) # invert lat

842

843 def _ll2g2(self, lat, lon):

844 # convert earth (lat, lon) to grid (lat, lon)

845 return -lat, _lon2lonE(lon) # invert lat

846

847 if _FOR_DOCS:

848 __call__ = _GeoidBase.__call__

849 height = _GeoidBase.height

850 height_ = _GeoidBase.height_

851

852

853class GeoidG2012B(_GeoidBase):

854 '''Geoid height interpolator for U{GEOID12B Model

855 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/>} grids U{CONUS

856 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_CONUS.shtml>},

857 U{Alaska<https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_AK.shtml>},

858 U{Hawaii<https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_HI.shtml>},

859 U{Guam and Northern Mariana Islands

860 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_GMNI.shtml>},

861 U{Puerto Rico and U.S. Virgin Islands

862 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_PRVI.shtml>} and

863 U{American Samoa<https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_AS.shtml>}

864 based on C{SciPy} interpolation U{RectBivariateSpline<https://docs.SciPy.org/doc/

865 scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>}, U{interp2d

866 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}

867 or U{bisplrep/-ev<https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.

868 bisplrep.html>}.

869 '''

870 _datum = Datums.NAD83

871

872 def __init__(self, g2012b_bin, datum=Datums.NAD83, kind=3, smooth=0, **name_crop):

873 '''New L{GeoidG2012B} interpolator.

874

875 @arg g2012b_bin: A C{GEOID12B} grid file name (C{.bin}, see B{note}).

876 @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

877 L{a_f2Tuple}), overriding C{NAD83}.

878 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline

879 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.

880 RectBivariateSpline.html>} or -1, -3 or -5 for U{bisplrep/-ev<https://

881 docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}

882 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.

883 interpolate.interp2d.html>} C{linear}, C{cubic} respectively C{quintic},

884 see note for more details.

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

886 @kwarg name_crop: Optional geoid C{B{name}=NN} (C{str}) and UNSUPPORTED keyword argument

887 C{B{crop}=None}.

888

889 @raise GeoidError: Invalid B{C{crop}}, B{C{kind}} or B{C{smooth}} or a B{C{g2012b_bin}}

890 grid file issue.

891

892 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.

893

894 @raise LenError: Grid file B{C{g2012b_bin}} axis mismatch.

895

896 @raise SciPyError: A C{scipy} issue.

897

898 @raise SciPyWarning: A C{scipy} warning as exception.

899

900 @raise TypeError: Invalid B{C{datum}}.

901

902 @note: Only use any of the C{le} (little-endian) or C{be} (big-endian) C{g2012b*.bin}

903 I{binary} grid files.

904

905 @note: Specify C{B{kind}=-1, -3 or -5} to use C{scipy.interpolate.interp2d} I{before}

906 or C{scipy.interpolate.bisplrep/-ev} I{since} C{Scipy} version 1.14.

907 '''

908 crop, name = _xkwds_pop2(name_crop, crop=None)

909 if crop is not None:

910 raise GeoidError(crop=crop, txt_not_=_supported_)

911

912 g = self._open(g2012b_bin, datum, kind, smooth, name)

913 _ = self.numpy # import numpy for ._load and

914

915 try:

916 p = _Gpars()

917 n = (self.sizeB // 4) - 11 # number of f4 heights

918 # U{numpy dtype formats are different from Python struct formats

919 # <https://docs.SciPy.org/doc/numpy-1.15.0/reference/arrays.dtypes.html>}

920 for en_ in ('<', '>'):

921 # skip 4xf8, get 3xi4

922 p.nlat, p.nlon, ien = map(int, self._load(g, en_+'i4', 3, 32))

923 if ien == 1: # correct endian

924 p.nots = p.nlat * p.nlon

925 if p.nots == n and 1 < p.nlat < n \

926 and 1 < p.nlon < n:

927 self._endian = en_+'f4'

928 break

929 else: # couldn't validate endian

930 raise GeoidError(_endian_)

931

932 # get the first 4xf8

933 p.slat, p.wlon, p.dlat, p.dlon = map(float, self._load(g, en_+'f8', 4))

934 # read all f4 heights, ignoring the first 4xf8 and 3xi4

935 hs = self._load(g, self._endian, n, 44).reshape(p.nlat, p.nlon)

936 p.wlon -= _360_0 # western-most lon XXX _lonE2lon

937 _GeoidBase.__init__(self, hs, p)

938

939 except Exception as x:

940 raise _SciPyIssue(x, _in_, repr(g2012b_bin))

941 finally:

942 g.close()

943

944 def _g2ll2(self, lat, lon):

945 # convert grid (lat, lon) to earth (lat, lon)

946 return lat, lon

947

948 def _ll2g2(self, lat, lon):

949 # convert earth (lat, lon) to grid (lat, lon)

950 return lat, lon

951

952 if _FOR_DOCS:

953 __call__ = _GeoidBase.__call__

954 height = _GeoidBase.height

955 height_ = _GeoidBase.height_

956

957

958class GeoidHeight5Tuple(_NamedTuple):

959 '''5-Tuple C{(lat, lon, egm84, egm96, egm2008)} for U{GeoidHeights.dat

960 <https://SourceForge.net/projects/geographiclib/files/testdata/>}

961 tests with the heights for 3 different EGM grids at C{degrees90}

962 and C{degrees180} degrees (after converting C{lon} from original

963 C{0 <= EasterLon <= 360}).

964 '''

965 _Names_ = (_lat_, _lon_, 'egm84', 'egm96', 'egm2008')

966 _Units_ = ( Lat, Lon, Height, Height, Height)

967

968

969def _I(i):

970 '''(INTERNAL) Cache a single C{int} constant.

971 '''

972 i = int(i)

973 return _intCs.setdefault(i, i) # noqa: F821 del

974

975

976def _T(cs):

977 '''(INTERNAL) Cache a tuple of single C{int} constants.

978 '''

979 return tuple(map(_I, _splituple(cs)))

980

981_T0s12 = (_I(0),) * 12 # PYCHOK _T('0, 0, ..., 0')

982

983

984class GeoidKarney(_GeoidBase):

985 '''Geoid height interpolator for I{Karney}'s U{GeographicLib Earth

986 Gravitational Model (EGM)<https://GeographicLib.SourceForge.io/C++/doc/

987 geoid.html>} geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/

988 C++/doc/geoid.html#geoidinst>} datasets using bilinear or U{cubic

989 <https://dl.ACM.org/citation.cfm?id=368443>} interpolation and U{caching

990 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidcache>}

991 in pure Python, transcoded from I{Karney}'s U{C++ class Geoid

992 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinterp>}.

993

994 Use any of the geoid U{egm84-, egm96- or egm2008-*.pgm

995 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>}

996 datasets.

997 '''

998 _C0 = _F(372), _F(240), _F(372) # n, _ and s common denominators

999 # matrices c3n_, c3, c3s_, transposed from GeographicLib/Geoid.cpp

1000 _C3 = ((_T('0, 0, 62, 124, 124, 62, 0, 0, 0, 0, 0, 0'),

1001 _T0s12,

1002 _T('-131, 7, -31, -62, -62, -31, 45, 216, 156, -45, -55, -7'),

1003 _T0s12,

1004 _T('138, -138, 0, 0, 0, 0, -183, 33, 153, -3, 48, -48'), # PYCHOK indent

1005 _T('144, 42, -62, -124, -124, -62, -9, 87, 99, 9, 42, -42'),

1006 _T0s12,

1007 _T('0, 0, 0, 0, 0, 0, 93, -93, -93, 93, 0, 0'),

1008 _T('-102, 102, 0, 0, 0, 0, 18, 12, -12, -18, -84, 84'),

1009 _T('-31, -31, 31, 62, 62, 31, 0, -93, -93, 0, 31, 31')), # PYCHOK indent

1010

1011 (_T('9, -9, 9, 186, 54, -9, -9, 54, -54, 9, -9, 9'),

1012 _T('-18, 18, -88, -42, 162, -32, 8, -78, 78, -8, 18, -18'),

1013 _T('-88, 8, -18, -42, -78, 18, 18, 162, 78, -18, -32, -8'),

1014 _T('0, 0, 90, -150, 30, 30, 30, -90, 90, -30, 0, 0'),

1015 _T('96, -96, 96, -96, -24, 24, -96, -24, 144, -24, 24, -24'), # PYCHOK indent

1016 _T('90, 30, 0, -150, -90, 0, 0, 30, 90, 0, 30, -30'),

1017 _T('0, 0, -20, 60, -60, 20, -20, 60, -60, 20, 0, 0'),

1018 _T('0, 0, -60, 60, 60, -60, 60, -60, -60, 60, 0, 0'),

1019 _T('-60, 60, 0, 60, -60, 0, 0, 60, -60, 0, -60, 60'),

1020 _T('-20, -20, 0, 60, 60, 0, 0, -60, -60, 0, 20, 20')),

1021

1022 (_T('18, -18, 36, 210, 162, -36, 0, 0, 0, 0, -18, 18'), # PYCHOK indent

1023 _T('-36, 36, -165, 45, 141, -21, 0, 0, 0, 0, 36, -36'),

1024 _T('-122, -2, -27, -111, -75, 27, 62, 124, 124, 62, -64, 2'),

1025 _T('0, 0, 93, -93, -93, 93, 0, 0, 0, 0, 0, 0'),

1026 _T('120, -120, 147, -57, -129, 39, 0, 0, 0, 0, 66, -66'), # PYCHOK indent

1027 _T('135, 51, -9, -192, -180, 9, 31, 62, 62, 31, 51, -51'),

1028 _T0s12,

1029 _T('0, 0, -93, 93, 93, -93, 0, 0, 0, 0, 0, 0'),

1030 _T('-84, 84, 18, 12, -12, -18, 0, 0, 0, 0, -102, 102'),

1031 _T('-31, -31, 0, 93, 93, 0, -31, -62, -62, -31, 31, 31')))

1032

1033 _BT = (_T('0, 0'), # bilinear 4-tuple [i, j] indices

1034 _T('1, 0'),

1035 _T('0, 1'),

1036 _T('1, 1'))

1037

1038 _CM = (_T(' 0, -1'), # 10x12 cubic matrix [i, j] indices

1039 _T(' 1, -1'),

1040 _T('-1, 0'),

1041 _T(' 0, 0'), # _BT[0]

1042 _T(' 1, 0'), # _BT[1]

1043 _T(' 2, 0'),

1044 _T('-1, 1'),

1045 _T(' 0, 1'), # _BT[2]

1046 _T(' 1, 1'), # _BT[3]

1047 _T(' 2, 1'),

1048 _T(' 0, 2'),

1049 _T(' 1, 2'))

1050

1051# _cropped = None

1052 _endian = '>H' # struct.unpack 1 ushort (big endian, unsigned short)

1053 _4endian = '>4H' # struct.unpack 4 ushorts

1054 _Rendian = NN # struct.unpack a row of ushorts

1055# _highest = (-8.4, 147.367, 85.839) if egm2008-1.pgm else (

1056# (-8.167, 147.25, 85.422) if egm96-5.pgm else

1057# (-4.5, 148.75, 81.33)) # egm84-15.pgm

1058 _iscipy = False

1059# _lowest = (4.7, 78.767, -106.911) if egm2008-1.pgm else (

1060# (4.667, 78.833, -107.043) if egm96-5.pgm else

1061# (4.75, 79.25, -107.34)) # egm84-15.pgm

1062 _mean = _F(-1.317) # from egm2008-1, -1.438 egm96-5, -0.855 egm84-15

1063 _nBytes = None # not applicable

1064 _nterms = len(_C3[0]) # columns length, number of rows

1065 _smooth = None # not applicable

1066 _stdev = _F(29.244) # from egm2008-1, 29.227 egm96-5, 29.183 egm84-15

1067 _u2B = _calcsize(_endian) # pixelsize_ in bytes

1068 _4u2B = _calcsize(_4endian) # 4 pixelsize_s in bytes

1069 _Ru2B = 0 # row of pixelsize_s in bytes

1070 _yx_hits = 0 # cache hits

1071 _yx_i = () # cached (y, x) indices

1072 _yx_t = () # cached 4- or 10-tuple for _ev2k resp. _ev3k

1073

1074 def __init__(self, egm_pgm, crop=None, datum=_WGS84, kind=3, **name_smooth):

1075 '''New L{GeoidKarney} interpolator.

1076

1077 @arg egm_pgm: An U{EGM geoid dataset<https://GeographicLib.SourceForge.io/

1078 C++/doc/geoid.html#geoidinst>} file name (C{egm*.pgm}), see

1079 note below.

1080 @kwarg crop: Optional box to limit geoid locations, a 4-tuple (C{south,

1081 west, north, east}), 2-tuple (C{(south, west), (north, east)})

1082 with 2 C{degrees90} lat- and C{degrees180} longitudes or as

1083 2-tuple (C{LatLonSW, LatLonNE}) of C{LatLon} instances.

1084 @kwarg datum: Optional grid datum (C{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

1085 L{a_f2Tuple}), overriding C{WGS84}.

1086 @kwarg kind: Interpolation order (C{int}), 2 for C{bilinear} or 3 for C{cubic}.

1087 @kwarg name_smooth: Optional geoid C{B{name}=NN} (C{str}) and UNSUPPORTED

1088 keyword argument C{B{smooth}}, use C{B{smooth}=None} to ignore.

1089

1090 @raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid B{C{crop}},

1091 B{C{kind}} or B{C{smooth}}.

1092

1093 @raise TypeError: Invalid B{C{datum}}.

1094

1095 @see: Class L{GeoidPGM} and function L{egmGeoidHeights}.

1096

1097 @note: Geoid file B{C{egm_pgm}} remains open and I{must be closed} by calling

1098 method C{close} or by using C{with B{GeoidKarney}(...) as ...:} context.

1099 '''

1100 smooth, name = _xkwds_pop2(name_smooth, smooth=None)

1101 if smooth is not None:

1102 raise GeoidError(smooth=smooth, txt_not_=_supported_)

1103

1104 if _isin(kind, 2):

1105 self._ev2d = self._ev2k # see ._ev_name

1106 elif not _isin(kind, 3):

1107 raise GeoidError(kind=kind)

1108

1109 self._egm = g = self._open(egm_pgm, datum, kind, None, name)

1110 self._pgm = p = _PGM(g, pgm=egm_pgm, itemsize=self.u2B, sizeB=self.sizeB)

1111

1112 self._Rendian = self._4endian.replace(_4_, str(p.nlon))

1113 self._Ru2B = _calcsize(self._Rendian)

1114

1115 self._lon_of = float(p.flon) # forward offset

1116 self._lon_og = float(p.glon) # reverse offset

1117 # set earth (lat, lon) limits (s, w, n, e)

1118 self._lat_lo, self._lon_lo, \

1119 self._lat_hi, self._lon_hi = self._swne(crop if crop else p.crop4)

1120 self._cropped = bool(crop)

1121 self._nots = p.nots # number of grid knots

1122

1123 def _c0c3v(self, y, x):

1124 # get the common denominator, the 10x12 cubic matrix and

1125 # the 12 cubic v-coefficients around geoid index (y, x)

1126 p = self._pgm

1127 if 0 < x < (p.nlon - 2) and 0 < y < (p.nlat - 2):

1128 # read 4x4 ushorts, drop the 4 corners

1129 S = _os.SEEK_SET

1130 e = self._4endian

1131 g = self._egm

1132 n = self._4u2B

1133 R = self._Ru2B

1134 b = self._seek(y - 1, x - 1)

1135 v = _unpack(e, g.read(n))[1:3]

1136 b += R

1137 g.seek(b, S)

1138 v += _unpack(e, g.read(n))

1139 b += R

1140 g.seek(b, S)

1141 v += _unpack(e, g.read(n))

1142 b += R

1143 g.seek(b, S)

1144 v += _unpack(e, g.read(n))[1:3]

1145 j = 1

1146

1147 else: # likely some wrapped y and/or x's

1148 v = self._raws(y, x, GeoidKarney._CM)

1149 j = 0 if y < 1 else (1 if y < (p.nlat - 2) else 2)

1150

1151 return GeoidKarney._C0[j], GeoidKarney._C3[j], v

1152

1153 @property_RO

1154 def dtype(self):

1155 '''Get the geoid's grid data type (C{str}).

1156 '''

1157 return 'ushort'

1158

1159 def _ev(self, lat, lon): # PYCHOK expected

1160 # interpolate the geoid height at grid (lat, lon)

1161 fy, fx = self._g2yx2(lat, lon)

1162 y, x = int(_floor(fy)), int(_floor(fx))

1163 fy -= y

1164 fx -= x

1165 H = self._ev2d(fy, fx, y, x) # PYCHOK ._ev3k or ._ev2k

1166 H *= self._pgm.Scale # H.fmul(self._pgm.Scale)

1167 H += self._pgm.Offset # H.fadd(self._pgm.Offset)

1168 return H.fsum() # float(H)

1169

1170 def _ev2k(self, fy, fx, *yx):

1171 # compute the bilinear 4-tuple and interpolate raw H

1172 if self._yx_i == yx:

1173 self._yx_hits += 1

1174 else:

1175 y, x = self._yx_i = yx

1176 self._yx_t = self._raws(y, x, GeoidKarney._BT)

1177 t = self._yx_t

1178 v = _1_0, (-fx), fx

1179 H = _Dotf(v, t[0], t[0], t[1]).fmul(_1_0 - fy) # c = a * (1 - fy)

1180 H += _Dotf(v, t[2], t[2], t[3]).fmul(fy) # c += b * fy

1181 return H # Fsum

1182

1183 def _ev3k(self, fy, fx, *yx):

1184 # compute the cubic 10-tuple and interpolate raw H

1185 if self._yx_i == yx:

1186 self._yx_hits += 1

1187 else:

1188 c0, c3, v = self._c0c3v(*yx)

1189 # assert len(c3) == self._nterms

1190 self._yx_t = tuple(_Dotf(v, *r3).fover(c0) for r3 in c3)

1191 self._yx_i = yx

1192 # GeographicLib/Geoid.cpp Geoid::height(lat, lon) ...

1193 # real h = t[0] + fx * (t[1] + fx * (t[3] + fx * t[6])) +

1194 # fy * (t[2] + fx * (t[4] + fx * t[7]) +

1195 # fy * (t[5] + fx * t[8] + fy * t[9]));

1196 t = self._yx_t

1197 v = _1_0, fx, fy

1198 H = _Dotf(v, t[5], t[8], t[9])

1199 H *= fy

1200 H += _Hornerf(fx, t[2], t[4], t[7])

1201 H *= fy

1202 H += _Hornerf(fx, t[0], t[1], t[3], t[6])

1203 return H # Fsum

1204

1205 _ev2d = _ev3k # overriden for kind=2, see ._ev_name

1206

1207 def _g2ll2(self, lat, lon):

1208 # convert grid (lat, lon) to earth (lat, lon), uncropped

1209 return lat, _lonE2lon(lon)

1210

1211 def _g2yx2(self, lat, lon):

1212 # convert grid (lat, lon) to grid (y, x) indices

1213 p = self._pgm

1214 # note, slat = +90, rlat < 0 makes y >=0

1215 return ((lat - p.slat) * p.rlat), ((lon - p.wlon) * p.rlon)

1216

1217 def _gyx2g2(self, y, x):

1218 # convert grid (y, x) indices to grid (lat, lon)

1219 p = self._pgm

1220 return (p.slat + p.dlat * y), (p.wlon + p.dlon * x)

1221

1222 @Property_RO

1223 def _highest_ltd(self):

1224 '''(INTERNAL) Cache for C{.highest}.

1225 '''

1226 return self._LL3T(self._llh3minmax(True, -12, -4), name__=self.highest)

1227

1228 def highest(self, LatLon=None, full=False): # PYCHOK full

1229 '''Return the location and largest height of this geoid.

1230

1231 @kwarg LatLon: Optional class to return the location and height

1232 (C{LatLon}) or C{None}.

1233 @kwarg full: Search the full or limited latitude range (C{bool}).

1234

1235 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,

1236 height)} otherwise a B{C{LatLon}} instance with the lat-,

1237 longitude and geoid height of the highest grid location.

1238 '''

1239 llh = self._highest if full or self.cropped else self._highest_ltd

1240 return self._llh3LL(llh, LatLon)

1241

1242 def _lat2y2(self, lat2):

1243 # convert earth lat(s) to min and max grid y indices

1244 ys, m = [], self._pgm.nlat - 1

1245 for lat in lat2:

1246 y, _ = self._g2yx2(*self._ll2g2(lat, 0))

1247 ys.append(max(min(int(y), m), 0))

1248 return min(ys), max(ys) + 1

1249

1250 def _ll2g2(self, lat, lon):

1251 # convert earth (lat, lon) to grid (lat, lon), uncropped

1252 return lat, _lon2lonE(lon)

1253

1254 def _llh3minmax(self, highest, *lat2):

1255 # find highest or lowest, takes 10+ secs for egm2008-1.pgm geoid

1256 # (Python 2.7.16, macOS 10.13.6 High Sierra, iMac 3 GHz Core i3)

1257 if highest:

1258 def _mt(r, h):

1259 m = max(r)

1260 return m, (m > h)

1261

1262 else: # lowest

1263 def _mt(r, h): # PYCHOK redef

1264 m = min(r)

1265 return m, (m < h)

1266

1267 y = x = 0

1268 h = self._raw(y, x)

1269 for j, r in self._raw2(*lat2):

1270 m, t = _mt(r, h)

1271 if t:

1272 h, y, x = m, j, r.index(m)

1273 h *= self._pgm.Scale

1274 h += self._pgm.Offset

1275 return self._g2ll2(*self._gyx2g2(y, x)) + (h,)

1276

1277 @Property_RO

1278 def _lowest_ltd(self):

1279 '''(INTERNAL) Cache for C{.lowest}.

1280 '''

1281 return self._LL3T(self._llh3minmax(False, 0, 8), name__=self.lowest)

1282

1283 def lowest(self, LatLon=None, full=False): # PYCHOK full

1284 '''Return the location and lowest height of this geoid.

1285

1286 @kwarg LatLon: Optional class to return the location and height

1287 (C{LatLon}) or C{None}.

1288 @kwarg full: Search the full or limited latitude range (C{bool}).

1289

1290 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,

1291 height)} otherwise a B{C{LatLon}} instance with the lat-,

1292 longitude and geoid height of the lowest grid location.

1293 '''

1294 llh = self._lowest if full or self.cropped else self._lowest_ltd

1295 return self._llh3LL(llh, LatLon)

1296

1297 def _raw(self, y, x):

1298 # get the ushort geoid height at geoid index (y, x),

1299 # like GeographicLib/Geoid.hpp real rawval(is, iy)

1300 p = self._pgm

1301 if x < 0:

1302 x += p.nlon

1303 elif x >= p.nlon:

1304 x -= p.nlon

1305 h = p.nlon // 2

1306 if y < 0:

1307 y = -y

1308 elif y >= p.nlat:

1309 y = (p.nlat - 1) * 2 - y

1310 else:

1311 h = 0

1312 x += h if x < h else -h

1313 self._seek(y, x)

1314 h = _unpack(self._endian, self._egm.read(self._u2B))

1315 return h[0]

1316

1317 def _raws(self, y, x, ijs):

1318 # get bilinear 4-tuple or 10x12 cubic matrix

1319 return tuple(self._raw(y + j, x + i) for i, j in ijs)

1320

1321 def _raw2(self, *lat2):

1322 # yield a 2-tuple (y, ushorts) for each row or for

1323 # the rows between two (or more) earth lat values

1324 p = self._pgm

1325 g = self._egm

1326 e = self._Rendian

1327 n = self._Ru2B

1328 # min(lat2) <= lat <= max(lat2) or 0 <= y < p.nlat

1329 s, t = self._lat2y2(lat2) if lat2 else (0, p.nlat)

1330 self._seek(s, 0) # to start of row s

1331 for y in range(s, t):

1332 yield y, _unpack(e, g.read(n))

1333

1334 def _seek(self, y, x):

1335 # position geoid to grid index (y, x)

1336 p, g = self._pgm, self._egm

1337 if g:

1338 b = p.skip + (y * p.nlon + x) * self._u2B

1339 g.seek(b, _os.SEEK_SET)

1340 return b # position

1341 raise GeoidError('closed file', txt=repr(p.egm)) # IOError

1342

1343 @property_RO

1344 def shape(self):

1345 '''Get the geoid shape (C{tuple}).

1346 '''

1347 p = self._pgm

1348 return p.nlat, p.nlon

1349

1350

1351class GeoidPGM(_GeoidBase):

1352 '''Geoid height interpolator for I{Karney}'s U{GeographicLib Earth

1353 Gravitational Model (EGM)<https://GeographicLib.SourceForge.io/C++/doc/geoid.html>}

1354 geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>}

1355 datasets but based on C{SciPy} U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/

1356 reference/generated/scipy.interpolate.RectBivariateSpline.html>}, U{bisplrep/-ev

1357 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}

1358 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.

1359 interp2d.html>} interpolation.

1360

1361 Use any of the U{egm84-, egm96- or egm2008-*.pgm <https://GeographicLib.SourceForge.io/

1362 C++/doc/geoid.html#geoidinst>} datasets. However, unless cropped, an entire C{egm*.pgm}

1363 dataset is loaded into the C{SciPy} interpolator and converted from 2-byte C{int} to

1364 8-byte C{dtype float64}. Therefore, internal memory usage is 4x the U{egm*.pgm

1365 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>} file size and may

1366 exceed the available memory, especially with 32-bit Python, see properties C{.nBytes}

1367 and C{.sizeB}.

1368 '''

1369 _cropped = False

1370 _endian = '>u2'

1371

1372 def __init__(self, egm_pgm, crop=None, datum=_WGS84, kind=3, smooth=0, **name):

1373 '''New L{GeoidPGM} interpolator.

1374

1375 @arg egm_pgm: An U{EGM geoid dataset<https://GeographicLib.SourceForge.io/

1376 C++/doc/geoid.html#geoidinst>} file name (C{egm*.pgm}).

1377 @kwarg crop: Optional box to crop B{C{egm_pgm}}, a 4-tuple (C{south, west,

1378 north, east}) or 2-tuple (C{(south, west), (north, east)}),

1379 in C{degrees90} lat- and C{degrees180} longitudes or a 2-tuple

1380 (C{LatLonSW, LatLonNE}) of C{LatLon} instances.

1381 @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

1382 L{a_f2Tuple}), overriding C{WGS84}.

1383 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline

1384 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.

1385 RectBivariateSpline.html>} or -1, -3 or -5 for U{bisplrep/-ev<https://

1386 docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}

1387 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.

1388 interpolate.interp2d.html>} C{linear}, C{cubic} respectively C{quintic},

1389 see note for more details.

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

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

1392

1393 @raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid B{C{crop}}, B{C{kind}}

1394 or B{C{smooth}}.

1395

1396 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.

1397

1398 @raise LenError: EGM dataset B{C{egm_pgm}} axis mismatch.

1399

1400 @raise SciPyError: A C{scipy} issue.

1401

1402 @raise SciPyWarning: A C{scipy} warning as exception.

1403

1404 @raise TypeError: Invalid B{C{datum}} or unexpected argument.

1405

1406 @note: Specify C{B{kind}=-1, -3 or -5} to use C{scipy.interpolate.interp2d} before

1407 or C{scipy.interpolate.bisplrep/-ev} since C{Scipy} version 1.14.

1408

1409 @note: The U{GeographicLib egm*.pgm<https://GeographicLib.SourceForge.io/C++/doc/

1410 geoid.html#geoidinst>} file sizes are based on a 2-byte C{int} height

1411 converted to 8-byte C{dtype float64} for C{scipy} interpolators. Therefore,

1412 internal memory usage is 4 times the C{egm*.pgm} file size and may exceed

1413 the available memory, especially with 32-bit Python. To reduce memory

1414 usage, use keyword argument B{C{crop}} to the region of interest. For

1415 example C{B{crop}=(20, -125, 50, -65)} covers the U{conterminous US

1416 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/maps/GEOID12B_CONUS_grids.png>}

1417 (CONUS), less than 3% of the entire C{egm2008-1.pgm} dataset.

1418

1419 @see: Class L{GeoidKarney} and function L{egmGeoidHeights}.

1420 '''

1421 np = self.numpy

1422 self._u2B = np.dtype(self.endian).itemsize

1423

1424 g = self._open(egm_pgm, datum, kind, smooth, name)

1425 self._pgm = p = _PGM(g, pgm=egm_pgm, itemsize=self.u2B, sizeB=self.sizeB)

1426 if crop:

1427 g = p._cropped(g, abs(kind) + 1, *self._swne(crop))

1428 if _MODS.internals._version2(np.__version__) < (1, 9):

1429 g = open(g.name, _rb_) # reopen tempfile for numpy 1.8.0-

1430 self._cropped = True

1431 try:

1432 # U{numpy dtype formats are different from Python struct formats

1433 # <https://docs.SciPy.org/doc/numpy-1.15.0/reference/arrays.dtypes.html>}

1434 # read all heights, skipping the PGM header lines, converted to float

1435 hs = self._load(g, self.endian, p.nots, p.skip).reshape(p.nlat, p.nlon) * p.Scale

1436 if p.Offset: # offset

1437 hs = p.Offset + hs

1438 if p.dlat < 0: # flip the rows

1439 hs = np.flipud(hs)

1440 _GeoidBase.__init__(self, hs, p)

1441 except Exception as x:

1442 raise _SciPyIssue(x, _in_, repr(egm_pgm))

1443 finally:

1444 g.close()

1445

1446 def _g2ll2(self, lat, lon):

1447 # convert grid (lat, lon) to earth (lat, lon), un-/cropped

1448 if self._cropped:

1449 lon -= self._lon_of

1450 else:

1451 lon = _lonE2lon(lon)

1452 return lat, lon

1453

1454 def _ll2g2(self, lat, lon):

1455 # convert earth (lat, lon) to grid (lat, lon), un-/cropped

1456 if self._cropped:

1457 lon += self._lon_of

1458 else:

1459 lon = _lon2lonE(lon)

1460 return lat, lon

1461

1462 if _FOR_DOCS:

1463 __call__ = _GeoidBase.__call__

1464 height = _GeoidBase.height

1465 height_ = _GeoidBase.height_

1466

1467

1468class GeoidQuasi(_GeoidBase): # PYCHOK no cover

1469 '''Quasi-geoid height interpolator for C{1-degree, whole Earth grids}, used in package

1470 C{PyAxQG}, for example.

1471 '''

1472

1473 def __init__(self, knots, dtype=float, kind=3, smooth=0, name=NN, **pars):

1474 '''New L{GeoidQuasi} interpolator.

1475

1476 @arg knots: Geoid heights in row-major order (C{iterator} over C{nlat} * C{nlon}

1477 scalars from C{slat} north and C{wlon} east).

1478 @kwarg dtype: NumPy C{dtype} to use for the B{C{knots}} (C{str} or C{numpy.dtype}).

1479 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline

1480 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.

1481 RectBivariateSpline.html>} or -1, -3 or -5 for U{bisplrep/-ev<https://

1482 docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}

1483 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.

1484 interpolate.interp2d.html>} C{linear}, C{cubic} respectively C{quintic},

1485 see note for more details.

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

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

1488 @kwarg pars: Optional geoid parameters, overriding the defaults (C{slat=-90.0,

1489 wlon=-180.0, dlat=1.0, dlon=1.0, nlat=181 and nlon=361}).

1490

1491 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.

1492

1493 @raise SciPyError: A C{scipy} issue.

1494

1495 @raise SciPyWarning: A C{scipy} warning as exception.

1496

1497 @note: Specify C{B{kind}=-1, -3 or -5} to use C{scipy.interpolate.interp2d} before or

1498 C{scipy.interpolate.bisplrep/-ev} since C{Scipy} version 1.14.

1499 '''

1500 self._kind_smooth(kind, smooth, name)

1501 _ = self.scipy

1502 p = _Gpars(slat=-90, wlon=-180,

1503 dlat=1, dlon=1, dtype=dtype,

1504 nlat=181, nlon=361).update(**pars)

1505

1506 np = self.numpy

1507 ks = np.fromiter(knots, p.dtype) # len(knots) == p.nots

1508 ks = ks.reshape(p.nlat, p.nlon)

1509 _GeoidBase.__init__(self, ks, p)

1510

1511 def _g2ll2(self, lat, lon):

1512 # convert grid (lat, lon) to earth (lat, lon)

1513 return lat, lon

1514

1515 def _Nterpolate(self, lat, lon): # like ._hGeoid

1516 # return geoid height for safe C{lat} and C{lon} from pyaxqg

1517 try:

1518 return _float0d(self._ev(lat, lon)) # scipy 1.18.0

1519 except Exception as x:

1520 if self._iscipy and self.scipy: # True

1521 raise _SciPyIssue(x, self._ev_name)

1522 else:

1523 raise GeoidError((lat, lon), cause=x)

1524

1525 def _ll2g2(self, lat, lon):

1526 # convert earth (lat, lon) to grid (lat, lon)

1527 return lat, lon

1528

1529 if _FOR_DOCS:

1530 __call__ = _GeoidBase.__call__

1531 height = _GeoidBase.height

1532 height_ = _GeoidBase.height_

1533

1534

1535class _Gpars(_Named):

1536 '''(INTERNAL) Basic geoid parameters.

1537 '''

1538 # interpolator parameters

1539 dlat = 0 # +/- latitude resolution in C{degrees}

1540 dlon = 0 # longitude resolution in C{degrees}

1541 nlat = 1 # number of latitude knots (C{int})

1542 nlon = 0 # number of longitude knots (C{int})

1543 nots = 0 # number of knots, nlat * nlon (C{int})

1544 rlat = 0 # +/- latitude resolution in C{float}, 1 / .dlat

1545 rlon = 0 # longitude resolution in C{float}, 1 / .dlon

1546 slat = 0 # nothern- or southern most latitude (C{degrees90})

1547 wlon = 0 # western-most longitude in Eastern lon (C{degrees360})

1548

1549 flon = 0 # forward, earth to grid longitude offset

1550 glon = 0 # reverse, grid to earth longitude offset

1551

1552 dtype = float # numpy.dtype, 'f4', 'f8' == float

1553 skip = 0 # header bytes to skip (C{int})

1554

1555 def __init__(self, **pars):

1556 if pars:

1557 self.update(**pars)

1558

1559 def __repr__(self):

1560 t = _COMMASPACE_.join(pairs((a, getattr(self, a)) for

1561 a in dir(self.__class__)

1562 if a[:1].isupper()))

1563 return _COLONSPACE_(self, t)

1564

1565 def __str__(self):

1566 return Fmt.PAREN(self.classname, repr(self.name))

1567

1568 def update(self, **pars):

1569 if pars:

1570 self.__dict__.update(pars)

1571 self.nots = self.nlat * self.nlon

1572 return self

1573

1574

1575class _PGM(_Gpars):

1576 '''(INTERNAL) Parse an C{egm*.pgm} geoid dataset file.

1577

1578 # Geoid file in PGM format for the GeographicLib::Geoid class

1579 # Description WGS84 EGM96, 5-minute grid

1580 # URL https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/egm96/egm96.html

1581 # DateTime 2009-08-29 18:45:03

1582 # MaxBilinearError 0.140

1583 # RMSBilinearError 0.005

1584 # MaxCubicError 0.003

1585 # RMSCubicError 0.001

1586 # Offset -108

1587 # Scale 0.003

1588 # Origin 90N 0E

1589 # AREA_OR_POINT Point

1590 # Vertical_Datum WGS84

1591 <width> <height>

1592 <pixel>

1593 ...

1594 '''

1595 crop4 = () # 4-tuple (C{south, west, north, east}).

1596 egm = None

1597 glon = 180 # reverse offset, uncropped

1598# pgm = NN # name

1599 sizeB = 0

1600 u2B = 2 # item size of grid height (C{int}).

1601

1602 @staticmethod

1603 def _llstr2floats(latlon):

1604 # llstr to (lat, lon) floats

1605 lat, lon = latlon.split()

1606 return _MODS.dms.parseDMS2(lat, lon)

1607

1608 # PGM file attributes, CamelCase but not .istitle()

1609 AREA_OR_POINT = str

1610 DateTime = str

1611 Description = str # 'WGS84 EGM96, 5-minute grid'

1612 Geoid = str # 'file in PGM format for the GeographicLib::Geoid class'

1613 MaxBilinearError = float

1614 MaxCubicError = float

1615 Offset = float

1616 Origin = _llstr2floats

1617 Pixel = 0

1618 RMSBilinearError = float

1619 RMSCubicError = float

1620 Scale = float

1621 URL = str # 'https://Earth-Info.NGA.mil/GandG/wgs84/...'

1622 Vertical_Datum = str

1623

1624 def __init__(self, g, pgm=NN, itemsize=0, sizeB=0): # MCCABE 22

1625 '''(INTERNAL) New C{_PGM} parsed C{egm*.pgm} geoid dataset.

1626 '''

1627 self.name = pgm # geoid file name

1628 if itemsize:

1629 self._u2B = itemsize

1630 if sizeB:

1631 self.sizeB = sizeB

1632

1633 t = g.readline() # make sure newline == '\n'

1634 if t != b'P5\n' and t.strip() != b'P5':

1635 raise self._Errorf(_format_, _header_, t)

1636

1637 while True: # read all # Attr ... lines,

1638 try: # ignore empty ones or comments

1639 t = g.readline().strip()

1640 if t.startswith(_bHASH_):

1641 t = t.lstrip(_bHASH_).lstrip()

1642 a, v = map(_ub2str, t.split(None, 1))

1643 f = getattr(_PGM, a, None)

1644 if callable(f) and a[:1].isupper():

1645 setattr(self, a, f(v))

1646 elif t:

1647 break

1648 except (TypeError, ValueError):

1649 raise self._Errorf(_format_, 'Attr', t)

1650 else: # should never get here

1651 raise self._Errorf(_format_, _header_, g.tell())

1652

1653 try: # must be (even) width and (odd) height

1654 nlon, nlat = map(int, t.split())

1655 if nlon < 2 or nlon > (360 * 60) or isodd(nlon) or \

1656 nlat < 2 or nlat > (181 * 60) or not isodd(nlat):

1657 raise ValueError

1658 except (TypeError, ValueError):

1659 raise self._Errorf(_format_, _SPACE_(_width_, _height_), t)

1660

1661 try: # must be 16 bit pixel height

1662 t = g.readline().strip()

1663 self.Pixel = int(t)

1664 if not 255 < self.Pixel < 65536: # >u2 or >H only

1665 raise ValueError

1666 except (TypeError, ValueError):

1667 raise self._Errorf(_format_, 'pixel', t)

1668

1669 for a in dir(_PGM): # set undefined # Attr ... to None

1670 if a[:1].isupper() and callable(getattr(self, a)):

1671 setattr(self, a, None)

1672

1673 if self.Origin is None:

1674 raise self._Errorf(_format_, 'Origin', self.Origin)

1675 if self.Offset is None or self.Offset > 0:

1676 raise self._Errorf(_format_, 'Offset', self.Offset)

1677 if self.Scale is None or self.Scale < EPS:

1678 raise self._Errorf(_format_, 'Scale', self.Scale)

1679

1680 self.skip = g.tell()

1681 self.nots = nlat * nlon

1682

1683 self.nlat, self.nlon = nlat, nlon

1684 self.slat, self.wlon = self.Origin

1685 # note, negative .dlat and .rlat since rows

1686 # are from .slat 90N down in decreasing lat

1687 self.dlat, self.dlon = (_180_0 / (1 - nlat)), (_360_0 / nlon)

1688 self.rlat, self.rlon = ((1 - nlat) / _180_0), (nlon / _360_0)

1689

1690 # grid corners in earth (lat, lon), .slat = 90, .dlat < 0

1691 n = float(self.slat)

1692 s = n + self.dlat * (nlat - 1)

1693 w = self.wlon - self.glon

1694 e = w + self.dlon * nlon

1695 self.crop4 = s, w, n, e

1696

1697 n = self.sizeB - self.skip

1698 if n > 0 and n != (self.nots * self.u2B):

1699 raise self._Errorf('%s(%s x %s != %s)', _assert_, nlat, nlon, n)

1700

1701 def _cropped(self, g, k1, south, west, north, east): # MCCABE 15

1702 '''Crop the geoid to (south, west, north, east) box.

1703 '''

1704 # flon offset for both west and east

1705 f = 360 if west < 0 else 0

1706 # earth (lat, lon) to grid indices (y, x),

1707 # note y is decreasing, i.e. n < s

1708 s, w = self._lle2yx2(south, west, f)

1709 n, e = self._lle2yx2(north, east, f)

1710 s += 1 # s > n

1711 e += 1 # e > w

1712

1713 hi, wi = self.nlat, self.nlon

1714 # handle special cases

1715 if (s - n) > hi:

1716 n, s = 0, hi # entire lat range

1717 if (e - w) > wi:

1718 w, e, f = 0, wi, 180 # entire lon range

1719 if s == hi and w == n == 0 and e == wi:

1720 return g # use entire geoid as-is

1721

1722 if (e - w) < k1 or (s - n) < (k1 + 1):

1723 raise self._Errorf(_format_, 'swne', (north - south, east - west))

1724

1725 if e > wi > w: # wrap around

1726 # read w..wi and 0..e

1727 r, p = (wi - w), (e - wi)

1728 elif e > w:

1729 r, p = (e - w), 0

1730 else:

1731 raise self._Errorf('%s(%s < %s)', _assert_, w, e)

1732

1733 # convert to bytes

1734 r *= self.u2B

1735 p *= self.u2B

1736 q = wi * self.u2B # stride

1737 # number of rows and cols to skip from

1738 # the original (.slat, .wlon) origin

1739 z = self.skip + (n * wi + w) * self.u2B

1740 # sanity check

1741 if r < 2 or p < 0 or q < 2 or z < self.skip \

1742 or z > self.sizeB:

1743 raise self._Errorf(_format_, _assert_, (r, p, q, z))

1744

1745 # can't use _BytesIO since numpy

1746 # needs .fileno attr in .fromfile

1747 t, c = 0, self._tmpfile()

1748 # reading (s - n) rows, forward

1749 for y in range(n, s): # PYCHOK y unused

1750 g.seek(z, _os.SEEK_SET)

1751 # Python 2 tmpfile.write returns None

1752 t += c.write(g.read(r)) or r

1753 if p: # wrap around to start of row

1754 g.seek(-q, _os.SEEK_CUR)

1755 # assert(g.tell() == (z - w * self.u2B))

1756 # Python 2 tmpfile.write returns None

1757 t += c.write(g.read(p)) or p

1758 z += q

1759 c.flush()

1760 g.close()

1761

1762 s -= n # nlat

1763 e -= w # nlon

1764 k = s * e # nots

1765 z = k * self.u2B

1766 if t != z:

1767 raise self._Errorf('%s(%s != %s) %s', _assert_, t, z, self)

1768

1769 # update the _Gpars accordingly, note attributes

1770 # .dlat, .dlon, .rlat and .rlon remain unchanged

1771 self.slat += n * self.dlat

1772 self.wlon += w * self.dlon

1773 self.nlat = s

1774 self.nlon = e

1775 self.flon = self.glon = f

1776

1777 self.crop4 = south, west, north, east

1778 self.nots = k

1779 self.skip = 0 # no header lines in c

1780

1781 c.seek(0, _os.SEEK_SET)

1782 # c = open(c.name, _rb_) # reopen for numpy 1.8.0-

1783 return c

1784

1785 def _Errorf(self, fmt, *args): # PYCHOK no cover

1786 t = fmt % args

1787 e = self.pgm or NN

1788 if e:

1789 t = _SPACE_(t, _in_, repr(e))

1790 return PGMError(t)

1791

1792 def _lle2yx2(self, lat, lon, flon):

1793 # earth (lat, lon) to grid indices (y, x)

1794 # with .dlat decreasing from 90N .slat

1795 lat -= self.slat

1796 lon += flon - self.wlon

1797 return (min(self.nlat - 1, max(0, int(lat * self.rlat))),

1798 max(0, int(lon * self.rlon)))

1799

1800 def _tmpfile(self):

1801 # create a tmpfile to hold the cropped geoid grid

1802 try:

1803 from tempfile import NamedTemporaryFile as tmpfile

1804 except ImportError: # Python 2.7.16-

1805 from _os import tmpfile # PYCHOK from

1806 t = _os.path.basename(self.pgm)

1807 t = _os.path.splitext(t)[0]

1808 f = tmpfile(mode='w+b', prefix=t or 'egm')

1809 f.seek(0, _os.SEEK_SET) # force overwrite

1810 return f

1811

1812 @Property_RO

1813 def pgm(self):

1814 '''Get the geoid file name (C{str}).

1815 '''

1816 return self.name

1817

1818

1819class PGMError(GeoidError):

1820 '''An issue while parsing or cropping an C{egm*.pgm} geoid dataset.

1821 '''

1822 pass

1823

1824

1825def egmGeoidHeights(GeoidHeights_dat):

1826 '''Generate geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/

1827 C++/doc/geoid.html#geoidinst>} height tests from U{GeoidHeights.dat

1828 <https://SourceForge.net/projects/geographiclib/files/testdata/>}

1829 U{Test data for Geoids<https://GeographicLib.SourceForge.io/C++/doc/

1830 geoid.html#testgeoid>}.

1831

1832 @arg GeoidHeights_dat: The un-gz-ed C{GeoidHeights.dat} file

1833 (C{str} or C{file} handle).

1834

1835 @return: For each test, yield a L{GeoidHeight5Tuple}C{(lat, lon,

1836 egm84, egm96, egm2008)}.

1837

1838 @raise GeoidError: Invalid B{C{GeoidHeights_dat}}.

1839

1840 @note: Function L{egmGeoidHeights} is used to test the geoids

1841 L{GeoidKarney} and L{GeoidPGM}, see PyGeodesy module

1842 C{test/testGeoids.py}.

1843 '''

1844 dat = GeoidHeights_dat

1845 if isinstance(dat, bytes):

1846 dat = _BytesIO(dat)

1847

1848 try:

1849 dat.seek(0, _os.SEEK_SET) # reset

1850 except AttributeError as x:

1851 raise GeoidError(GeoidHeights_dat=type(dat), cause=x)

1852

1853 for t in dat.readlines():

1854 t = t.strip()

1855 if t and not t.startswith(_bHASH_):

1856 lat, lon, egm84, egm96, egm2008 = map(float, t.split())

1857 lon = _lonE2lon(lon) # Eastern to earth lon

1858 yield GeoidHeight5Tuple(lat, lon, egm84, egm96, egm2008)

1859

1860

1861def _lonE2lon(lon):

1862 '''(INTERNAL) East to earth longitude.

1863 '''

1864 while lon > _180_0:

1865 lon -= _360_0

1866 return lon

1867

1868

1869def _lon2lonE(lon):

1870 '''(INTERNAL) Earth to East longitude.

1871 '''

1872 while lon < 0:

1873 lon += _360_0

1874 return lon

1875

1876

1877__all__ += _ALL_DOCS(_GeoidBase)

1878

1879if __name__ == _DMAIN_: # MCCABE 14

1880

1881 from pygeodesy.internals import printf, _secs2str, _versions, _sys

1882 from time import time

1883

1884 _crop = {}

1885 _GeoidEGM = GeoidKarney

1886 _kind = 3

1887

1888 geoids = _sys.argv[1:]

1889 while geoids:

1890 G = geoids.pop(0)

1891 g = G.lower()

1892

1893 if '-crop'.startswith(g):

1894 _crop = dict(crp=(20, -125, 50, -65)) # CONUS

1895

1896 elif '-egm96'.startswith(g):

1897 _GeoidEGM = GeoidEGM96

1898

1899 elif '-karney'.startswith(g):

1900 _GeoidEGM = GeoidKarney

1901

1902 elif '-kind'.startswith(g):

1903 _kind = int(geoids.pop(0))

1904

1905 elif '-pgm'.startswith(g):

1906 _GeoidEGM = GeoidPGM

1907

1908 elif _isin(g[-4:], '.pgm', '.grd'):

1909 g = _GeoidEGM(G, kind=_kind, **_crop)

1910 t = time()

1911 _ = g.highest()

1912 t = _secs2str(time() - t)

1913 printf('%s: %s (%s)', g.toStr(), t, _versions(), nl=1, nt=1)

1914 t = g.pgm

1915 if t:

1916 printf(repr(t), nt=1)

1917 # <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>:

1918 # The height of the EGM96 geoid at Timbuktu

1919 # echo 16:46:33N 3:00:34W | GeoidEval

1920 # => 28.7068 -0.02e-6 -1.73e-6

1921 # The 1st number is the height of the geoid, the 2nd and

1922 # 3rd are its slopes in northerly and easterly direction

1923 t = 'Timbuktu %s' % (g,)

1924 k = {'egm84-15.pgm': '31.2979',

1925 'egm96-5.pgm': '28.7067',

1926 'egm2008-1.pgm': '28.7880'}.get(g.name.lower(), '28.7880')

1927 ll = _MODS.dms.parseDMS2('16:46:33N', '3:00:34W', sep=':')

1928 for ll in (ll, (16.776, -3.009),):

1929 try:

1930 h, ll = g.height(*ll), fstr(ll, prec=6)

1931 printf('%s.height(%s): %.4F vs %s', t, ll, h, k)

1932 except (GeoidError, RangeError) as x:

1933 printf(_COLONSPACE_(t, str(x)))

1934

1935 elif _isin(g[-4:], '.bin'):

1936 g = GeoidG2012B(G, kind=_kind)

1937 printf(g.toStr())

1938

1939 else:

1940 raise GeoidError(grid=repr(G))

1941

1942_I = int # PYCHOK unused _I

1943del _intCs, _T, _T0s12 # trash ints cache and map

1944

1945

1946# <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>

1947# _lowerleft = -90, -179, -30.1500 # egm2008-1.pgm

1948# _lowerleft = -90, -179, -29.5350 # egm96-5.pgm

1949# _lowerleft = -90, -179, -29.7120 # egm84-15.pgm

1950

1951# _center = 0, 0, 17.2260 # egm2008-1.pgm

1952# _center = 0, 0, 17.1630 # egm96-5.pgm

1953# _center = 0, 0, 18.3296 # egm84-15.pgm

1954

1955# _upperright = 90, 180, 14.8980 # egm2008-1.pgm

1956# _upperright = 90, 180, 13.6050 # egm96-5.pgm

1957# _upperright = 90, 180, 13.0980 # egm84-15.pgm

1958

1959

1960# % python3.13 -m pygeodesy.geoids -egm96 ../testGeoids/WW15MGH.GRD

1961#

1962# GeoidEGM96('WW15MGH.GRD'): kind=3, smooth=0, dtype=dtype('float64'), endian='tbd', mean=-1.426, nBytes=8311688, nots=1038961, shape=(721, 1441), sizeB=8153505, stdev=29.223, numpy=2.5.0, scipy=1.18.0, center=(0.0, 0.125, 17.125), highest=(-8.25, -32.75, 85.391), lowest=(4.75, -101.25, -106.991), lowerleft=(-90.0, -180.0, -29.534), loweright=(-90.0, 180.0, -29.534), upperleft=(90.0, -180.0, 13.606), upperright=(90.0, 180.0, 13.606): 598.907 us (pygeodesy 26.8.26 Python 3.13.13 64bit arm64 macOS 26.6.2)

1963#

1964# Timbuktu GeoidEGM96('WW15MGH.GRD').height(16.775833, -3.009444): 28.7073 vs 28.7880

1965# Timbuktu GeoidEGM96('WW15MGH.GRD').height(16.776, -3.009): 28.7072 vs 28.7880

1966

1967

1968# % python3.13 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm

1969#

1970# GeoidKarney('egm2008-1.pgm'): kind=3, cropped=False, dtype='ushort', endian='>H', hits=0, mean=-1.317, nots=233301600, shape=(10801, 21600), sizeB=466603604, stdev=29.244, center=(0.0, 0.0, 17.226), highest=(-8.4, 147.367, 85.839), lowest=(4.7, 78.767, -106.911), lowerleft=(-90.0, -180.0, -30.15), loweright=(-90.0, 180.0, -30.15), upperleft=(90.0, -180.0, 14.898), upperright=(90.0, 180.0, 14.898): 112.521 ms (pygeodesy 26.8.26 Python 3.13.13 64bit arm64 macOS 26.6.2)

1971#

1972# _PGM('../testGeoids/egm2008-1.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-31 06:54:00', Description='WGS84 EGM2008, 1-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.025, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.001, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008', Vertical_Datum='WGS84'

1973#

1974# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880

1975# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880

1976#

1977# GeoidKarney('egm84-30.pgm'): kind=3, cropped=False, dtype='ushort', endian='>H', hits=0, mean=-1.317, nots=259920, shape=(361, 720), sizeB=520256, stdev=29.244, center=(0.0, 0.0, 18.327), highest=(-4.5, 149.0, 81.33), lowest=(5.0, 79.0, -107.232), lowerleft=(-90.0, -180.0, -29.711), loweright=(-90.0, 180.0, -29.711), upperleft=(90.0, -180.0, 13.098), upperright=(90.0, 180.0, 13.098): 179.052 us (pygeodesy 26.8.26 Python 3.13.13 64bit arm64 macOS 26.6.2)

1978#

1979# _PGM('../testGeoids/egm84-30.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 30-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=1.546, MaxCubicError=0.274, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.07, RMSCubicError=0.014, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84'

1980#

1981# Timbuktu GeoidKarney('egm84-30.pgm').height(16.775833, -3.009444): 31.3031 vs 28.7880

1982# Timbuktu GeoidKarney('egm84-30.pgm').height(16.776, -3.009): 31.3027 vs 28.7880

1983#

1984# GeoidKarney('egm96-5.pgm'): kind=3, cropped=False, dtype='ushort', endian='>H', hits=0, mean=-1.317, nots=9335520, shape=(2161, 4320), sizeB=18671448, stdev=29.244, center=(0.0, 0.0, 17.163), highest=(-8.167, 147.25, 85.422), lowest=(4.667, 78.833, -107.043), lowerleft=(-90.0, -180.0, -29.535), loweright=(-90.0, 180.0, -29.535), upperleft=(90.0, -180.0, 13.605), upperright=(90.0, 180.0, 13.605): 4.645 ms (pygeodesy 26.8.26 Python 3.13.13 64bit arm64 macOS 26.6.2)

1985#

1986# _PGM('../testGeoids/egm96-5.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:03', Description='WGS84 EGM96, 5-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.14, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.005, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html', Vertical_Datum='WGS84'

1987#

1988# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067

1989# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067

1990

1991

1992# % python3.13 -m pygeodesy.geoids -PGM ../testGeoids/egm*.pgm

1993#

1994# GeoidPGM('egm2008-1.pgm'): kind=3, smooth=0, cropped=False, dtype=dtype('float64'), endian='>u2', mean=-1.317, nBytes=1866412800, nots=233301600, shape=(10801, 21600), sizeB=466603604, stdev=29.244, numpy=2.5.0, scipy=1.18.0, center=(0.0, 0.0, 17.226), highest=(-8.4, -32.633, 85.839), lowest=(4.683, -101.25, -106.911), lowerleft=(-90.0, -180.0, -30.15), loweright=(-90.0, 179.983, -30.15), upperleft=(90.0, -180.0, 14.898), upperright=(90.0, 179.983, 14.898): 574.773 ms (pygeodesy 26.8.26 Python 3.13.13 64bit arm64 macOS 26.6.2)

1995#

1996# _PGM('../testGeoids/egm2008-1.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-31 06:54:00', Description='WGS84 EGM2008, 1-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.025, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.001, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008', Vertical_Datum='WGS84'

1997#

1998# Timbuktu GeoidPGM('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880

1999# Timbuktu GeoidPGM('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880

2000#

2001# GeoidPGM('egm84-30.pgm'): kind=3, smooth=0, cropped=False, dtype=dtype('float64'), endian='>u2', mean=-0.865, nBytes=2079360, nots=259920, shape=(361, 720), sizeB=520256, stdev=29.175, numpy=2.5.0, scipy=1.18.0, center=(0.0, 0.0, 18.33), highest=(-4.5, -31.0, 81.33), lowest=(5.0, -101.0, -107.232), lowerleft=(-90.0, -180.0, -29.712), loweright=(-90.0, 179.5, -29.712), upperleft=(90.0, -180.0, 13.098), upperright=(90.0, 179.5, 13.098): 248.909 us (pygeodesy 26.8.26 Python 3.13.13 64bit arm64 macOS 26.6.2)

2002#

2003# _PGM('../testGeoids/egm84-30.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 30-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=1.546, MaxCubicError=0.274, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.07, RMSCubicError=0.014, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84'

2004#

2005# Timbuktu GeoidPGM('egm84-30.pgm').height(16.775833, -3.009444): 31.3010 vs 28.7880

2006# Timbuktu GeoidPGM('egm84-30.pgm').height(16.776, -3.009): 31.3006 vs 28.7880

2007#

2008# GeoidPGM('egm96-5.pgm'): kind=3, smooth=0, cropped=False, dtype=dtype('float64'), endian='>u2', mean=-1.438, nBytes=74684160, nots=9335520, shape=(2161, 4320), sizeB=18671448, stdev=29.227, numpy=2.5.0, scipy=1.18.0, center=(0.0, -0.0, 17.179), highest=(-8.167, -32.75, 85.422), lowest=(4.667, -101.167, -107.043), lowerleft=(-90.0, -180.0, -29.535), loweright=(-90.0, 179.917, -29.535), upperleft=(90.0, -180.0, 13.605), upperright=(90.0, 179.917, 13.605): 7.860 ms (pygeodesy 26.8.26 Python 3.13.13 64bit arm64 macOS 26.6.2)

2009#

2010# _PGM('../testGeoids/egm96-5.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:03', Description='WGS84 EGM96, 5-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.14, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.005, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html', Vertical_Datum='WGS84'

2011#

2012# Timbuktu GeoidPGM('egm96-5.pgm').height(16.775833, -3.009444): 28.7065 vs 28.7067

2013# Timbuktu GeoidPGM('egm96-5.pgm').height(16.776, -3.009): 28.7064 vs 28.7067

2014

2015

2016# % python2 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm

2017#

2018# GeoidKarney('egm2008-1.pgm'): kind=3, cropped=False, dtype='ushort', endian='>H', hits=0, mean=-1.317, nots=233301600, shape=(10801, 21600), sizeB=466603604, stdev=29.244, center=(0.0, 0.0, 17.226), highest=(-8.4, 147.367, 85.839), lowest=(4.7, 78.767, -106.911), lowerleft=(-90.0, -180.0, -30.15), loweright=(-90.0, 180.0, -30.15), upperleft=(90.0, -180.0, 14.898), upperright=(90.0, 180.0, 14.898): 167.144 ms (pygeodesy 26.8.26 Python 2.7.18 64bit arm64_x86_64 macOS 26.6.2)

2019#

2020# _PGM('../testGeoids/egm2008-1.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-31 06:54:00', Description='WGS84 EGM2008, 1-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.025, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.001, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008', Vertical_Datum='WGS84'

2021#

2022# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880

2023# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880

2024#

2025# GeoidKarney('egm84-30.pgm'): kind=3, cropped=False, dtype='ushort', endian='>H', hits=0, mean=-1.317, nots=259920, shape=(361, 720), sizeB=520256, stdev=29.244, center=(0.0, 0.0, 18.327), highest=(-4.5, 149.0, 81.33), lowest=(5.0, 79.0, -107.232), lowerleft=(-90.0, -180.0, -29.711), loweright=(-90.0, 180.0, -29.711), upperleft=(90.0, -180.0, 13.098), upperright=(90.0, 180.0, 13.098): 227.928 us (pygeodesy 26.8.26 Python 2.7.18 64bit arm64_x86_64 macOS 26.6.2)

2026#

2027# _PGM('../testGeoids/egm84-30.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 30-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=1.546, MaxCubicError=0.274, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.07, RMSCubicError=0.014, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84'

2028#

2029# Timbuktu GeoidKarney('egm84-30.pgm').height(16.775833, -3.009444): 31.3031 vs 28.7880

2030# Timbuktu GeoidKarney('egm84-30.pgm').height(16.776, -3.009): 31.3027 vs 28.7880

2031#

2032# GeoidKarney('egm96-5.pgm'): kind=3, cropped=False, dtype='ushort', endian='>H', hits=0, mean=-1.317, nots=9335520, shape=(2161, 4320), sizeB=18671448, stdev=29.244, center=(0.0, 0.0, 17.163), highest=(-8.167, 147.25, 85.422), lowest=(4.667, 78.833, -107.043), lowerleft=(-90.0, -180.0, -29.535), loweright=(-90.0, 180.0, -29.535), upperleft=(90.0, -180.0, 13.605), upperright=(90.0, 180.0, 13.605): 6.938 ms (pygeodesy 26.8.26 Python 2.7.18 64bit arm64_x86_64 macOS 26.6.2)

2033#

2034# _PGM('../testGeoids/egm96-5.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:03', Description='WGS84 EGM96, 5-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.14, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.005, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html', Vertical_Datum='WGS84'

2035#

2036# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067

2037# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067

2038

2039

2040# **) MIT License

2041#

2043#

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

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

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

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

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

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

2050#

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

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

2053#

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

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

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

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

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

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

2060# OTHER DEALINGS IN THE SOFTWARE.

Coverage for pygeodesy / geoids.py: 96%

713 statements