Coverage for pygeodesy/geoids.py: 96%
688 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-03-25 15:01 -0400
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 _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 Property_RO, property_RO, property_ROver
112from pygeodesy.streprs import attrs, Fmt, fstr, pairs
113from pygeodesy.units import Height, Int_, Lat, Lon
114# from pygeodesy.utily import _Wrap # from .heights
116from math import floor as _floor
117# import os as _os # _MODS
118# import os.path # _os.path
119from struct import calcsize as _calcsize, unpack as _unpack
120try:
121 from StringIO import StringIO as _BytesIO # reads bytes
122 _ub2str = str # PYCHOK convert bytes to str for egm*.pgm text
124except ImportError: # Python 3+
125 from io import BytesIO as _BytesIO # PYCHOK expected
126 from pygeodesy.basics import ub2str as _ub2str
128__all__ = _ALL_LAZY.geoids
129__version__ = '26.02.02'
131_assert_ = 'assert'
132_bHASH_ = b'#'
133_endian_ = 'endian'
134_format_ = '%s %r'
135_header_ = 'header'
136_intCs = {} # cache int value, del below
137_lli_ = 'lli'
138_os = _MODS.os
139_rb_ = 'rb'
140_supported_ = 'supported'
142if __debug__:
143 from pygeodesy.fmath import Fdot as _Dotf, Fhorner as _Hornerf
145else: # -OO ... runs GeoidKarney 8+X faster (w/o Fwelford)
146 from pygeodesy.fsums import _fsum
147 from operator import mul as _mul
149 class _GKsum(Fsum): # for GeoidKarney only
151 def __iadd__(self, other):
152 xs = other._ps if isinstance(other, _GKsum) else (other,)
153 self._ps_acc(self._ps, xs, up=False)
154 return self
156 def __imul__(self, x):
157 self._ps[:] = (x * p for p in self._ps)
158 return self
160 fmul = __imul__ # like .fsums.Fsum
162 def fsum(self): # PYCHOK signature
163 return _fsum(self._ps)
165 class _Dotf(_GKsum): # PYCHOK redef
166 '''Fast precision dot product M{sum(a[i] * b[i] for i=0..len(a))}
167 but for C{float} C{a} and C{b} only.
168 '''
169 def __init__(self, a, *b):
170 self._ps = self._ps_acc([], map(_mul, a, b), up=False)
172 class _Hornerf(_GKsum): # PYCHOK redef
173 '''Fast polynomial evaluation M{sum(cs[i] * x**i for i=0..len(cs))}
174 but for C{float} C{x} and C{cs} and C{incx=True} only.
175 '''
176 def __init__(self, x, *cs): # incx=True
177 self._ps = []
178 for c in reversed(cs): # multiply-accumulate
179 self *= x # ps[:] = (x * p for p in ps)
180 self += c # self._ps_acc(ps, (c,), up=False)
183class GeoidError(HeightError):
184 '''Geoid interpolator C{Geoid...} or interpolation issue.
185 '''
186 pass
189class _GeoidBase(_HeightBase):
190 '''(INTERNAL) Base class for C{Geoid...}s.
191 '''
192# _center = None # (lat, lon, height)
193 _cropped = None
194# _datum = _WGS84 # from _HeightBase
195 _egm = None # open C{egm*.pgm} geoid file
196 _endian = _tbd_
197 _Error = GeoidError # in ._HeightBase._as_lls, ...
198 _geoid = _n_a_
199 _hs_y_x = None # numpy 2darray, row-major order
200 _iscipy = True # scipy or Karney's interpolation
201 _kind = 3 # order for interp2d, RectBivariateSpline
202# _kmin = 2 # min number of knots
203 _knots = 0 # nlat * nlon
204 _mean = None # fixed in GeoidKarney
205 _nBytes = 0 # numpy size in bytes, float64
206 _pgm = None # PGM attributes, C{_PGM} or C{None}
207 _sizeB = 0 # geoid file size in bytes
208 _smooth = 0 # used only for RectBivariateSpline
209 _stdev = None # fixed in GeoidKarney
210 _u2B = 0 # np.itemsize or undefined
211 _yx_hits = None # cache hits, ala Karney's
213# _lat_d = _0_0 # increment, +tive
214# _lat_lo = _0_0 # lower lat, south
215# _lat_hi = _0_0 # upper lat, noth
216# _lon_d = _0_0 # increment, +tive
217# _lon_lo = _0_0 # left lon, west
218# _lon_hi = _0_0 # right lon, east
219# _lon_of = _0_0 # forward lon offset
220# _lon_og = _0_0 # reverse lon offset
222 def __init__(self, hs, p):
223 '''(INTERNAL) Set up the grid axes, the C{SciPy} interpolator and
224 several internal geoid attributes.
226 @arg hs: Grid knots with known height (C{numpy 2darray}).
227 @arg p: The C{slat, wlon, nlat, nlon, dlat, dlon} and other
228 geoid parameters (C{INTERNAL}).
229 '''
230 spi = self.scipy_interpolate
231 # for 2d scipy.interpolate.interp2d(xs, ys, hs, ...) and
232 # scipy.interpolate.RectBivariateSpline(ys, xs, hs, ...)
233 # require the shape of hs to be (len(ys), len(xs)), note
234 # the different (xs, ys, ...) and (ys, xs, ...) orders
235 if (p.nlat, p.nlon) != hs.shape:
236 raise GeoidError(shape=hs.shape, txt=_incompatible((p.nlat, p.nlon)))
238 # both axes and bounding box
239 ys, self._lat_d = self._gaxis2(p.slat, p.dlat, p.nlat, _lat_ + _s_)
240 xs, self._lon_d = self._gaxis2(p.wlon, p.dlon, p.nlon, _lon_ + _s_)
242 bb = ys[0], ys[-1], xs[0], xs[-1] + p.dlon # fudge lon_hi
243 # geoid grids are typically stored in row-major order, some
244 # with rows (90..-90) reversed and columns (0..360) wrapped
245 # to Easten longitude, 0 <= east < 180 and 180 <= west < 360
246 k = self.kind
247 if k in self._k2interp2d: # see _HeightBase
248 self._interp2d(xs, ys, hs, kind=k)
249 else: # XXX order ys and xs, see HeightLSQBiSpline
250 k = self._kxky(k)
251 self._ev = spi.RectBivariateSpline(ys, xs, hs, bbox=bb, ky=k, kx=k,
252 s=self._smooth).ev
253 self._hs_y_x = hs # numpy 2darray, row-major
254 self._nBytes = hs.nbytes # numpy size in bytes
255 self._knots = p.knots # grid knots == len(hs)
256 self._lon_of = float(p.flon) # forward offset
257 self._lon_og = g = float(p.glon) # reverse offset
258 # shrink the bounding box by 1 unit on every side:
259 # +self._lat_d, -self._lat_d, +self._lon_d, -self._lon_d
260 self._lat_lo, \
261 self._lat_hi, \
262 self._lon_lo, \
263 self._lon_hi = map(float, bb)
264 self._lon_lo -= g
265 self._lon_hi -= g
267 def __call__(self, *llis, **wrap_H):
268 '''Interpolate the geoid (or orthometric) height for one or more locations.
270 @arg llis: One or several locations (each C{LatLon}), all positional.
271 @kwarg wrap_H: Keyword arguments C{B{wrap}=False} (C{bool}) and
272 C{B{H}=False} (C{bool}). Use C{B{wrap}=True} to wrap
273 or I{normalize} all B{C{llis}} locations. If C{B{H}
274 is True}, return the I{orthometric} height instead of
275 the I{geoid} height at each location.
277 @return: A single geoid (or orthometric) height (C{float}) or
278 a list or tuple of geoid (or orthometric) heights (each
279 C{float}).
281 @raise GeoidError: Insufficient number of B{C{llis}}, an invalid
282 B{C{lli}} or the C{egm*.pgm} geoid file is closed.
284 @raise RangeError: An B{C{lli}} is outside this geoid's lat- or
285 longitude range.
287 @raise SciPyError: A C{scipy} issue.
289 @raise SciPyWarning: A C{scipy} warning as exception.
291 @note: To obtain I{orthometric} heights, each B{C{llis}} location
292 must have an ellipsoid C{height} or C{h} attribute, otherwise
293 C{height=0} is used.
295 @see: Function L{pygeodesy.heightOrthometric}.
296 '''
297 return self._called(llis, **wrap_H)
299 def __enter__(self):
300 '''Open context.
301 '''
302 return self
304 def __exit__(self, *unused): # PYCHOK exc_type, exc_value, exc_traceback)
305 '''Close context.
306 '''
307 self.close()
308 # return None # XXX False
310 def __repr__(self):
311 return self.toStr()
313 def __str__(self):
314 return Fmt.PAREN(self.classname, repr(self.name))
316 def _called(self, llis, wrap=False, H=False):
317 # handle __call__
318 _H = self._heightOrthometric if H else None
319 _as, llis = _as_llis2(llis, Error=GeoidError)
320 _w, hs = _Wrap._latlonop(wrap), []
321 _h, _a = self._hGeoid, hs.append
322 try:
323 for i, lli in enumerate(llis):
324 N = _h(*_w(lli.lat, lli.lon))
325 # orthometric or geoid height
326 _a(_H(lli, N) if _H else N)
327 return _as(hs)
328 except (GeoidError, RangeError) as x:
329 # XXX avoid str(LatLon()) degree symbols
330 n = _lli_ if _as is _ascalar else Fmt.INDEX(llis=i)
331 t = fstr((lli.lat, lli.lon), strepr=repr)
332 raise type(x)(n, t, wrap=wrap, H=H, cause=x)
333 except Exception as x:
334 if self._iscipy and self.scipy:
335 raise _SciPyIssue(x, self._ev_name)
336 else:
337 raise
339 @Property_RO
340 def _center(self):
341 '''(INTERNAL) Cache for method L{center}.
342 '''
343 return self._llh3(favg(self._lat_lo, self._lat_hi),
344 favg(self._lon_lo, self._lon_hi))
346 def center(self, LatLon=None):
347 '''Return the center location and height of this geoid.
349 @kwarg LatLon: Optional class to return the location and height
350 (C{LatLon}) or C{None}.
352 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,
353 height)} otherwise a B{C{LatLon}} instance with the lat-,
354 longitude and geoid height of the center grid location.
355 '''
356 return self._llh3LL(self._center, LatLon)
358 def close(self):
359 '''Close the C{egm*.pgm} geoid file if open (and applicable).
360 '''
361 if not self.closed:
362 self._egm.close()
363 self._egm = None
365 @property_RO
366 def closed(self):
367 '''Get the C{egm*.pgm} geoid file status.
368 '''
369 return self._egm is None
371 @Property_RO
372 def cropped(self):
373 '''Is geoid cropped (C{bool} or C{None} if crop not supported).
374 '''
375 return self._cropped
377 @Property_RO
378 def dtype(self):
379 '''Get the grid C{scipy} U{dtype<https://docs.SciPy.org/doc/numpy/
380 reference/generated/numpy.ndarray.dtype.html>} (C{numpy.dtype}).
381 '''
382 return self._hs_y_x.dtype
384 @Property_RO
385 def endian(self):
386 '''Get the geoid endianess and U{dtype<https://docs.SciPy.org/
387 doc/numpy/reference/generated/numpy.dtype.html>} (C{str}).
388 '''
389 return self._endian
391 def _ev(self, y, x): # PYCHOK overwritten with .RectBivariateSpline.ev
392 # see methods _HeightBase._ev and -._interp2d
393 return self._ev2d(x, y) # (y, x) flipped!
395 def _gaxis2(self, lo, d, n, name):
396 # build grid axis, hi = lo + (n - 1) * d
397 m, a = len2(frange(lo, n, d))
398 if m != n:
399 raise LenError(type(self), grid=m, **{name: n})
400 if d < 0:
401 d, a = -d, list(reversed(a))
402 a = self.numpy.array(a)
403 m, i = min2(*map(float, a[1:] - a[:-1]))
404 if m < EPS: # non-increasing axis
405 i = Fmt.INDEX(name, i + 1)
406 raise GeoidError(i, m, txt_not_='increasing')
407 return a, d
409 def _g2ll2(self, lat, lon): # PYCHOK no cover
410 '''(INTERNAL) I{Must be overloaded}.'''
411 self._notOverloaded(lat, lon)
413 def _gyx2g2(self, y, x):
414 # convert grid (y, x) indices to grid (lat, lon)
415 return ((self._lat_lo + self._lat_d * y),
416 (self._lon_lo + self._lon_of + self._lon_d * x))
418 def height(self, lats, lons, **wrap):
419 '''Interpolate the geoid height for one or several lat-/longitudes.
421 @arg lats: Latitude or latitudes (each C{degrees}).
422 @arg lons: Longitude or longitudes (each C{degrees}).
423 @kwarg wrap: Use C{B{wrap}=True} to wrap or I{normalize} all
424 B{C{lats}} and B{C{lons}}.
426 @return: A single geoid height (C{float}) or a list of geoid
427 heights (each C{float}).
429 @raise GeoidError: Insufficient or unequal number of B{C{lats}}
430 and B{C{lons}}.
432 @raise RangeError: A B{C{lat}} or B{C{lon}} is outside this geoid's
433 lat- or longitude range.
435 @raise SciPyError: A C{scipy} issue.
437 @raise SciPyWarning: A C{scipy} warning as exception.
438 '''
439 return _HeightBase.height(self, lats, lons, **wrap)
441 def height_(self, *latlons, **wrap):
442 '''Interpolate the geoid height for each M{(latlons[i], latlons[i+1])
443 pair for i in range(0, len(latlons), B{2})}.
445 @arg latlons: Alternating lat-/longitude pairs (each C{degrees}),
446 all positional.
448 @see: Method L{height} for further details.
450 @return: A tuple of geoid heights (each C{float}).
451 '''
452 lls = tuple(self._as_lls(latlons[0::2], *latlons[1::2]))
453 return self._called(lls, **wrap)
455 @property_ROver
456 def _heightOrthometric(self):
457 return _MODS.formy.heightOrthometric # overwrite property_ROver
459 def _hGeoid(self, lat, lon):
460 out = self.outside(lat, lon)
461 if out: # XXX avoid str(LatLon()) degree symbols
462 t = fstr((lat, lon), strepr=repr)
463 raise RangeError(lli=t, txt=_SPACE_(_outside_, _on_, out))
464 return float(self._ev(*self._ll2g2(lat, lon)))
466 @Property_RO
467 def _highest(self):
468 '''(INTERNAL) Cache for C{.highest}.
469 '''
470 return self._LL3T(self._llh3minmax(True), name__=self.highest)
472 def highest(self, LatLon=None, **unused):
473 '''Return the location and largest height of this geoid.
475 @kwarg LatLon: Optional class to return the location and height
476 (C{LatLon}) or C{None}.
478 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,
479 height)} otherwise a B{C{LatLon}} instance with the lat-,
480 longitude and geoid height of the highest grid location.
481 '''
482 return self._llh3LL(self._highest, LatLon)
484 @Property_RO
485 def hits(self):
486 '''Get the number of cache hits (C{int} or C{None}).
487 '''
488 return self._yx_hits
490 @Property_RO
491 def kind(self):
492 '''Get the interpolator kind and order (C{int}).
493 '''
494 return self._kind
496 @Property_RO
497 def knots(self):
498 '''Get the number of grid knots (C{int}).
499 '''
500 return self._knots
502 def _ll2g2(self, lat, lon): # PYCHOK no cover
503 '''(INTERNAL) I{Must be overloaded}.'''
504 self._notOverloaded(lat, lon)
506 @property_ROver
507 def _LL3T(self):
508 '''(INTERNAL) Get L{LatLon3Tuple}, I{once}.
509 '''
510 return _MODS.namedTuples.LatLon3Tuple # overwrite property_ROver
512 def _llh3(self, lat, lon):
513 return self._LL3T(lat, lon, self._hGeoid(lat, lon), name=self.name)
515 def _llh3LL(self, llh, LatLon):
516 return llh if LatLon is None else self._xnamed(LatLon(*llh))
518 def _llh3minmax(self, highest, *unused):
519 hs, np = self._hs_y_x, self.numpy
520 # <https://docs.SciPy.org/doc/numpy/reference/generated/
521 # numpy.argmin.html#numpy.argmin>
522 arg = np.argmax if highest else np.argmin
523 y, x = np.unravel_index(arg(hs, axis=None), hs.shape)
524 return self._g2ll2(*self._gyx2g2(y, x)) + (float(hs[y, x]),)
526 def _load(self, g, dtype=float, n=-1, offset=0, **sep): # sep=NN
527 # numpy.fromfile, like .frombuffer
528 g.seek(offset, _os.SEEK_SET)
529 return self.numpy.fromfile(g, dtype, count=n, **sep)
531 @Property_RO
532 def _lowerleft(self):
533 '''(INTERNAL) Cache for C{.lowerleft}.
534 '''
535 return self._llh3(self._lat_lo, self._lon_lo)
537 def lowerleft(self, LatLon=None):
538 '''Return the lower-left location and height of this geoid.
540 @kwarg LatLon: Optional class to return the location
541 (C{LatLon}) and height or C{None}.
543 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon, height)}
544 otherwise a B{C{LatLon}} instance with the lat-, longitude and
545 geoid height of the lower-left, SW grid corner.
546 '''
547 return self._llh3LL(self._lowerleft, LatLon)
549 @Property_RO
550 def _loweright(self):
551 '''(INTERNAL) Cache for C{.loweright}.
552 '''
553 return self._llh3(self._lat_lo, self._lon_hi)
555 def loweright(self, LatLon=None):
556 '''Return the lower-right location and height of this geoid.
558 @kwarg LatLon: Optional class to return the location and height
559 (C{LatLon}) or C{None}.
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-right, SE grid corner.
564 '''
565 return self._llh3LL(self._loweright, LatLon)
567 lowerright = loweright # synonymous
569 @Property_RO
570 def _lowest(self):
571 '''(INTERNAL) Cache for C{.lowest}.
572 '''
573 return self._LL3T(self._llh3minmax(False), name__=self.lowest)
575 def lowest(self, LatLon=None, **unused):
576 '''Return the location and lowest height of this geoid.
578 @kwarg LatLon: Optional class to return the location and height
579 (C{LatLon}) or C{None}.
581 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,
582 height)} otherwise a B{C{LatLon}} instance with the lat-,
583 longitude and geoid height of the lowest grid location.
584 '''
585 return self._llh3LL(self._lowest, LatLon)
587 @Property_RO
588 def mean(self):
589 '''Get the mean of this geoid's heights (C{float}).
590 '''
591 if self._mean is None: # see GeoidKarney
592 self._mean = float(self.numpy.mean(self._hs_y_x))
593 return self._mean
595 @property_RO
596 def name(self):
597 '''Get the name of this geoid (C{str}).
598 '''
599 return _HeightBase.name.fget(self) or self._geoid # recursion
601 @Property_RO
602 def nBytes(self):
603 '''Get the grid in-memory size in bytes (C{int}).
604 '''
605 return self._nBytes
607 def _open(self, geoid, datum, kind, name, smooth):
608 # open the geoid file
609 try:
610 self._geoid = _os.path.basename(geoid)
611 self._sizeB = _os.path.getsize(geoid)
612 g = open(geoid, _rb_)
613 except (IOError, OSError) as x:
614 raise GeoidError(geoid=geoid, cause=x)
616 if not _isin(datum, None, self._datum):
617 self._datum = _ellipsoidal_datum(datum, name=name)
618 self._kind = int(kind)
619 if name:
620 _HeightBase.name.fset(self, name) # rename
621 if smooth:
622 self._smooth = Int_(smooth=smooth, Error=GeoidError, low=0)
624 return g
626 def outside(self, lat, lon):
627 '''Check whether a location is outside this geoid's lat-/longitude
628 or crop range.
630 @arg lat: The latitude (C{degrees}).
631 @arg lon: The longitude (C{degrees}).
633 @return: A 1- or 2-character C{str} if outside, an empty C{str} otherwise.
634 '''
635 lat = _S_ if lat < self._lat_lo else (_N_ if lat > self._lat_hi else NN)
636 lon = _W_ if lon < self._lon_lo else (_E_ if lon > self._lon_hi else NN)
637 return NN(lat, lon) if lat and lon else (lat or lon)
639 @Property_RO
640 def pgm(self):
641 '''Get the PGM attributes (C{_PGM} or C{None} if not available/applicable).
642 '''
643 return self._pgm
645 @Property_RO
646 def sizeB(self):
647 '''Get the geoid grid file size in bytes (C{int}).
648 '''
649 return self._sizeB
651 @Property_RO
652 def smooth(self):
653 '''Get the C{RectBivariateSpline} smoothing (C{int}).
654 '''
655 return self._smooth
657 @Property_RO
658 def stdev(self):
659 '''Get the standard deviation of this geoid's heights (C{float}) or C{None}.
660 '''
661 if self._stdev is None: # see GeoidKarney
662 self._stdev = float(self.numpy.std(self._hs_y_x))
663 return self._stdev
665 def _swne(self, crop):
666 # crop box to 4-tuple (s, w, n, e)
667 try:
668 if len(crop) == 2:
669 try: # sw, ne LatLons
670 swne = (crop[0].lat, crop[0].lon,
671 crop[1].lat, crop[1].lon)
672 except AttributeError: # (s, w), (n, e)
673 swne = tuple(crop[0]) + tuple(crop[1])
674 else: # (s, w, n, e)
675 swne = crop
676 if len(swne) == 4:
677 s, w, n, e = map(float, swne)
678 if _N_90_0 <= s <= (n - _1_0) <= 89.0 and \
679 _N_180_0 <= w <= (e - _1_0) <= 179.0:
680 return s, w, n, e
681 except (IndexError, TypeError, ValueError):
682 pass
683 raise GeoidError(crop=crop)
685 def toStr(self, prec=3, sep=_COMMASPACE_): # PYCHOK signature
686 '''This geoid and all geoid attributes as a string.
688 @kwarg prec: Number of decimal digits (0..9 or C{None} for
689 default). Trailing zero decimals are stripped
690 for B{C{prec}} values of 1 and above, but kept
691 for negative B{C{prec}} values.
692 @kwarg sep: Separator to join (C{str}).
694 @return: Geoid name and attributes (C{str}).
695 '''
696 s = 1 if self.kind < 0 else 2
697 t = _MODS.internals.typename
698 t = tuple(Fmt.PAREN(t(m), fstr(m(), prec=prec)) for m in
699 (self.lowerleft, self.upperright,
700 self.center,
701 self.highest, self.lowest)) + \
702 attrs( _mean_, _stdev_, prec=prec, Nones=False) + \
703 attrs((_kind_, 'smooth')[:s], prec=prec, Nones=False) + \
704 attrs( 'cropped', 'dtype', _endian_, 'hits', 'knots', 'nBytes',
705 'sizeB', _scipy_, _numpy_, prec=prec, Nones=False)
706 return _COLONSPACE_(self, sep.join(t))
708 @Property_RO
709 def u2B(self):
710 '''Get the PGM itemsize in bytes (C{int}).
711 '''
712 return self._u2B
714 @Property_RO
715 def _upperleft(self):
716 '''(INTERNAL) Cache for C{.upperleft}.
717 '''
718 return self._llh3(self._lat_hi, self._lon_lo)
720 def upperleft(self, LatLon=None):
721 '''Return the upper-left location and height of this geoid.
723 @kwarg LatLon: Optional class to return the location and height
724 (C{LatLon}) or C{None}.
726 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon, height)}
727 otherwise a B{C{LatLon}} instance with the lat-, longitude and
728 geoid height of the upper-left, NW grid corner.
729 '''
730 return self._llh3LL(self._upperleft, LatLon)
732 @Property_RO
733 def _upperright(self):
734 '''(INTERNAL) Cache for C{.upperright}.
735 '''
736 return self._llh3(self._lat_hi, self._lon_hi)
738 def upperright(self, LatLon=None):
739 '''Return the upper-right location and height of this geoid.
741 @kwarg LatLon: Optional class to return the location and height
742 (C{LatLon}) or C{None}.
744 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon, height)}
745 otherwise a B{C{LatLon}} instance with the lat-, longitude and
746 geoid height of the upper-right, NE grid corner.
747 '''
748 return self._llh3LL(self._upperright, LatLon)
751class GeoidEGM96(_GeoidBase):
752 '''Geoid height interpolator for the EGM96 U{15 Minute Interpolation Grid<https://earth-info.NGA.mil>}
753 based on C{SciPy} interpolation U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/reference/
754 generated/scipy.interpolate.RectBivariateSpline.html>}, U{interp2d<https://docs.SciPy.org/doc/scipy/
755 reference/generated/scipy.interpolate.interp2d.html>} or U{bisplrep/-ev<https://docs.scipy.org/doc/
756 scipy/reference/generated/scipy.interpolate.bisplrep.html>}.
758 Use only the C{WW15MGH.GRD} file, unzipped from the EGM96 U{15 Minute Interpolation Grid
759 <https://earth-info.NGA.mil/index.php?dir=wgs84&action=wgs84>} download.
760 '''
761 def __init__(self, EGM96_grd, datum=_WGS84, kind=3, smooth=0, **name_crop):
762 '''New L{GeoidEGM96} interpolator.
764 @arg EGM96_grd: An C{EGM96_grd} grid file name (C{.GRD}).
765 @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}),
766 overriding C{WGS84}.
767 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline
768 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.
769 RectBivariateSpline.html>} or -1, -3 or -5 for U{bisplrep/-ev<https://
770 docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}
771 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
772 interpolate.interp2d.html>} C{linear}, C{cubic} respectively C{quintic},
773 see note for more details.
774 @kwarg smooth: Smoothing factor for C{B{kind}=1..5} only (C{int}).
775 @kwarg name_crop: Optional geoid C{B{name}=NN} (C{str}) and UNSUPPORTED keyword argument
776 C{B{crop}=None}.
778 @raise GeoidError: Invalid B{C{crop}}, B{C{kind}} or B{C{smooth}} or a ECM96 grid file
779 B{C{ECM96_grd}} issue.
781 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.
783 @raise LenError: Grid file B{C{EGM96_grd}} axis mismatch.
785 @raise SciPyError: A C{scipy} issue.
787 @raise SciPyWarning: A C{scipy} warning as exception.
789 @raise TypeError: Invalid B{C{datum}}.
791 @note: Specify C{B{kind}=-1, -3 or -5} to use C{scipy.interpolate.interp2d}
792 before or C{scipy.interpolate.bisplrep/-ev} since C{Scipy} version 1.14.
793 '''
794 crop, name = _xkwds_pop2(name_crop, crop=None)
795 if crop is not None:
796 raise GeoidError(crop=crop, txt_not_=_supported_)
798 g = self._open(EGM96_grd, datum, kind, _name__(**name), smooth)
799 _ = self.numpy # import numpy for .fromfile, .reshape
801 try:
802 p, hs = _Gpars(), self._load(g, sep=_SPACE_) # text
803 p.slat, n, p.wlon, e, p.dlat, p.dlon = hs[:6] # n-s, 0-E
804 p.nlat = int((n - p.slat) / p.dlat) + 1 # include S
805 p.nlon = int((e - p.wlon) / p.dlon) + 1 # include W
806 p.knots = p.nlat * p.nlon # inverted lats N downto S
807 p.glon = _180_0 # Eastern lons 0-360
808 hs = hs[6:].reshape(p.nlat, p.nlon)
809 _GeoidBase.__init__(self, hs, p)
811 except Exception as x:
812 raise _SciPyIssue(x, _in_, repr(EGM96_grd))
813 finally:
814 g.close()
816 def _g2ll2(self, lat, lon):
817 # convert grid (lat, lon) to earth (lat, lon)
818 return -lat, _lonE2lon(lon) # invert lat
820 def _ll2g2(self, lat, lon):
821 # convert earth (lat, lon) to grid (lat, lon)
822 return -lat, _lon2lonE(lon) # invert lat
824 if _FOR_DOCS:
825 __call__ = _GeoidBase.__call__
826 height = _GeoidBase.height
829class GeoidG2012B(_GeoidBase):
830 '''Geoid height interpolator for U{GEOID12B Model
831 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/>} grids U{CONUS
832 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_CONUS.shtml>},
833 U{Alaska<https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_AK.shtml>},
834 U{Hawaii<https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_HI.shtml>},
835 U{Guam and Northern Mariana Islands
836 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_GMNI.shtml>},
837 U{Puerto Rico and U.S. Virgin Islands
838 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_PRVI.shtml>} and
839 U{American Samoa<https://Geodesy.NOAA.gov/GEOID/GEOID12B/GEOID12B_AS.shtml>}
840 based on C{SciPy} interpolation U{RectBivariateSpline<https://docs.SciPy.org/doc/
841 scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>}, U{interp2d
842 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
843 or U{bisplrep/-ev<https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.
844 bisplrep.html>}.
845 '''
846 _datum = Datums.NAD83
848 def __init__(self, g2012b_bin, datum=Datums.NAD83, kind=3, smooth=0, **name_crop):
849 '''New L{GeoidG2012B} interpolator.
851 @arg g2012b_bin: A C{GEOID12B} grid file name (C{.bin}, see B{note}).
852 @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
853 L{a_f2Tuple}), overriding C{NAD83}.
854 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline
855 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.
856 RectBivariateSpline.html>} or -1, -3 or -5 for U{bisplrep/-ev<https://
857 docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}
858 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
859 interpolate.interp2d.html>} C{linear}, C{cubic} respectively C{quintic},
860 see note for more details.
861 @kwarg smooth: Smoothing factor for C{B{kind}=1..5} only (C{int}).
862 @kwarg name_crop: Optional geoid C{B{name}=NN} (C{str}) and UNSUPPORTED keyword argument
863 C{B{crop}=None}.
865 @raise GeoidError: Invalid B{C{crop}}, B{C{kind}} or B{C{smooth}} or a B{C{g2012b_bin}}
866 grid file issue.
868 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.
870 @raise LenError: Grid file B{C{g2012b_bin}} axis mismatch.
872 @raise SciPyError: A C{scipy} issue.
874 @raise SciPyWarning: A C{scipy} warning as exception.
876 @raise TypeError: Invalid B{C{datum}}.
878 @note: Only use any of the C{le} (little-endian) or C{be} (big-endian) C{g2012b*.bin}
879 I{binary} grid files.
881 @note: Specify C{B{kind}=-1, -3 or -5} to use C{scipy.interpolate.interp2d} I{before}
882 or C{scipy.interpolate.bisplrep/-ev} I{since} C{Scipy} version 1.14.
883 '''
884 crop, name = _xkwds_pop2(name_crop, crop=None)
885 if crop is not None:
886 raise GeoidError(crop=crop, txt_not_=_supported_)
888 g = self._open(g2012b_bin, datum, kind, _name__(**name), smooth)
889 _ = self.numpy # import numpy for ._load and
891 try:
892 p = _Gpars()
893 n = (self.sizeB // 4) - 11 # number of f4 heights
894 # U{numpy dtype formats are different from Python struct formats
895 # <https://docs.SciPy.org/doc/numpy-1.15.0/reference/arrays.dtypes.html>}
896 for en_ in ('<', '>'):
897 # skip 4xf8, get 3xi4
898 p.nlat, p.nlon, ien = map(int, self._load(g, en_+'i4', 3, 32))
899 if ien == 1: # correct endian
900 p.knots = p.nlat * p.nlon
901 if p.knots == n and 1 < p.nlat < n \
902 and 1 < p.nlon < n:
903 self._endian = en_+'f4'
904 break
905 else: # couldn't validate endian
906 raise GeoidError(_endian_)
908 # get the first 4xf8
909 p.slat, p.wlon, p.dlat, p.dlon = map(float, self._load(g, en_+'f8', 4))
910 # read all f4 heights, ignoring the first 4xf8 and 3xi4
911 hs = self._load(g, self._endian, n, 44).reshape(p.nlat, p.nlon)
912 p.wlon -= _360_0 # western-most lon XXX _lonE2lon
913 _GeoidBase.__init__(self, hs, p)
915 except Exception as x:
916 raise _SciPyIssue(x, _in_, repr(g2012b_bin))
917 finally:
918 g.close()
920 def _g2ll2(self, lat, lon):
921 # convert grid (lat, lon) to earth (lat, lon)
922 return lat, lon
924 def _ll2g2(self, lat, lon):
925 # convert earth (lat, lon) to grid (lat, lon)
926 return lat, lon
928 if _FOR_DOCS:
929 __call__ = _GeoidBase.__call__
930 height = _GeoidBase.height
933class GeoidHeight5Tuple(_NamedTuple): # .geoids.py
934 '''5-Tuple C{(lat, lon, egm84, egm96, egm2008)} for U{GeoidHeights.dat
935 <https://SourceForge.net/projects/geographiclib/files/testdata/>}
936 tests with the heights for 3 different EGM grids at C{degrees90}
937 and C{degrees180} degrees (after converting C{lon} from original
938 C{0 <= EasterLon <= 360}).
939 '''
940 _Names_ = (_lat_, _lon_, 'egm84', 'egm96', 'egm2008')
941 _Units_ = ( Lat, Lon, Height, Height, Height)
944def _I(i):
945 '''(INTERNAL) Cache a single C{int} constant.
946 '''
947 i = int(i)
948 return _intCs.setdefault(i, i) # noqa: F821 del
951def _T(cs):
952 '''(INTERNAL) Cache a tuple of single C{int} constants.
953 '''
954 return tuple(map(_I, _splituple(cs)))
956_T0s12 = (_I(0),) * 12 # PYCHOK _T('0, 0, ..., 0')
959class GeoidKarney(_GeoidBase):
960 '''Geoid height interpolator for I{Karney}'s U{GeographicLib Earth
961 Gravitational Model (EGM)<https://GeographicLib.SourceForge.io/C++/doc/
962 geoid.html>} geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/
963 C++/doc/geoid.html#geoidinst>} datasets using bilinear or U{cubic
964 <https://dl.ACM.org/citation.cfm?id=368443>} interpolation and U{caching
965 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidcache>}
966 in pure Python, transcoded from I{Karney}'s U{C++ class Geoid
967 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinterp>}.
969 Use any of the geoid U{egm84-, egm96- or egm2008-*.pgm
970 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>}
971 datasets.
972 '''
973 _C0 = _F(372), _F(240), _F(372) # n, _ and s common denominators
974 # matrices c3n_, c3, c3s_, transposed from GeographicLib/Geoid.cpp
975 _C3 = ((_T('0, 0, 62, 124, 124, 62, 0, 0, 0, 0, 0, 0'),
976 _T0s12,
977 _T('-131, 7, -31, -62, -62, -31, 45, 216, 156, -45, -55, -7'),
978 _T0s12,
979 _T('138, -138, 0, 0, 0, 0, -183, 33, 153, -3, 48, -48'), # PYCHOK indent
980 _T('144, 42, -62, -124, -124, -62, -9, 87, 99, 9, 42, -42'),
981 _T0s12,
982 _T('0, 0, 0, 0, 0, 0, 93, -93, -93, 93, 0, 0'),
983 _T('-102, 102, 0, 0, 0, 0, 18, 12, -12, -18, -84, 84'),
984 _T('-31, -31, 31, 62, 62, 31, 0, -93, -93, 0, 31, 31')), # PYCHOK indent
986 (_T('9, -9, 9, 186, 54, -9, -9, 54, -54, 9, -9, 9'),
987 _T('-18, 18, -88, -42, 162, -32, 8, -78, 78, -8, 18, -18'),
988 _T('-88, 8, -18, -42, -78, 18, 18, 162, 78, -18, -32, -8'),
989 _T('0, 0, 90, -150, 30, 30, 30, -90, 90, -30, 0, 0'),
990 _T('96, -96, 96, -96, -24, 24, -96, -24, 144, -24, 24, -24'), # PYCHOK indent
991 _T('90, 30, 0, -150, -90, 0, 0, 30, 90, 0, 30, -30'),
992 _T('0, 0, -20, 60, -60, 20, -20, 60, -60, 20, 0, 0'),
993 _T('0, 0, -60, 60, 60, -60, 60, -60, -60, 60, 0, 0'),
994 _T('-60, 60, 0, 60, -60, 0, 0, 60, -60, 0, -60, 60'),
995 _T('-20, -20, 0, 60, 60, 0, 0, -60, -60, 0, 20, 20')),
997 (_T('18, -18, 36, 210, 162, -36, 0, 0, 0, 0, -18, 18'), # PYCHOK indent
998 _T('-36, 36, -165, 45, 141, -21, 0, 0, 0, 0, 36, -36'),
999 _T('-122, -2, -27, -111, -75, 27, 62, 124, 124, 62, -64, 2'),
1000 _T('0, 0, 93, -93, -93, 93, 0, 0, 0, 0, 0, 0'),
1001 _T('120, -120, 147, -57, -129, 39, 0, 0, 0, 0, 66, -66'), # PYCHOK indent
1002 _T('135, 51, -9, -192, -180, 9, 31, 62, 62, 31, 51, -51'),
1003 _T0s12,
1004 _T('0, 0, -93, 93, 93, -93, 0, 0, 0, 0, 0, 0'),
1005 _T('-84, 84, 18, 12, -12, -18, 0, 0, 0, 0, -102, 102'),
1006 _T('-31, -31, 0, 93, 93, 0, -31, -62, -62, -31, 31, 31')))
1008 _BT = (_T('0, 0'), # bilinear 4-tuple [i, j] indices
1009 _T('1, 0'),
1010 _T('0, 1'),
1011 _T('1, 1'))
1013 _CM = (_T(' 0, -1'), # 10x12 cubic matrix [i, j] indices
1014 _T(' 1, -1'),
1015 _T('-1, 0'),
1016 _T(' 0, 0'), # _BT[0]
1017 _T(' 1, 0'), # _BT[1]
1018 _T(' 2, 0'),
1019 _T('-1, 1'),
1020 _T(' 0, 1'), # _BT[2]
1021 _T(' 1, 1'), # _BT[3]
1022 _T(' 2, 1'),
1023 _T(' 0, 2'),
1024 _T(' 1, 2'))
1026# _cropped = None
1027 _endian = '>H' # struct.unpack 1 ushort (big endian, unsigned short)
1028 _4endian = '>4H' # struct.unpack 4 ushorts
1029 _Rendian = NN # struct.unpack a row of ushorts
1030# _highest = (-8.4, 147.367, 85.839) if egm2008-1.pgm else (
1031# (-8.167, 147.25, 85.422) if egm96-5.pgm else
1032# (-4.5, 148.75, 81.33)) # egm84-15.pgm
1033 _iscipy = False
1034# _lowest = (4.7, 78.767, -106.911) if egm2008-1.pgm else (
1035# (4.667, 78.833, -107.043) if egm96-5.pgm else
1036# (4.75, 79.25, -107.34)) # egm84-15.pgm
1037 _mean = _F(-1.317) # from egm2008-1, -1.438 egm96-5, -0.855 egm84-15
1038 _nBytes = None # not applicable
1039 _nterms = len(_C3[0]) # columns length, number of rows
1040 _smooth = None # not applicable
1041 _stdev = _F(29.244) # from egm2008-1, 29.227 egm96-5, 29.183 egm84-15
1042 _u2B = _calcsize(_endian) # pixelsize_ in bytes
1043 _4u2B = _calcsize(_4endian) # 4 pixelsize_s in bytes
1044 _Ru2B = 0 # row of pixelsize_s in bytes
1045 _yx_hits = 0 # cache hits
1046 _yx_i = () # cached (y, x) indices
1047 _yx_t = () # cached 4- or 10-tuple for _ev2k resp. _ev3k
1049 def __init__(self, egm_pgm, crop=None, datum=_WGS84, kind=3, **name_smooth):
1050 '''New L{GeoidKarney} interpolator.
1052 @arg egm_pgm: An U{EGM geoid dataset<https://GeographicLib.SourceForge.io/
1053 C++/doc/geoid.html#geoidinst>} file name (C{egm*.pgm}), see
1054 note below.
1055 @kwarg crop: Optional box to limit geoid locations, a 4-tuple (C{south,
1056 west, north, east}), 2-tuple (C{(south, west), (north, east)})
1057 with 2 C{degrees90} lat- and C{degrees180} longitudes or as
1058 2-tuple (C{LatLonSW, LatLonNE}) of C{LatLon} instances.
1059 @kwarg datum: Optional grid datum (C{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1060 L{a_f2Tuple}), overriding C{WGS84}.
1061 @kwarg kind: Interpolation order (C{int}), 2 for C{bilinear} or 3 for C{cubic}.
1062 @kwarg name_smooth: Optional geoid C{B{name}=NN} (C{str}) and UNSUPPORTED
1063 keyword argument C{B{smooth}}, use C{B{smooth}=None} to ignore.
1065 @raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid B{C{crop}},
1066 B{C{kind}} or B{C{smooth}}.
1068 @raise TypeError: Invalid B{C{datum}}.
1070 @see: Class L{GeoidPGM} and function L{egmGeoidHeights}.
1072 @note: Geoid file B{C{egm_pgm}} remains open and I{must be closed} by calling
1073 method C{close} or by using C{with B{GeoidKarney}(...) as ...:} context.
1074 '''
1075 smooth, name = _xkwds_pop2(name_smooth, smooth=None)
1076 if smooth is not None:
1077 raise GeoidError(smooth=smooth, txt_not_=_supported_)
1079 if _isin(kind, 2):
1080 self._ev2d = self._ev2k # see ._ev_name
1081 elif not _isin(kind, 3):
1082 raise GeoidError(kind=kind)
1084 self._egm = g = self._open(egm_pgm, datum, kind, _name__(**name), None)
1085 self._pgm = p = _PGM(g, pgm=egm_pgm, itemsize=self.u2B, sizeB=self.sizeB)
1087 self._Rendian = self._4endian.replace(_4_, str(p.nlon))
1088 self._Ru2B = _calcsize(self._Rendian)
1090 self._knots = p.knots # grid knots
1091 self._lon_of = float(p.flon) # forward offset
1092 self._lon_og = float(p.glon) # reverse offset
1093 # set earth (lat, lon) limits (s, w, n, e)
1094 self._lat_lo, self._lon_lo, \
1095 self._lat_hi, self._lon_hi = self._swne(crop if crop else p.crop4)
1096 self._cropped = bool(crop)
1098 def _c0c3v(self, y, x):
1099 # get the common denominator, the 10x12 cubic matrix and
1100 # the 12 cubic v-coefficients around geoid index (y, x)
1101 p = self._pgm
1102 if 0 < x < (p.nlon - 2) and 0 < y < (p.nlat - 2):
1103 # read 4x4 ushorts, drop the 4 corners
1104 S = _os.SEEK_SET
1105 e = self._4endian
1106 g = self._egm
1107 n = self._4u2B
1108 R = self._Ru2B
1109 b = self._seek(y - 1, x - 1)
1110 v = _unpack(e, g.read(n))[1:3]
1111 b += R
1112 g.seek(b, S)
1113 v += _unpack(e, g.read(n))
1114 b += R
1115 g.seek(b, S)
1116 v += _unpack(e, g.read(n))
1117 b += R
1118 g.seek(b, S)
1119 v += _unpack(e, g.read(n))[1:3]
1120 j = 1
1122 else: # likely some wrapped y and/or x's
1123 v = self._raws(y, x, GeoidKarney._CM)
1124 j = 0 if y < 1 else (1 if y < (p.nlat - 2) else 2)
1126 return GeoidKarney._C0[j], GeoidKarney._C3[j], v
1128 @Property_RO
1129 def dtype(self):
1130 '''Get the geoid's grid data type (C{str}).
1131 '''
1132 return 'ushort'
1134 def _ev(self, lat, lon): # PYCHOK expected
1135 # interpolate the geoid height at grid (lat, lon)
1136 fy, fx = self._g2yx2(lat, lon)
1137 y, x = int(_floor(fy)), int(_floor(fx))
1138 fy -= y
1139 fx -= x
1140 H = self._ev2d(fy, fx, y, x) # PYCHOK ._ev3k or ._ev2k
1141 H *= self._pgm.Scale # H.fmul(self._pgm.Scale)
1142 H += self._pgm.Offset # H.fadd(self._pgm.Offset)
1143 return H.fsum() # float(H)
1145 def _ev2k(self, fy, fx, *yx):
1146 # compute the bilinear 4-tuple and interpolate raw H
1147 if self._yx_i == yx:
1148 self._yx_hits += 1
1149 else:
1150 y, x = self._yx_i = yx
1151 self._yx_t = self._raws(y, x, GeoidKarney._BT)
1152 t = self._yx_t
1153 v = _1_0, (-fx), fx
1154 H = _Dotf(v, t[0], t[0], t[1]).fmul(_1_0 - fy) # c = a * (1 - fy)
1155 H += _Dotf(v, t[2], t[2], t[3]).fmul(fy) # c += b * fy
1156 return H # Fsum
1158 def _ev3k(self, fy, fx, *yx):
1159 # compute the cubic 10-tuple and interpolate raw H
1160 if self._yx_i == yx:
1161 self._yx_hits += 1
1162 else:
1163 c0, c3, v = self._c0c3v(*yx)
1164 # assert len(c3) == self._nterms
1165 self._yx_t = tuple(_Dotf(v, *r3).fover(c0) for r3 in c3)
1166 self._yx_i = yx
1167 # GeographicLib/Geoid.cpp Geoid::height(lat, lon) ...
1168 # real h = t[0] + fx * (t[1] + fx * (t[3] + fx * t[6])) +
1169 # fy * (t[2] + fx * (t[4] + fx * t[7]) +
1170 # fy * (t[5] + fx * t[8] + fy * t[9]));
1171 t = self._yx_t
1172 v = _1_0, fx, fy
1173 H = _Dotf(v, t[5], t[8], t[9])
1174 H *= fy
1175 H += _Hornerf(fx, t[2], t[4], t[7])
1176 H *= fy
1177 H += _Hornerf(fx, t[0], t[1], t[3], t[6])
1178 return H # Fsum
1180 _ev2d = _ev3k # overriden for kind=2, see ._ev_name
1182 def _g2ll2(self, lat, lon):
1183 # convert grid (lat, lon) to earth (lat, lon), uncropped
1184 return lat, _lonE2lon(lon)
1186 def _g2yx2(self, lat, lon):
1187 # convert grid (lat, lon) to grid (y, x) indices
1188 p = self._pgm
1189 # note, slat = +90, rlat < 0 makes y >=0
1190 return ((lat - p.slat) * p.rlat), ((lon - p.wlon) * p.rlon)
1192 def _gyx2g2(self, y, x):
1193 # convert grid (y, x) indices to grid (lat, lon)
1194 p = self._pgm
1195 return (p.slat + p.dlat * y), (p.wlon + p.dlon * x)
1197 @Property_RO
1198 def _highest_ltd(self):
1199 '''(INTERNAL) Cache for C{.highest}.
1200 '''
1201 return self._LL3T(self._llh3minmax(True, -12, -4), name__=self.highest)
1203 def highest(self, LatLon=None, full=False): # PYCHOK full
1204 '''Return the location and largest height of this geoid.
1206 @kwarg LatLon: Optional class to return the location and height
1207 (C{LatLon}) or C{None}.
1208 @kwarg full: Search the full or limited latitude range (C{bool}).
1210 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,
1211 height)} otherwise a B{C{LatLon}} instance with the lat-,
1212 longitude and geoid height of the highest grid location.
1213 '''
1214 llh = self._highest if full or self.cropped else self._highest_ltd
1215 return self._llh3LL(llh, LatLon)
1217 def _lat2y2(self, lat2):
1218 # convert earth lat(s) to min and max grid y indices
1219 ys, m = [], self._pgm.nlat - 1
1220 for lat in lat2:
1221 y, _ = self._g2yx2(*self._ll2g2(lat, 0))
1222 ys.append(max(min(int(y), m), 0))
1223 return min(ys), max(ys) + 1
1225 def _ll2g2(self, lat, lon):
1226 # convert earth (lat, lon) to grid (lat, lon), uncropped
1227 return lat, _lon2lonE(lon)
1229 def _llh3minmax(self, highest, *lat2):
1230 # find highest or lowest, takes 10+ secs for egm2008-1.pgm geoid
1231 # (Python 2.7.16, macOS 10.13.6 High Sierra, iMac 3 GHz Core i3)
1232 if highest:
1233 def _mt(r, h):
1234 m = max(r)
1235 return m, (m > h)
1237 else: # lowest
1238 def _mt(r, h): # PYCHOK redef
1239 m = min(r)
1240 return m, (m < h)
1242 y = x = 0
1243 h = self._raw(y, x)
1244 for j, r in self._raw2(*lat2):
1245 m, t = _mt(r, h)
1246 if t:
1247 h, y, x = m, j, r.index(m)
1248 h *= self._pgm.Scale
1249 h += self._pgm.Offset
1250 return self._g2ll2(*self._gyx2g2(y, x)) + (h,)
1252 @Property_RO
1253 def _lowest_ltd(self):
1254 '''(INTERNAL) Cache for C{.lowest}.
1255 '''
1256 return self._LL3T(self._llh3minmax(False, 0, 8), name__=self.lowest)
1258 def lowest(self, LatLon=None, full=False): # PYCHOK full
1259 '''Return the location and lowest height of this geoid.
1261 @kwarg LatLon: Optional class to return the location and height
1262 (C{LatLon}) or C{None}.
1263 @kwarg full: Search the full or limited latitude range (C{bool}).
1265 @return: If C{B{LatLon} is None}, a L{LatLon3Tuple}C{(lat, lon,
1266 height)} otherwise a B{C{LatLon}} instance with the lat-,
1267 longitude and geoid height of the lowest grid location.
1268 '''
1269 llh = self._lowest if full or self.cropped else self._lowest_ltd
1270 return self._llh3LL(llh, LatLon)
1272 def _raw(self, y, x):
1273 # get the ushort geoid height at geoid index (y, x),
1274 # like GeographicLib/Geoid.hpp real rawval(is, iy)
1275 p = self._pgm
1276 if x < 0:
1277 x += p.nlon
1278 elif x >= p.nlon:
1279 x -= p.nlon
1280 h = p.nlon // 2
1281 if y < 0:
1282 y = -y
1283 elif y >= p.nlat:
1284 y = (p.nlat - 1) * 2 - y
1285 else:
1286 h = 0
1287 x += h if x < h else -h
1288 self._seek(y, x)
1289 h = _unpack(self._endian, self._egm.read(self._u2B))
1290 return h[0]
1292 def _raws(self, y, x, ijs):
1293 # get bilinear 4-tuple or 10x12 cubic matrix
1294 return tuple(self._raw(y + j, x + i) for i, j in ijs)
1296 def _raw2(self, *lat2):
1297 # yield a 2-tuple (y, ushorts) for each row or for
1298 # the rows between two (or more) earth lat values
1299 p = self._pgm
1300 g = self._egm
1301 e = self._Rendian
1302 n = self._Ru2B
1303 # min(lat2) <= lat <= max(lat2) or 0 <= y < p.nlat
1304 s, t = self._lat2y2(lat2) if lat2 else (0, p.nlat)
1305 self._seek(s, 0) # to start of row s
1306 for y in range(s, t):
1307 yield y, _unpack(e, g.read(n))
1309 def _seek(self, y, x):
1310 # position geoid to grid index (y, x)
1311 p, g = self._pgm, self._egm
1312 if g:
1313 b = p.skip + (y * p.nlon + x) * self._u2B
1314 g.seek(b, _os.SEEK_SET)
1315 return b # position
1316 raise GeoidError('closed file', txt=repr(p.egm)) # IOError
1319class GeoidPGM(_GeoidBase):
1320 '''Geoid height interpolator for I{Karney}'s U{GeographicLib Earth
1321 Gravitational Model (EGM)<https://GeographicLib.SourceForge.io/C++/doc/geoid.html>}
1322 geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>}
1323 datasets but based on C{SciPy} U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/
1324 reference/generated/scipy.interpolate.RectBivariateSpline.html>}, U{bisplrep/-ev
1325 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}
1326 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.
1327 interp2d.html>} interpolation.
1329 Use any of the U{egm84-, egm96- or egm2008-*.pgm <https://GeographicLib.SourceForge.io/
1330 C++/doc/geoid.html#geoidinst>} datasets. However, unless cropped, an entire C{egm*.pgm}
1331 dataset is loaded into the C{SciPy} interpolator and converted from 2-byte C{int} to
1332 8-byte C{dtype float64}. Therefore, internal memory usage is 4x the U{egm*.pgm
1333 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>} file size and may
1334 exceed the available memory, especially with 32-bit Python, see properties C{.nBytes}
1335 and C{.sizeB}.
1336 '''
1337 _cropped = False
1338 _endian = '>u2'
1340 def __init__(self, egm_pgm, crop=None, datum=_WGS84, kind=3, smooth=0, **name):
1341 '''New L{GeoidPGM} interpolator.
1343 @arg egm_pgm: An U{EGM geoid dataset<https://GeographicLib.SourceForge.io/
1344 C++/doc/geoid.html#geoidinst>} file name (C{egm*.pgm}).
1345 @kwarg crop: Optional box to crop B{C{egm_pgm}}, a 4-tuple (C{south, west,
1346 north, east}) or 2-tuple (C{(south, west), (north, east)}),
1347 in C{degrees90} lat- and C{degrees180} longitudes or a 2-tuple
1348 (C{LatLonSW, LatLonNE}) of C{LatLon} instances.
1349 @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
1350 L{a_f2Tuple}), overriding C{WGS84}.
1351 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline
1352 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.
1353 RectBivariateSpline.html>} or -1, -3 or -5 for U{bisplrep/-ev<https://
1354 docs.SciPy.org/doc/scipy/reference/generated/scipy.interpolate.bisplrep.html>}
1355 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
1356 interpolate.interp2d.html>} C{linear}, C{cubic} respectively C{quintic},
1357 see note for more details.
1358 @kwarg smooth: Smoothing factor for C{B{kind}=1..5} only (C{int}).
1359 @kwarg name: Optional geoid C{B{name}=NN} (C{str}).
1361 @raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid B{C{crop}}, B{C{kind}}
1362 or B{C{smooth}}.
1364 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.
1366 @raise LenError: EGM dataset B{C{egm_pgm}} axis mismatch.
1368 @raise SciPyError: A C{scipy} issue.
1370 @raise SciPyWarning: A C{scipy} warning as exception.
1372 @raise TypeError: Invalid B{C{datum}} or unexpected argument.
1374 @note: Specify C{B{kind}=-1, -3 or -5} to use C{scipy.interpolate.interp2d}
1375 before or C{scipy.interpolate.bisplrep/-ev} since C{Scipy} version 1.14.
1377 @note: The U{GeographicLib egm*.pgm<https://GeographicLib.SourceForge.io/C++/doc/
1378 geoid.html#geoidinst>} file sizes are based on a 2-byte C{int} height
1379 converted to 8-byte C{dtype float64} for C{scipy} interpolators. Therefore,
1380 internal memory usage is 4 times the C{egm*.pgm} file size and may exceed
1381 the available memory, especially with 32-bit Python. To reduce memory
1382 usage, set keyword argument B{C{crop}} to the region of interest. For
1383 example C{B{crop}=(20, -125, 50, -65)} covers the U{conterminous US
1384 <https://Geodesy.NOAA.gov/GEOID/GEOID12B/maps/GEOID12B_CONUS_grids.png>}
1385 (CONUS), less than 3% of the entire C{egm2008-1.pgm} dataset.
1387 @see: Class L{GeoidKarney} and function L{egmGeoidHeights}.
1388 '''
1389 np = self.numpy
1390 self._u2B = np.dtype(self.endian).itemsize
1392 g = self._open(egm_pgm, datum, kind, _name__(**name), smooth)
1393 self._pgm = p = _PGM(g, pgm=egm_pgm, itemsize=self.u2B, sizeB=self.sizeB)
1394 if crop:
1395 g = p._cropped(g, abs(kind) + 1, *self._swne(crop))
1396 if _MODS.internals._version2(np.__version__) < (1, 9):
1397 g = open(g.name, _rb_) # reopen tempfile for numpy 1.8.0-
1398 self._cropped = True
1399 try:
1400 # U{numpy dtype formats are different from Python struct formats
1401 # <https://docs.SciPy.org/doc/numpy-1.15.0/reference/arrays.dtypes.html>}
1402 # read all heights, skipping the PGM header lines, converted to float
1403 hs = self._load(g, self.endian, p.knots, p.skip).reshape(p.nlat, p.nlon) * p.Scale
1404 if p.Offset: # offset
1405 hs = p.Offset + hs
1406 if p.dlat < 0: # flip the rows
1407 hs = np.flipud(hs)
1408 _GeoidBase.__init__(self, hs, p)
1409 except Exception as x:
1410 raise _SciPyIssue(x, _in_, repr(egm_pgm))
1411 finally:
1412 g.close()
1414 def _g2ll2(self, lat, lon):
1415 # convert grid (lat, lon) to earth (lat, lon), un-/cropped
1416 if self._cropped:
1417 lon -= self._lon_of
1418 else:
1419 lon = _lonE2lon(lon)
1420 return lat, lon
1422 def _ll2g2(self, lat, lon):
1423 # convert earth (lat, lon) to grid (lat, lon), un-/cropped
1424 if self._cropped:
1425 lon += self._lon_of
1426 else:
1427 lon = _lon2lonE(lon)
1428 return lat, lon
1430 if _FOR_DOCS:
1431 __call__ = _GeoidBase.__call__
1432 height = _GeoidBase.height
1435class _Gpars(_Named):
1436 '''(INTERNAL) Basic geoid parameters.
1437 '''
1438 # interpolator parameters
1439 dlat = 0 # +/- latitude resolution in C{degrees}
1440 dlon = 0 # longitude resolution in C{degrees}
1441 nlat = 1 # number of latitude knots (C{int})
1442 nlon = 0 # number of longitude knots (C{int})
1443 rlat = 0 # +/- latitude resolution in C{float}, 1 / .dlat
1444 rlon = 0 # longitude resolution in C{float}, 1 / .dlon
1445 slat = 0 # nothern- or southern most latitude (C{degrees90})
1446 wlon = 0 # western-most longitude in Eastern lon (C{degrees360})
1448 flon = 0 # forward, earth to grid longitude offset
1449 glon = 0 # reverse, grid to earth longitude offset
1451 knots = 0 # number of knots, nlat * nlon (C{int})
1452 skip = 0 # header bytes to skip (C{int})
1454 def __repr__(self):
1455 t = _COMMASPACE_.join(pairs((a, getattr(self, a)) for
1456 a in dir(self.__class__)
1457 if a[:1].isupper()))
1458 return _COLONSPACE_(self, t)
1460 def __str__(self):
1461 return Fmt.PAREN(self.classname, repr(self.name))
1464class _PGM(_Gpars):
1465 '''(INTERNAL) Parse an C{egm*.pgm} geoid dataset file.
1467 # Geoid file in PGM format for the GeographicLib::Geoid class
1468 # Description WGS84 EGM96, 5-minute grid
1469 # URL https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/egm96/egm96.html
1470 # DateTime 2009-08-29 18:45:03
1471 # MaxBilinearError 0.140
1472 # RMSBilinearError 0.005
1473 # MaxCubicError 0.003
1474 # RMSCubicError 0.001
1475 # Offset -108
1476 # Scale 0.003
1477 # Origin 90N 0E
1478 # AREA_OR_POINT Point
1479 # Vertical_Datum WGS84
1480 <width> <height>
1481 <pixel>
1482 ...
1483 '''
1484 crop4 = () # 4-tuple (C{south, west, north, east}).
1485 egm = None
1486 glon = 180 # reverse offset, uncropped
1487# pgm = NN # name
1488 sizeB = 0
1489 u2B = 2 # item size of grid height (C{int}).
1491 @staticmethod
1492 def _llstr2floats(latlon):
1493 # llstr to (lat, lon) floats
1494 lat, lon = latlon.split()
1495 return _MODS.dms.parseDMS2(lat, lon)
1497 # PGM file attributes, CamelCase but not .istitle()
1498 AREA_OR_POINT = str
1499 DateTime = str
1500 Description = str # 'WGS84 EGM96, 5-minute grid'
1501 Geoid = str # 'file in PGM format for the GeographicLib::Geoid class'
1502 MaxBilinearError = float
1503 MaxCubicError = float
1504 Offset = float
1505 Origin = _llstr2floats
1506 Pixel = 0
1507 RMSBilinearError = float
1508 RMSCubicError = float
1509 Scale = float
1510 URL = str # 'https://Earth-Info.NGA.mil/GandG/wgs84/...'
1511 Vertical_Datum = str
1513 def __init__(self, g, pgm=NN, itemsize=0, sizeB=0): # MCCABE 22
1514 '''(INTERNAL) New C{_PGM} parsed C{egm*.pgm} geoid dataset.
1515 '''
1516 self.name = pgm # geoid file name
1517 if itemsize:
1518 self._u2B = itemsize
1519 if sizeB:
1520 self.sizeB = sizeB
1522 t = g.readline() # make sure newline == '\n'
1523 if t != b'P5\n' and t.strip() != b'P5':
1524 raise self._Errorf(_format_, _header_, t)
1526 while True: # read all # Attr ... lines,
1527 try: # ignore empty ones or comments
1528 t = g.readline().strip()
1529 if t.startswith(_bHASH_):
1530 t = t.lstrip(_bHASH_).lstrip()
1531 a, v = map(_ub2str, t.split(None, 1))
1532 f = getattr(_PGM, a, None)
1533 if callable(f) and a[:1].isupper():
1534 setattr(self, a, f(v))
1535 elif t:
1536 break
1537 except (TypeError, ValueError):
1538 raise self._Errorf(_format_, 'Attr', t)
1539 else: # should never get here
1540 raise self._Errorf(_format_, _header_, g.tell())
1542 try: # must be (even) width and (odd) height
1543 nlon, nlat = map(int, t.split())
1544 if nlon < 2 or nlon > (360 * 60) or isodd(nlon) or \
1545 nlat < 2 or nlat > (181 * 60) or not isodd(nlat):
1546 raise ValueError
1547 except (TypeError, ValueError):
1548 raise self._Errorf(_format_, _SPACE_(_width_, _height_), t)
1550 try: # must be 16 bit pixel height
1551 t = g.readline().strip()
1552 self.Pixel = int(t)
1553 if not 255 < self.Pixel < 65536: # >u2 or >H only
1554 raise ValueError
1555 except (TypeError, ValueError):
1556 raise self._Errorf(_format_, 'pixel', t)
1558 for a in dir(_PGM): # set undefined # Attr ... to None
1559 if a[:1].isupper() and callable(getattr(self, a)):
1560 setattr(self, a, None)
1562 if self.Origin is None:
1563 raise self._Errorf(_format_, 'Origin', self.Origin)
1564 if self.Offset is None or self.Offset > 0:
1565 raise self._Errorf(_format_, 'Offset', self.Offset)
1566 if self.Scale is None or self.Scale < EPS:
1567 raise self._Errorf(_format_, 'Scale', self.Scale)
1569 self.skip = g.tell()
1570 self.knots = nlat * nlon
1572 self.nlat, self.nlon = nlat, nlon
1573 self.slat, self.wlon = self.Origin
1574 # note, negative .dlat and .rlat since rows
1575 # are from .slat 90N down in decreasing lat
1576 self.dlat, self.dlon = (_180_0 / (1 - nlat)), (_360_0 / nlon)
1577 self.rlat, self.rlon = ((1 - nlat) / _180_0), (nlon / _360_0)
1579 # grid corners in earth (lat, lon), .slat = 90, .dlat < 0
1580 n = float(self.slat)
1581 s = n + self.dlat * (nlat - 1)
1582 w = self.wlon - self.glon
1583 e = w + self.dlon * nlon
1584 self.crop4 = s, w, n, e
1586 n = self.sizeB - self.skip
1587 if n > 0 and n != (self.knots * self.u2B):
1588 raise self._Errorf('%s(%s x %s != %s)', _assert_, nlat, nlon, n)
1590 def _cropped(self, g, k1, south, west, north, east): # MCCABE 15
1591 '''Crop the geoid to (south, west, north, east) box.
1592 '''
1593 # flon offset for both west and east
1594 f = 360 if west < 0 else 0
1595 # earth (lat, lon) to grid indices (y, x),
1596 # note y is decreasing, i.e. n < s
1597 s, w = self._lle2yx2(south, west, f)
1598 n, e = self._lle2yx2(north, east, f)
1599 s += 1 # s > n
1600 e += 1 # e > w
1602 hi, wi = self.nlat, self.nlon
1603 # handle special cases
1604 if (s - n) > hi:
1605 n, s = 0, hi # entire lat range
1606 if (e - w) > wi:
1607 w, e, f = 0, wi, 180 # entire lon range
1608 if s == hi and w == n == 0 and e == wi:
1609 return g # use entire geoid as-is
1611 if (e - w) < k1 or (s - n) < (k1 + 1):
1612 raise self._Errorf(_format_, 'swne', (north - south, east - west))
1614 if e > wi > w: # wrap around
1615 # read w..wi and 0..e
1616 r, p = (wi - w), (e - wi)
1617 elif e > w:
1618 r, p = (e - w), 0
1619 else:
1620 raise self._Errorf('%s(%s < %s)', _assert_, w, e)
1622 # convert to bytes
1623 r *= self.u2B
1624 p *= self.u2B
1625 q = wi * self.u2B # stride
1626 # number of rows and cols to skip from
1627 # the original (.slat, .wlon) origin
1628 z = self.skip + (n * wi + w) * self.u2B
1629 # sanity check
1630 if r < 2 or p < 0 or q < 2 or z < self.skip \
1631 or z > self.sizeB:
1632 raise self._Errorf(_format_, _assert_, (r, p, q, z))
1634 # can't use _BytesIO since numpy
1635 # needs .fileno attr in .fromfile
1636 t, c = 0, self._tmpfile()
1637 # reading (s - n) rows, forward
1638 for y in range(n, s): # PYCHOK y unused
1639 g.seek(z, _os.SEEK_SET)
1640 # Python 2 tmpfile.write returns None
1641 t += c.write(g.read(r)) or r
1642 if p: # wrap around to start of row
1643 g.seek(-q, _os.SEEK_CUR)
1644 # assert(g.tell() == (z - w * self.u2B))
1645 # Python 2 tmpfile.write returns None
1646 t += c.write(g.read(p)) or p
1647 z += q
1648 c.flush()
1649 g.close()
1651 s -= n # nlat
1652 e -= w # nlon
1653 k = s * e # knots
1654 z = k * self.u2B
1655 if t != z:
1656 raise self._Errorf('%s(%s != %s) %s', _assert_, t, z, self)
1658 # update the _Gpars accordingly, note attributes
1659 # .dlat, .dlon, .rlat and .rlon remain unchanged
1660 self.slat += n * self.dlat
1661 self.wlon += w * self.dlon
1662 self.nlat = s
1663 self.nlon = e
1664 self.flon = self.glon = f
1666 self.crop4 = south, west, north, east
1667 self.knots = k
1668 self.skip = 0 # no header lines in c
1670 c.seek(0, _os.SEEK_SET)
1671 # c = open(c.name, _rb_) # reopen for numpy 1.8.0-
1672 return c
1674 def _Errorf(self, fmt, *args): # PYCHOK no cover
1675 t = fmt % args
1676 e = self.pgm or NN
1677 if e:
1678 t = _SPACE_(t, _in_, repr(e))
1679 return PGMError(t)
1681 def _lle2yx2(self, lat, lon, flon):
1682 # earth (lat, lon) to grid indices (y, x)
1683 # with .dlat decreasing from 90N .slat
1684 lat -= self.slat
1685 lon += flon - self.wlon
1686 return (min(self.nlat - 1, max(0, int(lat * self.rlat))),
1687 max(0, int(lon * self.rlon)))
1689 def _tmpfile(self):
1690 # create a tmpfile to hold the cropped geoid grid
1691 try:
1692 from tempfile import NamedTemporaryFile as tmpfile
1693 except ImportError: # Python 2.7.16-
1694 from _os import tmpfile # PYCHOK from
1695 t = _os.path.basename(self.pgm)
1696 t = _os.path.splitext(t)[0]
1697 f = tmpfile(mode='w+b', prefix=t or 'egm')
1698 f.seek(0, _os.SEEK_SET) # force overwrite
1699 return f
1701 @Property_RO
1702 def pgm(self):
1703 '''Get the geoid file name (C{str}).
1704 '''
1705 return self.name
1708class PGMError(GeoidError):
1709 '''An issue while parsing or cropping an C{egm*.pgm} geoid dataset.
1710 '''
1711 pass
1714def egmGeoidHeights(GeoidHeights_dat):
1715 '''Generate geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/
1716 C++/doc/geoid.html#geoidinst>} height tests from U{GeoidHeights.dat
1717 <https://SourceForge.net/projects/geographiclib/files/testdata/>}
1718 U{Test data for Geoids<https://GeographicLib.SourceForge.io/C++/doc/
1719 geoid.html#testgeoid>}.
1721 @arg GeoidHeights_dat: The un-gz-ed C{GeoidHeights.dat} file
1722 (C{str} or C{file} handle).
1724 @return: For each test, yield a L{GeoidHeight5Tuple}C{(lat, lon,
1725 egm84, egm96, egm2008)}.
1727 @raise GeoidError: Invalid B{C{GeoidHeights_dat}}.
1729 @note: Function L{egmGeoidHeights} is used to test the geoids
1730 L{GeoidKarney} and L{GeoidPGM}, see PyGeodesy module
1731 C{test/testGeoids.py}.
1732 '''
1733 dat = GeoidHeights_dat
1734 if isinstance(dat, bytes):
1735 dat = _BytesIO(dat)
1737 try:
1738 dat.seek(0, _os.SEEK_SET) # reset
1739 except AttributeError as x:
1740 raise GeoidError(GeoidHeights_dat=type(dat), cause=x)
1742 for t in dat.readlines():
1743 t = t.strip()
1744 if t and not t.startswith(_bHASH_):
1745 lat, lon, egm84, egm96, egm2008 = map(float, t.split())
1746 lon = _lonE2lon(lon) # Eastern to earth lon
1747 yield GeoidHeight5Tuple(lat, lon, egm84, egm96, egm2008)
1750def _lonE2lon(lon):
1751 '''(INTERNAL) East to earth longitude.
1752 '''
1753 while lon > _180_0:
1754 lon -= _360_0
1755 return lon
1758def _lon2lonE(lon):
1759 '''(INTERNAL) Earth to East longitude.
1760 '''
1761 while lon < 0:
1762 lon += _360_0
1763 return lon
1766__all__ += _ALL_DOCS(_GeoidBase)
1768if __name__ == _DMAIN_: # MCCABE 14
1770 from pygeodesy.internals import printf, _secs2str, _versions, _sys
1771 from time import time
1773 _crop = {}
1774 _GeoidEGM = GeoidKarney
1775 _kind = 3
1777 geoids = _sys.argv[1:]
1778 while geoids:
1779 G = geoids.pop(0)
1780 g = G.lower()
1782 if '-crop'.startswith(g):
1783 _crop = dict(crp=(20, -125, 50, -65)) # CONUS
1785 elif '-egm96'.startswith(g):
1786 _GeoidEGM = GeoidEGM96
1788 elif '-karney'.startswith(g):
1789 _GeoidEGM = GeoidKarney
1791 elif '-kind'.startswith(g):
1792 _kind = int(geoids.pop(0))
1794 elif '-pgm'.startswith(g):
1795 _GeoidEGM = GeoidPGM
1797 elif _isin(g[-4:], '.pgm', '.grd'):
1798 g = _GeoidEGM(G, kind=_kind, **_crop)
1799 t = time()
1800 _ = g.highest()
1801 t = _secs2str(time() - t)
1802 printf('%s: %s (%s)', g.toStr(), t, _versions(), nl=1, nt=1)
1803 t = g.pgm
1804 if t:
1805 printf(repr(t), nt=1)
1806 # <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>:
1807 # The height of the EGM96 geoid at Timbuktu
1808 # echo 16:46:33N 3:00:34W | GeoidEval
1809 # => 28.7068 -0.02e-6 -1.73e-6
1810 # The 1st number is the height of the geoid, the 2nd and
1811 # 3rd are its slopes in northerly and easterly direction
1812 t = 'Timbuktu %s' % (g,)
1813 k = {'egm84-15.pgm': '31.2979',
1814 'egm96-5.pgm': '28.7067',
1815 'egm2008-1.pgm': '28.7880'}.get(g.name.lower(), '28.7880')
1816 ll = _MODS.dms.parseDMS2('16:46:33N', '3:00:34W', sep=':')
1817 for ll in (ll, (16.776, -3.009),):
1818 try:
1819 h, ll = g.height(*ll), fstr(ll, prec=6)
1820 printf('%s.height(%s): %.4F vs %s', t, ll, h, k)
1821 except (GeoidError, RangeError) as x:
1822 printf(_COLONSPACE_(t, str(x)))
1824 elif _isin(g[-4:], '.bin'):
1825 g = GeoidG2012B(G, kind=_kind)
1826 printf(g.toStr())
1828 else:
1829 raise GeoidError(grid=repr(G))
1831_I = int # PYCHOK unused _I
1832del _intCs, _T, _T0s12 # trash ints cache and map
1835# <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>
1836# _lowerleft = -90, -179, -30.1500 # egm2008-1.pgm
1837# _lowerleft = -90, -179, -29.5350 # egm96-5.pgm
1838# _lowerleft = -90, -179, -29.7120 # egm84-15.pgm
1840# _center = 0, 0, 17.2260 # egm2008-1.pgm
1841# _center = 0, 0, 17.1630 # egm96-5.pgm
1842# _center = 0, 0, 18.3296 # egm84-15.pgm
1844# _upperright = 90, 180, 14.8980 # egm2008-1.pgm
1845# _upperright = 90, 180, 13.6050 # egm96-5.pgm
1846# _upperright = 90, 180, 13.0980 # egm84-15.pgm
1849# % python3.12 -m pygeodesy.geoids -egm96 ../testGeoids/WW15MGH.GRD
1850#
1851# GeoidEGM96('WW15MGH.GRD'): lowerleft(-90.0, -180.0, -29.534), upperright(90.0, 180.25, 13.606), center(0.0, 0.125, 17.125), highest(-8.25, -32.75, 85.391), lowest(4.75, -101.25, -106.991): 1.267 ms (pygeodesy 24.12.24 Python 3.12.7 64bit arm64 macOS 14.6.1)
1852#
1853# Timbuktu GeoidEGM96('WW15MGH.GRD').height(16.775833, -3.009444): 28.7073 vs 28.7880
1854# Timbuktu GeoidEGM96('WW15MGH.GRD').height(16.776, -3.009): 28.7072 vs 28.7880
1857# % python3.12 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm
1858#
1859# GeoidKarney('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, 147.367, 85.839), lowest(4.7, 78.767, -106.911): 204.334 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)
1860#
1861# _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'
1862#
1863# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880
1864# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880
1865#
1866# GeoidKarney('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, 148.75, 81.33), lowest(4.75, 79.25, -107.34): 1.007 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)
1867#
1868# _PGM('../testGeoids/egm84-15.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 15-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.413, MaxCubicError=0.02, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84'
1869#
1870# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979
1871# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979
1872#
1873# GeoidKarney('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, 0.0, 17.163), highest(-8.167, 147.25, 85.422), lowest(4.667, 78.833, -107.043): 8.509 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)
1874#
1875# _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'
1876#
1877# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067
1878# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067
1881# % python3.8 -m pygeodesy.geoids -egm96 ../testGeoids/WW15MGH.GRD
1882#
1883# GeoidEGM96('WW15MGH.GRD'): lowerleft(-90.0, -180.0, -29.534), upperright(90.0, 180.25, 13.606), center(0.0, 0.125, 17.125), highest(-8.25, -32.75, 85.391), lowest(4.75, -101.25, -106.991): 1.267 ms (pygeodesy 24.12.24 Python 3.8.10 64bit arm64_x86_64 macOS 10.16)
1884#
1885# Timbuktu GeoidEGM96('WW15MGH.GRD').height(16.775833, -3.009444): 28.7073 vs 28.7880
1886# Timbuktu GeoidEGM96('WW15MGH.GRD').height(16.776, -3.009): 28.7072 vs 28.7880
1889# % python3.8 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm
1890#
1891# GeoidKarney('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, 147.367, 85.839), lowest(4.7, 78.767, -106.911): 353.050 ms (pygeodesy 24.8.24 Python 3.8.10 64bit arm64_x86_64 macOS 10.16)
1892#
1893# _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'
1894#
1895# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880
1896# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880
1897#
1898# GeoidKarney('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, 148.75, 81.33), lowest(4.75, 79.25, -107.34): 1.727 ms (pygeodesy 24.8.24 Python 3.8.10 64bit arm64_x86_64 macOS 10.16)
1899#
1900# _PGM('../testGeoids/egm84-15.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 15-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.413, MaxCubicError=0.02, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84'
1901#
1902# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979
1903# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979
1904#
1905# GeoidKarney('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, 0.0, 17.163), highest(-8.167, 147.25, 85.422), lowest(4.667, 78.833, -107.043): 14.807 ms (pygeodesy 24.8.24 Python 3.8.10 64bit arm64_x86_64 macOS 10.16)
1906#
1907# _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'
1908#
1909# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067
1910# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067
1913# % python2 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm
1914#
1915# GeoidKarney('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, 147.367, 85.839), lowest(4.7, 78.767, -106.911): 283.362 ms (pygeodesy 24.8.24 Python 2.7.18 64bit arm64_x86_64 macOS 10.16)
1916#
1917# _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'
1918#
1919# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880
1920# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880
1921#
1922# GeoidKarney('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, 148.75, 81.33), lowest(4.75, 79.25, -107.34): 1.378 ms (pygeodesy 24.8.24 Python 2.7.18 64bit arm64_x86_64 macOS 10.16)
1923#
1924# _PGM('../testGeoids/egm84-15.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 15-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.413, MaxCubicError=0.02, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84'
1925#
1926# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979
1927# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979
1928#
1929# GeoidKarney('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, 0.0, 17.163), highest(-8.167, 147.25, 85.422), lowest(4.667, 78.833, -107.043): 11.612 ms (pygeodesy 24.8.24 Python 2.7.18 64bit arm64_x86_64 macOS 10.16)
1930#
1931# _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'
1932#
1933# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067
1934# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067
1937# % python3.12 -m pygeodesy.geoids -PGM ../testGeoids/egm*.pgm
1938#
1939# GeoidPGM('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, -32.633, 85.839), lowest(4.683, -101.25, -106.911): 543.148 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)
1940#
1941# _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'
1942#
1943# Timbuktu GeoidPGM('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880
1944# Timbuktu GeoidPGM('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880
1945#
1946# GeoidPGM('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, -31.25, 81.33), lowest(4.75, -100.75, -107.34): 1.762 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)
1947#
1948# _PGM('../testGeoids/egm84-15.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 15-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.413, MaxCubicError=0.02, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84'
1949#
1950# Timbuktu GeoidPGM('egm84-15.pgm').height(16.775833, -3.009444): 31.2979 vs 31.2979
1951# Timbuktu GeoidPGM('egm84-15.pgm').height(16.776, -3.009): 31.2975 vs 31.2979
1952#
1953# GeoidPGM('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, -0.0, 17.179), highest(-8.167, -32.75, 85.422), lowest(4.667, -101.167, -107.043): 12.594 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)
1954#
1955# _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'
1956#
1957# Timbuktu GeoidPGM('egm96-5.pgm').height(16.775833, -3.009444): 28.7065 vs 28.7067
1958# Timbuktu GeoidPGM('egm96-5.pgm').height(16.776, -3.009): 28.7064 vs 28.7067
1960# **) MIT License
1961#
1962# Copyright (C) 2016-2026 -- mrJean1 at Gmail -- All Rights Reserved.
1963#
1964# Permission is hereby granted, free of charge, to any person obtaining a
1965# copy of this software and associated documentation files (the "Software"),
1966# to deal in the Software without restriction, including without limitation
1967# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1968# and/or sell copies of the Software, and to permit persons to whom the
1969# Software is furnished to do so, subject to the following conditions:
1970#
1971# The above copyright notice and this permission notice shall be included
1972# in all copies or substantial portions of the Software.
1973#
1974# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1975# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1976# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1977# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1978# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1979# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1980# OTHER DEALINGS IN THE SOFTWARE.