Coverage for pygeodesy/geoids.py: 96%
664 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-06-07 08:37 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-06-07 08:37 -0400
2# -*- coding: utf-8 -*-
4u'''Geoid models and geoid height interpolations.
6Classes L{GeoidG2012B}, L{GeoidKarney} and L{GeoidPGM} to interpolate the
7height of various U{geoid<https://WikiPedia.org/wiki/Geoid>}s at C{LatLon}
8locations or separate lat-/longitudes using different interpolation methods
9and 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.
13The L{GeoidG2012B} and L{GeoidPGM} interpolators both depend on U{scipy
14<https://SciPy.org>} and U{numpy<https://PyPI.org/project/numpy>} and
15require those packages to be installed.
17In addition, each geoid interpolator needs C{grid knots} (down)loaded from
18a C{geoid} model file, I{specific to the interpolator}, more details below
19and in the documentation of the interpolator class. For each interpolator,
20there are several interpolation choices, like I{linear}, I{cubic}, etc.
22Typical usage
23=============
251. Choose one of the interpolator classes L{GeoidG2012B}, L{GeoidKarney}
26or L{GeoidPGM} and download a C{geoid} model file, containing locations with
27known heights also referred to as the C{grid knots}. See the documentation
28of interpolator class for references to available C{grid} models.
30C{>>> from pygeodesy import GeoidG2012B # or -Karney or -PGM as GeoidXyz}
322. Instantiate an interpolator with the C{geoid} model file and use keyword
33arguments to select different interpolation options
35C{>>> ginterpolator = GeoidXyz(geoid_model_file, **options)}
373. Get the interpolated geoid height of other C{LatLon} location(s) with
39C{>>> ll = LatLon(1, 2, ...)}
40C{>>> h = ginterpolator(ll)}
42or
44C{>>> h0, h1, h2, ... = ginterpolator(ll0, ll1, ll2, ...)}
46or a list, tuple, generator, etc. of C{LatLon}s
48C{>>> hs = ginterpolator(lls)}
504. For separate lat- and longitudes invoke the C{.height} method as
52C{>>> h = ginterpolator.height(lat, lon)}
54or as 2 lists, 2 tuples, etc.
56C{>>> hs = ginterpolator.height(lats, lons)}
585. An example is in U{issue #64<https://GitHub.com/mrJean1/PyGeodesy/issues/64>},
59courtesy of SBFRF.
61@note: Classes L{GeoidG2012B} and L{GeoidPGM} require both U{numpy
62 <https://PyPI.org/project/numpy>} and U{scipy<https://PyPI.org/project/scipy>}
63 to be installed.
65@note: Errors from C{scipy} are raised as L{SciPyError}s. Warnings issued by
66 C{scipy} can be thrown as L{SciPyWarning} exceptions, provided Python
67 C{warnings} are filtered accordingly, see L{SciPyWarning}.
69@see: I{Karney}'s U{GeographicLib<https://GeographicLib.SourceForge.io/html/index.html>},
70 U{Geoid height<https://GeographicLib.SourceForge.io/html/geoid.html>} and U{Installing
71 the Geoid datasets<https://GeographicLib.SourceForge.io/html/geoid.html#geoidinst>},
72 U{SciPy<https://docs.SciPy.org/doc/scipy/reference/interpolate.html>} interpolation
73 U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
74 interpolate.RectBivariateSpline.html>} and U{interp2d<https://docs.SciPy.org/doc/scipy/
75 reference/generated/scipy.interpolate.interp2d.html>}, functions L{elevations.elevation2}
76 and L{elevations.geoidHeight2}, U{I{Ellispoid vs Orthometric Elevations}<https://
77 www.YouTube.com/watch?v=dX6a6kCk3Po>} and U{I{Pitfalls Related to Ellipsoid Height
78 and Height Above Mean Sea Level (AMSL)}<https://Wiki.ROS.org/mavros#mavros.2FPlugins.
79 Avoiding_Pitfalls_Related_to_Ellipsoid_Height_and_Height_Above_Mean_Sea_Level>}.
80'''
81# make sure int/int division yields float quotient, see .basics
82from __future__ import division as _; del _ # PYCHOK semicolon
84from pygeodesy.basics import len2, map1, map2, isodd, ub2str as _ub2str
85from pygeodesy.constants import EPS, _float as _F, _0_0, _1_0, _180_0, _360_0
86# from pygeodesy.datums import _ellipsoidal_datum # from .heights
87from pygeodesy.dms import parseDMS2
88from pygeodesy.errors import _incompatible, LenError, RangeError, _SciPyIssue
89from pygeodesy.fmath import favg, Fdot, fdot, Fhorner, frange
90from pygeodesy.heights import _as_llis2, _ascalar, _height_called, HeightError, \
91 _HeightsBase, _ellipsoidal_datum, _Wrap
92from pygeodesy.interns import NN, _4_, _COLONSPACE_, _COMMASPACE_, _cubic_, \
93 _DOT_, _E_, _height_, _in_, _kind_, _knots_, \
94 _lat_, _linear_, _lon_, _mean_, _N_, _n_a_, _not_, \
95 _numpy_, _on_, _outside_, _S_, _s_, _scipy_, \
96 _SPACE_, _stdev_, _supported_, _tbd_, _W_, _width_
97from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, _FOR_DOCS
98from pygeodesy.named import _Named, _NamedTuple, notOverloaded
99# from pygeodesy.namedTuples import LatLon3Tuple # _MODS
100from pygeodesy.props import Property_RO, property_RO
101from pygeodesy.streprs import attrs, Fmt, fstr, pairs
102from pygeodesy.units import Height, Int_, Lat, Lon
103# from pygeodesy.utils import _Wrap # from .heights
105from math import floor
106import os.path as _os_path
107from os import SEEK_CUR as _SEEK_CUR, SEEK_SET as _SEEK_SET
108from struct import calcsize as _calcsize, unpack as _unpack
110try:
111 from StringIO import StringIO as _BytesIO # reads bytes
112 _ub2str = str # PYCHOK convert bytes to str for egm*.pgm text
114except ImportError: # Python 3+
115 from io import BytesIO as _BytesIO # PYCHOK expected
117__all__ = _ALL_LAZY.geoids
118__version__ = '23.05.26'
120_assert_ = 'assert'
121_bHASH_ = b'#'
122_endian_ = 'endian'
123_format_ = '%s %r'
124_header_ = 'header'
125# temporarily hold a single instance for each int value
126_intCs = {}
127_interp2d_ks = {-2: _linear_,
128 -3: _cubic_,
129 -5: 'quintic'}
130_lli_ = 'lli'
131_non_increasing_ = 'non-increasing'
132_rb_ = 'rb'
135class _GeoidBase(_HeightsBase):
136 '''(INTERNAL) Base class for C{Geoid...}s.
137 '''
138 _cropped = None
139# _datum = _WGS84 # from _HeightsBase
140 _egm = None # open C{egm*.pgm} geoid file
141 _endian = _tbd_
142 _geoid = _n_a_
143 _hs_y_x = None # numpy 2darray, row-major order
144 _interp2d = None # interp2d interpolation
145 _kind = 3 # order for interp2d, RectBivariateSpline
146# _kmin = 2 # min number of knots
147 _knots = 0 # nlat * nlon
148 _mean = None # fixed in GeoidKarney
149# _name = NN # _Named
150 _nBytes = 0 # numpy size in bytes, float64
151 _pgm = None # PGM attributes, C{_PGM} or C{None}
152 _sizeB = 0 # geoid file size in bytes
153 _smooth = 0 # used only for RectBivariateSpline
154 _stdev = None # fixed in GeoidKarney
155 _u2B = 0 # np.itemsize or undefined
157 _lat_d = _0_0 # increment, +tive
158 _lat_lo = _0_0 # lower lat, south
159 _lat_hi = _0_0 # upper lat, noth
160 _lon_d = _0_0 # increment, +tive
161 _lon_lo = _0_0 # left lon, west
162 _lon_hi = _0_0 # right lon, east
163 _lon_of = _0_0 # forward lon offset
164 _lon_og = _0_0 # reverse lon offset
166 _center = None # (lat, lon, height)
167 _yx_hits = None # cache hits, ala Karney
169 def __init__(self, hs, p):
170 '''(INTERNAL) Set up the grid axes, the C{SciPy} interpolator
171 and several internal geoid attributes.
173 @arg hs: Grid knots with known height (C{numpy 2darray}).
174 @arg p: The C{slat, wlon, nlat, nlon, dlat, dlon} and
175 other geoid parameters (C{INTERNAL}).
177 @raise GeoidError: Incompatible grid B{C{hs}} shape or
178 invalid B{C{kind}}.
180 @raise LenError: Mismatch grid B{C{hs}} axis.
182 @raise SciPyError: A C{scipy.interpolate.inter2d} or
183 C{-.RectBivariateSpline} issue.
185 @raise SciPyWarning: A C{scipy.interpolate.inter2d} or
186 C{-.RectBivariateSpline} warning as
187 exception.
188 '''
189 spi = self.scipy_interpolate
190 # for 2d scipy.interpolate.interp2d(xs, ys, hs, ...) and
191 # scipy.interpolate.RectBivariateSpline(ys, xs, hs, ...)
192 # require the shape of hs to be (len(ys), len(xs)), note
193 # the different (xs, ys, ...) and (ys, xs, ...) orders
194 if (p.nlat, p.nlon) != hs.shape:
195 raise GeoidError(shape=hs.shape, txt=_incompatible((p.nlat, p.nlon)))
197 # both axes and bounding box
198 ys, self._lat_d = self._gaxis2(p.slat, p.dlat, p.nlat, _lat_ + _s_)
199 xs, self._lon_d = self._gaxis2(p.wlon, p.dlon, p.nlon, _lon_ + _s_)
201 bb = ys[0], ys[-1], xs[0], xs[-1] + p.dlon # fudge lon_hi
202 # geoid grids are typically stored in row-major order, some
203 # with rows (90..-90) reversed and columns (0..360) wrapped
204 # to Easten longitude, 0 <= east < 180 and 180 <= west < 360
205 k = self.kind
206 if k in _interp2d_ks:
207 self._interp2d = spi.interp2d(xs, ys, hs, kind=_interp2d_ks[k])
208 elif 1 <= k <= 5:
209 self._ev = spi.RectBivariateSpline(ys, xs, hs, bbox=bb, ky=k, kx=k,
210 s=self._smooth).ev
211 else:
212 raise GeoidError(kind=k)
214 self._hs_y_x = hs # numpy 2darray, row-major
215 self._nBytes = hs.nbytes # numpy size in bytes
216 self._knots = p.knots # grid knots
217 self._lon_of = float(p.flon) # forward offset
218 self._lon_og = float(p.glon) # reverse offset
219 # shrink the box by 1 unit on every side
220 # bb += self._lat_d, -self._lat_d, self._lon_d, -self._lon_d
221 self._lat_lo = float(bb[0])
222 self._lat_hi = float(bb[1])
223 self._lon_lo = float(bb[2] - p.glon)
224 self._lon_hi = float(bb[3] - p.glon)
226 def __call__(self, *llis, **wrap):
227 '''Interpolate the geoid height for one or several locations.
229 @arg llis: The location or locations (C{LatLon}, ... or
230 C{LatLon}s).
231 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
232 locations (C{bool}).
234 @return: A single interpolated geoid height (C{float}) or
235 a list or tuple of interpolated geoid heights
236 (C{float}s).
238 @raise GeoidError: Insufficient number of B{C{llis}}, an
239 invalid B{C{lli}} or the C{egm*.pgm}
240 geoid file is closed.
242 @raise RangeError: An B{C{lli}} is outside this geoid's lat-
243 or longitude range.
245 @raise SciPyError: A C{scipy.interpolate.inter2d} or
246 C{-.RectBivariateSpline} issue.
248 @raise SciPyWarning: A C{scipy.interpolate.inter2d} or
249 C{-.RectBivariateSpline} warning as
250 exception.
251 '''
252 return self._called(llis, True, **wrap)
254 def __enter__(self):
255 '''Open context.
256 '''
257 return self
259 def __exit__(self, *unused): # PYCHOK exc_type, exc_value, exc_traceback)
260 '''Close context.
261 '''
262 self.close()
263 # return None # XXX False
265 def __repr__(self):
266 return self.toStr()
268 def __str__(self):
269 return Fmt.PAREN(self.classname, repr(self.name))
271 def _called(self, llis, scipy, wrap=False):
272 # handle __call__
273 _as, llis = _as_llis2(llis, Error=GeoidError)
274 try:
275 hs, _w = [], _Wrap._latlonop(wrap)
276 for i, lli in enumerate(llis):
277 hs.append(self._hGeoid(*_w(lli.lat, lli.lon)))
278 return _as(hs)
280 except (GeoidError, RangeError) as x:
281 # XXX avoid str(LatLon()) degree symbols
282 t = _lli_ if _as is _ascalar else Fmt.SQUARE(llis=i)
283 lli = fstr((lli.lat, lli.lon), strepr=repr)
284 raise type(x)(t, lli, wrap=wrap, cause=x)
285 except Exception as x:
286 if scipy and self.scipy:
287 raise _SciPyIssue(x)
288 else:
289 raise
291 @Property_RO
292 def _center(self):
293 ''' Cache for method L{center}.
294 '''
295 return self._llh3(favg(self._lat_lo, self._lat_hi),
296 favg(self._lon_lo, self._lon_hi))
298 def center(self, LatLon=None):
299 '''Return the center location and height of this geoid.
301 @kwarg LatLon: Optional class to return the location and height
302 (C{LatLon}) or C{None}.
304 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
305 lon, height)} otherwise a B{C{LatLon}} instance
306 with the lat-, longitude and height of the grid
307 center location.
308 '''
309 return self._llh3LL(self._center, LatLon)
311 def close(self):
312 '''Close the C{egm*.pgm} geoid file if open (and applicable).
313 '''
314 if not self.closed:
315 self._egm.close()
316 self._egm = None
318 @property_RO
319 def closed(self):
320 '''Get the C{egm*.pgm} geoid file status.
321 '''
322 return self._egm is None
324 @Property_RO
325 def cropped(self):
326 '''Is geoid cropped (C{bool} or C{None} if crop not supported).
327 '''
328 return self._cropped
330 @Property_RO
331 def dtype(self):
332 '''Get the grid C{scipy} U{dtype<https://docs.SciPy.org/doc/numpy/
333 reference/generated/numpy.ndarray.dtype.html>} (C{numpy.dtype}).
334 '''
335 return self._hs_y_x.dtype
337 @Property_RO
338 def endian(self):
339 '''Get the geoid endianess and U{dtype<https://docs.SciPy.org/
340 doc/numpy/reference/generated/numpy.dtype.html>} (C{str}).
341 '''
342 return self._endian
344 def _ev(self, y, x): # PYCHOK expected
345 # only used for .interpolate.interp2d, but
346 # overwritten for .RectBivariateSpline,
347 # note (y, x) must be flipped!
348 return self._interp2d(x, y)
350 def _gaxis2(self, lo, d, n, name):
351 # build grid axis, hi = lo + (n - 1) * d
352 m, a = len2(frange(lo, n, d))
353 if m != n:
354 raise LenError(self.__class__, grid=m, **{name: n})
355 if d < 0:
356 d, a = -d, list(reversed(a))
357 for i in range(1, m):
358 e = a[i] - a[i-1]
359 if e < EPS: # non-increasing axis
360 i = Fmt.SQUARE(name, i)
361 raise GeoidError(i, e, txt=_non_increasing_)
362 return self.numpy.array(a), d
364 def _g2ll2(self, lat, lon): # PYCHOK no cover
365 notOverloaded(self, lat, lon)
367 def _gyx2g2(self, y, x):
368 # convert grid (y, x) indices to grid (lat, lon)
369 return ((self._lat_lo + self._lat_d * y),
370 (self._lon_lo + self._lon_of + self._lon_d * x))
372 def height(self, lats, lons, **wrap):
373 '''Interpolate the geoid height for one or several lat-/longitudes.
375 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
376 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
377 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
378 and B{C{lons}} locations (C{bool}).
380 @return: A single interpolated geoid height (C{float}) or a
381 list of interpolated geoid heights (C{float}s).
383 @raise GeoidError: Insufficient or non-matching number of
384 B{C{lats}} and B{C{lons}}.
386 @raise RangeError: A B{C{lat}} or B{C{lon}} is outside this
387 geoid's lat- or longitude range.
389 @raise SciPyError: A C{scipy.interpolate.inter2d} or
390 C{-.RectBivariateSpline} issue.
392 @raise SciPyWarning: A C{scipy.interpolate.inter2d} or
393 C{-.RectBivariateSpline} warning as
394 exception.
395 '''
396 return _height_called(self, lats, lons, Error=GeoidError, **wrap)
398 def _hGeoid(self, lat, lon):
399 out = self.outside(lat, lon)
400 if out:
401 lli = fstr((lat, lon), strepr=repr)
402 raise RangeError(lli=lli, txt=_SPACE_(_outside_, _on_, out))
403 return float(self._ev(*self._ll2g2(lat, lon)))
405 @Property_RO
406 def _highest(self):
407 '''(INTERNAL) Cache for L{highest} method.
408 '''
409 return self._llh3minmax(True)
411 def highest(self, LatLon=None, **unused):
412 '''Return the location and largest height of this geoid.
414 @kwarg LatLon: Optional class to return the location and height
415 (C{LatLon}) or C{None}.
417 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
418 lon, height)} otherwise a B{C{LatLon}} instance
419 with the lat-, longitude and height of the highest
420 grid location.
421 '''
422 return self._llh3LL(self._highest, LatLon)
424 @Property_RO
425 def hits(self):
426 '''Get the number of cache hits (C{int} or C{None}).
427 '''
428 return self._yx_hits
430 @Property_RO
431 def kind(self):
432 '''Get the interpolator kind and order (C{int}).
433 '''
434 return self._kind
436 @Property_RO
437 def knots(self):
438 '''Get the number of grid knots (C{int}).
439 '''
440 return self._knots
442 def _ll2g2(self, lat, lon): # PYCHOK no cover
443 notOverloaded(self, lat, lon)
445 @property_RO
446 def _LL3T(self):
447 '''(INTERNAL) Get L{LatLon3Tuple} once.
448 '''
449 T = _MODS.namedTuples.LatLon3Tuple
450 _GeoidBase._LL3T = T # overwrite poperty_RO
451 return T
453 def _llh3(self, lat, lon):
454 return self._LL3T(lat, lon, self._hGeoid(lat, lon), name=self.name)
456 def _llh3LL(self, llh, LatLon):
457 return llh if LatLon is None else self._xnamed(LatLon(*llh))
459 def _llh3minmax(self, highest=True, *unused):
460 hs, np = self._hs_y_x, self.numpy
461 # <https://docs.SciPy.org/doc/numpy/reference/generated/
462 # numpy.argmin.html#numpy.argmin>
463 arg = np.argmax if highest else np.argmin
464 y, x = np.unravel_index(arg(hs, axis=None), hs.shape)
465 return self._g2ll2(*self._gyx2g2(y, x)) + (float(hs[y, x]),)
467 def _load(self, g, dtype, n, offset=0):
468 # numpy.fromfile, like .frombuffer
469 g.seek(offset, _SEEK_SET)
470 return self.numpy.fromfile(g, dtype, n)
472 @Property_RO
473 def _lowerleft(self):
474 '''(INTERNAL) Cache for L{lowerleft}.
475 '''
476 return self._llh3(self._lat_lo, self._lon_lo)
478 def lowerleft(self, LatLon=None):
479 '''Return the lower-left location and height of this geoid.
481 @kwarg LatLon: Optional class to return the location
482 (C{LatLon}) and height or C{None}.
484 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
485 lon, height)} otherwise a B{C{LatLon}} instance
486 with the lat-, longitude and height of the lower-left,
487 SW grid corner.
488 '''
489 return self._llh3LL(self._lowerleft, LatLon)
491 @Property_RO
492 def _loweright(self):
493 '''(INTERNAL) Cache for L{loweright}.
494 '''
495 return self._llh3(self._lat_lo, self._lon_hi)
497 def loweright(self, LatLon=None):
498 '''Return the lower-right location and height of this geoid.
500 @kwarg LatLon: Optional class to return the location and height
501 (C{LatLon}) or C{None}.
503 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
504 lon, height)} otherwise a B{C{LatLon}} instance
505 with the lat-, longitude and height of the lower-right,
506 SE grid corner.
507 '''
509 return self._llh3LL(self._loweright, LatLon)
511 lowerright = loweright # synonymous
513 @Property_RO
514 def _lowest(self):
515 '''(INTERNAL) Cache for L{lowest}.
516 '''
517 return self._llh3minmax(False)
519 def lowest(self, LatLon=None, **unused):
520 '''Return the location and lowest height of this geoid.
522 @kwarg LatLon: Optional class to return the location and height
523 (C{LatLon}) or C{None}.
525 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
526 lon, height)} otherwise a B{C{LatLon}} instance
527 with the lat-, longitude and height of the lowest
528 grid location.
529 '''
530 return self._llh3LL(self._lowest, LatLon)
532 @Property_RO
533 def mean(self):
534 '''Get the mean of this geoid's heights (C{float}).
535 '''
536 if self._mean is None: # see GeoidKarney
537 self._mean = float(self.numpy.mean(self._hs_y_x))
538 return self._mean
540 @property_RO
541 def name(self):
542 '''Get the name of this geoid (C{str}).
543 '''
544 return _HeightsBase.name.fget(self) or self._geoid # recursion
546 @Property_RO
547 def nBytes(self):
548 '''Get the grid in-memory size in bytes (C{int}).
549 '''
550 return self._nBytes
552 def _open(self, geoid, datum, kind, name, smooth):
553 # open the geoid file
554 try:
555 self._geoid = _os_path.basename(geoid)
556 self._sizeB = _os_path.getsize(geoid)
557 g = open(geoid, _rb_)
558 except (IOError, OSError) as x:
559 raise GeoidError(geoid=geoid, cause=x)
561 if datum not in (None, self._datum):
562 self._datum = _ellipsoidal_datum(datum, name=name)
563 self._kind = int(kind)
564 if name:
565 _HeightsBase.name.fset(self, name) # rename
566 if smooth:
567 self._smooth = Int_(smooth=smooth, Error=GeoidError, low=0)
569 return g
571 def outside(self, lat, lon):
572 '''Check whether a location is outside this geoid's
573 lat-/longitude or crop range.
575 @arg lat: The latitude (C{degrees}).
576 @arg lon: The longitude (C{degrees}).
578 @return: A 1- or 2-character C{str} if outside or an
579 empty C{str} if inside.
580 '''
581 return (_S_ if lat < self._lat_lo else
582 (_N_ if lat > self._lat_hi else NN)) + \
583 (_W_ if lon < self._lon_lo else
584 (_E_ if lon > self._lon_hi else NN))
586 @Property_RO
587 def pgm(self):
588 '''Get the PGM attributes (C{_PGM} or C{None} if not available/applicable).
589 '''
590 return self._pgm
592 @Property_RO
593 def sizeB(self):
594 '''Get the geoid grid file size in bytes (C{int}).
595 '''
596 return self._sizeB
598 @Property_RO
599 def smooth(self):
600 '''Get the C{RectBivariateSpline} smoothing (C{int}).
601 '''
602 return self._smooth
604 @Property_RO
605 def stdev(self):
606 '''Get the standard deviation of this geoid's heights (C{float}) or C{None}.
607 '''
608 if self._stdev is None: # see GeoidKarney
609 self._stdev = float(self.numpy.std(self._hs_y_x))
610 return self._stdev
612 def _swne(self, crop):
613 # crop box to 4-tuple (s, w, n, e)
614 try:
615 if len(crop) == 2:
616 try: # sw, ne LatLons
617 swne = (crop[0].lat, crop[0].lon,
618 crop[1].lat, crop[1].lon)
619 except AttributeError: # (s, w), (n, e)
620 swne = tuple(crop[0]) + tuple(crop[1])
621 else: # (s, w, n, e)
622 swne = crop
623 if len(swne) == 4:
624 s, w, n, e = map(float, swne)
625 if -90 <= s <= (n - _1_0) <= 89 and \
626 -180 <= w <= (e - _1_0) <= 179:
627 return s, w, n, e
628 except (IndexError, TypeError, ValueError):
629 pass
630 raise GeoidError(crop=crop)
632 def toStr(self, prec=3, sep=_COMMASPACE_): # PYCHOK signature
633 '''This geoid and all geoid attributes as a string.
635 @kwarg prec: Number of decimal digits (0..9 or C{None} for
636 default). Trailing zero decimals are stripped
637 for B{C{prec}} values of 1 and above, but kept
638 for negative B{C{prec}} values.
639 @kwarg sep: Separator to join (C{str}).
641 @return: Geoid name and attributes (C{str}).
642 '''
643 s = 1 if self.kind < 0 else 2
644 t = tuple(Fmt.PAREN(m.__name__, fstr(m(), prec=prec)) for m in
645 (self.lowerleft, self.upperright,
646 self.center,
647 self.highest, self.lowest)) + \
648 attrs( _mean_, _stdev_, prec=prec, Nones=False) + \
649 attrs((_kind_, 'smooth')[:s], prec=prec, Nones=False) + \
650 attrs( 'cropped', 'dtype', _endian_, 'hits', _knots_, 'nBytes',
651 'sizeB', _scipy_, _numpy_, prec=prec, Nones=False)
652 return _COLONSPACE_(self, sep.join(t))
654 @Property_RO
655 def u2B(self):
656 '''Get the PGM itemsize in bytes (C{int}).
657 '''
658 return self._u2B
660 @Property_RO
661 def _upperleft(self):
662 '''(INTERNAL) Cache for method L{upperleft}.
663 '''
664 return self._llh3(self._lat_hi, self._lon_lo)
666 def upperleft(self, LatLon=None):
667 '''Return the upper-left location and height of this geoid.
669 @kwarg LatLon: Optional class to return the location and height
670 (C{LatLon}) or C{None}.
672 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
673 lon, height)} otherwise a B{C{LatLon}} instance
674 with the lat-, longitude and height of the upper-left,
675 NW grid corner.
676 '''
677 return self._llh3LL(self._upperleft, LatLon)
679 @Property_RO
680 def _upperright(self):
681 '''(INTERNAL) Cache for method L{upperright}.
682 '''
683 return self._llh3(self._lat_hi, self._lon_hi)
685 def upperright(self, LatLon=None):
686 '''Return the upper-right location and height of this geoid.
688 @kwarg LatLon: Optional class to return the location and height
689 (C{LatLon}) or C{None}.
691 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
692 lon, height)} otherwise a B{C{LatLon}} instance
693 with the lat-, longitude and height of the upper-right,
694 NE grid corner.
695 '''
696 return self._llh3LL(self._upperright, LatLon)
699class GeoidError(HeightError):
700 '''Geoid interpolator C{Geoid...} or interpolation issue.
701 '''
702 pass
705class GeoidG2012B(_GeoidBase):
706 '''Geoid height interpolator for U{GEOID12B Model
707 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/>} grids U{CONUS
708 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_CONUS.shtml>},
709 U{Alaska<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_AK.shtml>},
710 U{Hawaii<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_HI.shtml>},
711 U{Guam and Northern Mariana Islands
712 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_GMNI.shtml>},
713 U{Puerto Rico and U.S. Virgin Islands
714 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_PRVI.shtml>} and
715 U{American Samoa<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_AS.shtml>}
716 based on C{SciPy} U{RectBivariateSpline<https://docs.SciPy.org/doc/
717 scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>}
718 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/
719 scipy.interpolate.interp2d.html>} interpolation.
721 Use any of the binary C{le} (little endian) or C{be} (big endian)
722 C{g2012b*.bin} grid files.
723 '''
724 def __init__(self, g2012b_bin, crop=None, datum=None, # NAD 83 Ellipsoid
725 kind=3, name=NN, smooth=0):
726 '''New L{GeoidG2012B} interpolator.
728 @arg g2012b_bin: A C{GEOID12B} grid file name (C{.bin}).
729 @kwarg crop: Optional crop box, not supported (C{None}).
730 @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2}
731 or L{a_f2Tuple}), default C{WGS84}.
732 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for
733 U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/
734 reference/generated/scipy.interpolate.RectBivariateSpline.html>},
735 -2 for U{interp2d linear<https://docs.SciPy.org/doc/scipy/
736 reference/generated/scipy.interpolate.interp2d.html>}, -3
737 for C{interp2d cubic} or -5 for C{interp2d quintic}.
738 @kwarg name: Optional geoid name (C{str}).
739 @kwarg smooth: Smoothing factor for U{RectBivariateSpline
740 <https://docs.SciPy.org/doc/scipy/reference/generated/
741 scipy.interpolate.RectBivariateSpline.html>}
742 only (C{int}).
744 @raise GeoidError: G2012B grid file B{C{g2012b_bin}} issue or invalid
745 B{C{crop}}, B{C{kind}} or B{C{smooth}}.
747 @raise ImportError: Package C{numpy} or C{scipy} not found or
748 not installed.
750 @raise LenError: Grid file B{C{g2012b_bin}} axis mismatch.
752 @raise SciPyError: A C{RectBivariateSpline} or C{inter2d} issue.
754 @raise SciPyWarning: A C{RectBivariateSpline} or C{inter2d}
755 warning as exception.
757 @raise TypeError: Invalid B{C{datum}}.
758 '''
759 if crop is not None:
760 raise GeoidError(crop=crop, txt=_not_(_supported_))
762 g = self._open(g2012b_bin, datum, kind, name, smooth)
763 _ = self.numpy # import numpy for ._load and
765 try:
766 p = _Gpars()
767 n = (self.sizeB // 4) - 11 # number of f4 heights
768 # U{numpy dtype formats are different from Python struct formats
769 # <https://docs.SciPy.org/doc/numpy-1.15.0/reference/arrays.dtypes.html>}
770 for en_ in ('<', '>'):
771 # skip 4xf8, get 3xi4
772 p.nlat, p.nlon, ien = map(int, self._load(g, en_+'i4', 3, 32))
773 if ien == 1: # correct endian
774 p.knots = p.nlat * p.nlon
775 if p.knots == n and 1 < p.nlat < n \
776 and 1 < p.nlon < n:
777 self._endian = en_+'f4'
778 break
779 else: # couldn't validate endian
780 raise GeoidError(_endian_)
782 # get the first 4xf8
783 p.slat, p.wlon, p.dlat, p.dlon = map(float, self._load(g, en_+'f8', 4))
784 # read all f4 heights, ignoring the first 4xf8 and 3xi4
785 hs = self._load(g, self._endian, n, 44).reshape(p.nlat, p.nlon)
786 p.wlon -= _360_0 # western-most East longitude to earth (..., lon)
787 _GeoidBase.__init__(self, hs, p)
789 except Exception as x:
790 raise _SciPyIssue(x, _in_, repr(g2012b_bin))
791 finally:
792 g.close()
794 def _g2ll2(self, lat, lon):
795 # convert grid (lat, lon) to earth (lat, lon)
796 return lat, lon
798 def _ll2g2(self, lat, lon):
799 # convert earth (lat, lon) to grid (lat, lon)
800 return lat, lon
802 if _FOR_DOCS:
803 __call__ = _GeoidBase.__call__
804 height = _GeoidBase.height
807class GeoidHeight5Tuple(_NamedTuple): # .geoids.py
808 '''5-Tuple C{(lat, lon, egm84, egm96, egm2008)} for U{GeoidHeights.dat
809 <https://SourceForge.net/projects/geographiclib/files/testdata/>}
810 tests with the heights for 3 different EGM grids at C{degrees90}
811 and C{degrees180} degrees (after converting C{lon} from original
812 C{0 <= EasterLon <= 360}).
813 '''
814 _Names_ = (_lat_, _lon_, 'egm84', 'egm96', 'egm2008')
815 _Units_ = ( Lat, Lon, Height, Height, Height)
818def _I(i):
819 '''(INTERNAL) Cache a single C{int} constant.
820 '''
821 return _intCs.setdefault(i, i) # PYCHOK undefined due to del _intCs
824def _T(*cs):
825 '''(INTERNAL) Cache a tuple of single C{int} constants.
826 '''
827 return map1(_I, *cs)
830class GeoidKarney(_GeoidBase):
831 '''Geoid height interpolator for I{Karney}'s U{GeographicLib Earth
832 Gravitational Model (EGM)<https://GeographicLib.SourceForge.io/C++/doc/
833 geoid.html>} geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/
834 html/geoid.html#geoidinst>} datasets using bilinear or U{cubic
835 <https://dl.ACM.org/citation.cfm?id=368443>} interpolation and U{caching
836 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidcache>}
837 in pure Python, transcoded from I{Karney}'s U{C++ class Geoid
838 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinterp>}.
840 Use any of the geoid U{egm84-, egm96- or egm2008-*.pgm
841 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>}
842 datasets.
843 '''
844 _C0 = _F(372), _F(240), _F(372) # n, _ and s common denominators
845 # matrices c3n_, c3, c3s_, transposed from GeographicLib/Geoid.cpp
846 _C3 = ((_T(0, 0, 62, 124, 124, 62, 0, 0, 0, 0, 0, 0),
847 _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
848 _T(-131, 7, -31, -62, -62, -31, 45, 216, 156, -45, -55, -7),
849 _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
850 _T(138, -138, 0, 0, 0, 0, -183, 33, 153, -3, 48, -48), # PYCHOK indent
851 _T(144, 42, -62, -124, -124, -62, -9, 87, 99, 9, 42, -42),
852 _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
853 _T(0, 0, 0, 0, 0, 0, 93, -93, -93, 93, 0, 0),
854 _T(-102, 102, 0, 0, 0, 0, 18, 12, -12, -18, -84, 84),
855 _T(-31, -31, 31, 62, 62, 31, 0, -93, -93, 0, 31, 31)), # PYCHOK indent
857 (_T(9, -9, 9, 186, 54, -9, -9, 54, -54, 9, -9, 9),
858 _T(-18, 18, -88, -42, 162, -32, 8, -78, 78, -8, 18, -18),
859 _T(-88, 8, -18, -42, -78, 18, 18, 162, 78, -18, -32, -8),
860 _T(0, 0, 90, -150, 30, 30, 30, -90, 90, -30, 0, 0),
861 _T(96, -96, 96, -96, -24, 24, -96, -24, 144, -24, 24, -24), # PYCHOK indent
862 _T(90, 30, 0, -150, -90, 0, 0, 30, 90, 0, 30, -30),
863 _T(0, 0, -20, 60, -60, 20, -20, 60, -60, 20, 0, 0),
864 _T(0, 0, -60, 60, 60, -60, 60, -60, -60, 60, 0, 0),
865 _T(-60, 60, 0, 60, -60, 0, 0, 60, -60, 0, -60, 60),
866 _T(-20, -20, 0, 60, 60, 0, 0, -60, -60, 0, 20, 20)),
868 (_T(18, -18, 36, 210, 162, -36, 0, 0, 0, 0, -18, 18), # PYCHOK indent
869 _T(-36, 36, -165, 45, 141, -21, 0, 0, 0, 0, 36, -36),
870 _T(-122, -2, -27, -111, -75, 27, 62, 124, 124, 62, -64, 2),
871 _T(0, 0, 93, -93, -93, 93, 0, 0, 0, 0, 0, 0),
872 _T(120, -120, 147, -57, -129, 39, 0, 0, 0, 0, 66, -66), # PYCHOK indent
873 _T(135, 51, -9, -192, -180, 9, 31, 62, 62, 31, 51, -51),
874 _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
875 _T(0, 0, -93, 93, 93, -93, 0, 0, 0, 0, 0, 0),
876 _T(-84, 84, 18, 12, -12, -18, 0, 0, 0, 0, -102, 102),
877 _T(-31, -31, 0, 93, 93, 0, -31, -62, -62, -31, 31, 31)))
879 _BT = (_T(0, 0), # bilinear 4-tuple [i, j] indices
880 _T(1, 0),
881 _T(0, 1),
882 _T(1, 1))
884 _CM = (_T( 0, -1), # 10x12 cubic matrix [i, j] indices
885 _T( 1, -1),
886 _T(-1, 0),
887 _T( 0, 0),
888 _T( 1, 0),
889 _T( 2, 0),
890 _T(-1, 1),
891 _T( 0, 1),
892 _T( 1, 1),
893 _T( 2, 1),
894 _T( 0, 2),
895 _T( 1, 2))
897 _endian = '>H' # struct.unpack 1 ushort (big endian, unsigned short)
898 _4endian = '>4H' # struct.unpack 4 ushorts
899 _Rendian = NN # struct.unpack a row of ushorts
900# _highest = (-8.4, 147.367, 85.839) if egm2008-1.pgm else (
901# (-8.167, 147.25, 85.422) if egm96-5.pgm else
902# (-4.5, 148.75, 81.33)) # egm84-15.pgm
903# _lowest = (4.7, 78.767, -106.911) if egm2008-1.pgm else (
904# (4.667, 78.833, -107.043) if egm96-5.pgm else
905# (4.75, 79.25, -107.34)) # egm84-15.pgm
906 _mean = _F(-1.317) # from egm2008-1, -1.438 egm96-5, -0.855 egm84-15
907 _nBytes = None # not applicable
908 _nterms = len(_C3[0]) # columns length, number of row
909 _smooth = None # not applicable
910 _stdev = _F(29.244) # from egm2008-1, 29.227 egm96-5, 29.183 egm84-15
911 _u2B = _calcsize(_endian) # pixelsize_ in bytes
912 _4u2B = _calcsize(_4endian) # 4 pixelsize_s in bytes
913 _Ru2B = 0 # row of pixelsize_s in bytes
914 _yxH = () # cache (y, x) indices
915 _yxHt = () # cached 4- or 10-tuple for _ev2H resp. _ev3H
916 _yx_hits = 0 # cache hits
918 def __init__(self, egm_pgm, crop=None, datum=None, # WGS84
919 kind=3, name=NN, smooth=None):
920 '''New L{GeoidKarney} interpolator.
922 @arg egm_pgm: An U{EGM geoid dataset<https://GeographicLib.SourceForge.io/
923 html/geoid.html#geoidinst>} file name (C{egm*.pgm}), see
924 note below.
925 @kwarg crop: Optional box to limit geoid locations, a 4-tuple (C{south,
926 west, north, east}), 2-tuple (C{(south, west), (north,
927 east)}) or 2, in C{degrees90} lat- and C{degrees180}
928 longitudes or a 2-tuple (C{LatLonSW, LatLonNE}) of
929 C{LatLon} instances.
930 @kwarg datum: Optional grid datum (C{Datum}, L{Ellipsoid}, L{Ellipsoid2}
931 or L{a_f2Tuple}), default C{WGS84}.
932 @kwarg kind: Interpolation order (C{int}), 2 for C{bilinear} or 3
933 for C{cubic}.
934 @kwarg name: Optional geoid name (C{str}).
935 @kwarg smooth: Smoothing factor, unsupported (C{None}).
937 @raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid
938 B{C{crop}}, B{C{kind}} or B{C{smooth}}.
940 @raise TypeError: Invalid B{C{datum}}.
942 @see: Class L{GeoidPGM} and function L{egmGeoidHeights}.
944 @note: Geoid file B{C{egm_pgm}} remains open and must be closed
945 by calling the C{close} method or by using this instance
946 in a C{with B{GeoidKarney}(...) as ...} context.
947 '''
948 if smooth is not None:
949 raise GeoidError(smooth=smooth, txt=_not_(_supported_))
951 if kind in (2,):
952 self._evH = self._ev2H
953 elif kind not in (3,):
954 raise GeoidError(kind=kind)
956 self._egm = g = self._open(egm_pgm, datum, kind, name, smooth)
957 self._pgm = p = _PGM(g, pgm=egm_pgm, itemsize=self.u2B, sizeB=self.sizeB)
959 self._Rendian = self._4endian.replace(_4_, str(p.nlon))
960 self._Ru2B = _calcsize(self._Rendian)
962 self._knots = p.knots # grid knots
963 self._lon_of = float(p.flon) # forward offset
964 self._lon_og = float(p.glon) # reverse offset
965 # set earth (lat, lon) limits (s, w, n, e)
966 self._lat_lo, \
967 self._lon_lo, \
968 self._lat_hi, \
969 self._lon_hi = self._swne(crop if crop else p.crop4)
970 self._cropped = True if crop else False
972 def __call__(self, *llis, **wrap):
973 '''Interpolate the geoid height for one or several locations.
975 @arg llis: The location or locations (C{LatLon}, ... or
976 C{LatLon}s).
977 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
978 locations (C{bool}).
980 @return: A single interpolated geoid height (C{float}) or
981 a list or tuple of interpolated geoid heights
982 (C{float}s).
984 @raise GeoidError: Insufficient number of B{C{llis}}, an
985 invalid B{C{lli}} or the C{egm*.pgm}
986 geoid file is closed.
988 @raise RangeError: An B{C{lli}} is outside this geoid's lat-
989 or longitude range.
990 '''
991 return self._called(llis, False, **wrap)
993 def _c0c3v(self, y, x):
994 # get the common denominator, the 10x12 cubic matrix and
995 # the 12 cubic v-coefficients around geoid index (y, x)
996 p = self._pgm
997 if 0 < x < (p.nlon - 2) and 0 < y < (p.nlat - 2):
998 # read 4x4 ushorts, drop the 4 corners
999 g = self._egm
1000 e = self._4endian
1001 n = self._4u2B
1002 R = self._Ru2B
1004 b = self._seek(y - 1, x - 1)
1005 v = _unpack(e, g.read(n))[1:3]
1006 b += R
1007 g.seek(b, _SEEK_SET)
1008 v += _unpack(e, g.read(n))
1009 b += R
1010 g.seek(b, _SEEK_SET)
1011 v += _unpack(e, g.read(n))
1012 b += R
1013 g.seek(b, _SEEK_SET)
1014 v += _unpack(e, g.read(n))[1:3]
1015 j = 1
1017 else: # likely some wrapped y and/or x's
1018 v = self._raws(y, x, GeoidKarney._CM)
1019 j = 0 if y < 1 else (1 if y < (p.nlat - 2) else 2)
1021 return GeoidKarney._C0[j], GeoidKarney._C3[j], v
1023 @Property_RO
1024 def dtype(self):
1025 '''Get the geoid's grid data type (C{str}).
1026 '''
1027 return 'ushort'
1029 def _ev(self, lat, lon): # PYCHOK expected
1030 # interpolate the geoid height at grid (lat, lon)
1031 fy, fx = self._g2yx2(lat, lon)
1032 y, x = int(floor(fy)), int(floor(fx))
1033 fy -= y
1034 fx -= x
1035 H = self._evH(fy, fx, y, x) # ._ev3H or ._ev2H
1036 H *= self._pgm.Scale # H.fmul(self._pgm.Scale)
1037 H += self._pgm.Offset # H.fadd(self._pgm.Offset)
1038 return H.fsum()
1040 def _ev2H(self, fy, fx, *yx):
1041 # compute the bilinear 4-tuple and interpolate raw H
1042 if self._yxH == yx:
1043 t = self._yxHt
1044 self._yx_hits += 1
1045 else:
1046 y, x = self._yxH = yx
1047 self._yxHt = t = self._raws(y, x, GeoidKarney._BT)
1048 v = _1_0, -fx, fx
1049 H = Fdot(v, t[0], t[0], t[1]).fmul(_1_0 - fy) # c = a * (1 - fy)
1050 H += Fdot(v, t[2], t[2], t[3]).fmul(fy) # c += b * fy
1051 return H
1053 def _ev3H(self, fy, fx, *yx):
1054 # compute the cubic 10-tuple and interpolate raw H
1055 if self._yxH == yx:
1056 t = self._yxHt
1057 self._yx_hits += 1
1058 else:
1059 self._yxH = yx
1060 c0, c3, v = self._c0c3v(*yx)
1061 t = [fdot(v, *c3[i]) / c0 for i in range(self._nterms)]
1062 self._yxHt = t = tuple(t)
1063 # GeographicLib/Geoid.cpp Geoid::height(lat, lon) ...
1064 # real h = t[0] + fx * (t[1] + fx * (t[3] + fx * t[6])) +
1065 # fy * (t[2] + fx * (t[4] + fx * t[7]) +
1066 # fy * (t[5] + fx * t[8] + fy * t[9]));
1067 v = _1_0, fx, fy
1068 H = Fdot(v, t[5], t[8], t[9])
1069 H *= fy
1070 H += Fhorner(fx, t[2], t[4], t[7])
1071 H *= fy
1072 H += Fhorner(fx, t[0], t[1], t[3], t[6])
1073 return H
1075 _evH = _ev3H # overriden for kind == 2
1077 def _g2ll2(self, lat, lon):
1078 # convert grid (lat, lon) to earth (lat, lon), uncropped
1079 while lon > _180_0:
1080 lon -= _360_0
1081 return lat, lon
1083 def _g2yx2(self, lat, lon):
1084 # convert grid (lat, lon) to grid (y, x) indices
1085 p = self._pgm
1086 # note, slat = +90, rlat < 0 makes y >=0
1087 return ((lat - p.slat) * p.rlat), ((lon - p.wlon) * p.rlon)
1089 def _gyx2g2(self, y, x):
1090 # convert grid (y, x) indices to grid (lat, lon)
1091 p = self._pgm
1092 return (p.slat + p.dlat * y), (p.wlon + p.dlon * x)
1094 def height(self, lats, lons, **wrap):
1095 '''Interpolate the geoid height for one or several lat-/longitudes.
1097 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
1098 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
1099 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
1100 and B{C{lons}} locations (C{bool}).
1102 @return: A single interpolated geoid height (C{float}) or a
1103 list of interpolated geoid heights (C{float}s).
1105 @raise GeoidError: Insufficient or non-matching number of
1106 B{C{lats}} and B{C{lons}} or the C{egm*.pgm}
1107 geoid file is closed.
1109 @raise RangeError: A B{C{lat}} or B{C{lon}} is outside this
1110 geoid's lat- or longitude range.
1111 '''
1112 return _height_called(self, lats, lons, Error=GeoidError, **wrap)
1114 @Property_RO
1115 def _highest_ltd(self):
1116 '''(INTERNAL) Cache for L{highest} mesthod.
1117 '''
1118 return self._llh3minmax(True, -12, -4)
1120 def highest(self, LatLon=None, full=False): # PYCHOK full
1121 '''Return the location and largest height of this geoid.
1123 @kwarg LatLon: Optional class to return the location and height
1124 (C{LatLon}) or C{None}.
1125 @kwarg full: Search the full or limited latitude range (C{bool}).
1127 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
1128 lon, height)} otherwise a B{C{LatLon}} instance
1129 with the lat-, longitude and height of the highest
1130 grid location.
1131 '''
1132 llh = self._highest if full or self.cropped else self._highest_ltd
1133 return self._llh3LL(llh, LatLon)
1135 def _lat2y2(self, lat2):
1136 # convert earth lat(s) to min and max grid y indices
1137 ys, m = [], self._pgm.nlat - 1
1138 for lat in lat2:
1139 y, _ = self._g2yx2(*self._ll2g2(lat, 0))
1140 ys.append(max(min(int(y), m), 0))
1141 return min(ys), max(ys) + 1
1143 def _ll2g2(self, lat, lon):
1144 # convert earth (lat, lon) to grid (lat, lon), uncropped
1145 while lon < 0:
1146 lon += _360_0
1147 return lat, lon
1149 def _llh3minmax(self, highest=True, *lat2):
1150 # find highest or lowest, takes 10+ secs for egm2008-1.pgm geoid
1151 # (Python 2.7.16, macOS 10.13.6 High Sierra, iMac 3 GHz Core i3)
1152 y = x = 0
1153 h = self._raw(y, x)
1154 if highest:
1155 for j, r in self._raw2(*lat2):
1156 m = max(r)
1157 if m > h:
1158 h, y, x = m, j, r.index(m)
1159 else: # lowest
1160 for j, r in self._raw2(*lat2):
1161 m = min(r)
1162 if m < h:
1163 h, y, x = m, j, r.index(m)
1164 h *= self._pgm.Scale
1165 h += self._pgm.Offset
1166 return self._g2ll2(*self._gyx2g2(y, x)) + (h,)
1168 @Property_RO
1169 def _lowest_ltd(self):
1170 '''(INTERNAL) Cache for L{lowest}.
1171 '''
1172 return self._llh3minmax(False, 0, 8)
1174 def lowest(self, LatLon=None, full=False): # PYCHOK full
1175 '''Return the location and lowest height of this geoid.
1177 @kwarg LatLon: Optional class to return the location and height
1178 (C{LatLon}) or C{None}.
1179 @kwarg full: Search the full or limited latitude range (C{bool}).
1181 @return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat,
1182 lon, height)} otherwise a B{C{LatLon}} instance
1183 with the lat-, longitude and height of the lowest
1184 grid location.
1185 '''
1186 llh = self._lowest if full or self.cropped else self._lowest_ltd
1187 return self._llh3LL(llh, LatLon)
1189 def _raw(self, y, x):
1190 # get the ushort geoid height at geoid index (y, x),
1191 # like GeographicLib/Geoid.hpp real rawval(is, iy)
1192 p = self._pgm
1193 if x < 0:
1194 x += p.nlon
1195 elif x >= p.nlon:
1196 x -= p.nlon
1197 h = p.nlon // 2
1198 if y < 0:
1199 y = -y
1200 elif y >= p.nlat:
1201 y = (p.nlat - 1) * 2 - y
1202 else:
1203 h = 0
1204 x += h if x < h else -h
1205 self._seek(y, x)
1206 h = _unpack(self._endian, self._egm.read(self._u2B))
1207 return h[0]
1209 def _raws(self, y, x, ijs):
1210 # get bilinear 4-tuple or 10x12 cubic matrix
1211 return tuple(self._raw(y + j, x + i) for i, j in ijs)
1213 def _raw2(self, *lat2):
1214 # yield a 2-tuple (y, ushorts) for each row or for
1215 # the rows between two (or more) earth lat values
1216 p = self._pgm
1217 g = self._egm
1218 e = self._Rendian
1219 n = self._Ru2B
1220 # min(lat2) <= lat <= max(lat2) or 0 <= y < p.nlat
1221 s, t = self._lat2y2(lat2) if lat2 else (0, p.nlat)
1222 self._seek(s, 0) # to start of row s
1223 for y in range(s, t):
1224 yield y, _unpack(e, g.read(n))
1226 def _seek(self, y, x):
1227 # position geoid to grid index (y, x)
1228 p, g = self._pgm, self._egm
1229 if g:
1230 b = p.skip + (y * p.nlon + x) * self._u2B
1231 g.seek(b, _SEEK_SET)
1232 return b # position
1233 raise GeoidError('closed file: %r' % (p.egm,)) # IOError
1236class GeoidPGM(_GeoidBase):
1237 '''Geoid height interpolator for I{Karney}'s U{GeographicLib Earth
1238 Gravitational Model (EGM)<https://GeographicLib.SourceForge.io/C++/doc/
1239 geoid.html>} geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/
1240 html/geoid.html#geoidinst>} datasets but based on C{SciPy}
1241 U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/reference/
1242 generated/scipy.interpolate.RectBivariateSpline.html>} or
1243 U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/
1244 scipy.interpolate.interp2d.html>} interpolation.
1246 Use any of the U{egm84-, egm96- or egm2008-*.pgm
1247 <https://GeographicLib.SourceForge.io/C++/doc/geoid.html#geoidinst>}
1248 datasets. However, unless cropped, an entire C{egm*.pgm} dataset
1249 is loaded into the C{SciPy} U{RectBivariateSpline<https://docs.SciPy.org/
1250 doc/scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>}
1251 or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/
1252 scipy.interpolate.interp2d.html>} interpolator and converted from
1253 2-byte C{int} to 8-byte C{dtype float64}. Therefore, internal memory
1254 usage is 4x the U{egm*.pgm<https://GeographicLib.SourceForge.io/C++/doc/
1255 geoid.html#geoidinst>} file size and may exceed the available memory,
1256 especially with 32-bit Python, see properties C{.nBytes} and C{.sizeB}.
1257 '''
1258 _endian = '>u2'
1260 def __init__(self, egm_pgm, crop=None, datum=None, # WGS84
1261 kind=3, name=NN, smooth=0):
1262 '''New L{GeoidPGM} interpolator.
1264 @arg egm_pgm: An U{EGM geoid dataset<https://GeographicLib.SourceForge.io/
1265 html/geoid.html#geoidinst>} file name (C{egm*.pgm}).
1266 @kwarg crop: Optional box to crop B{C{egm_pgm}}, a 4-tuple (C{south, west,
1267 north, east}) or 2-tuple (C{(south, west), (north, east)}),
1268 in C{degrees90} lat- and C{degrees180} longitudes or a
1269 2-tuple (C{LatLonSW, LatLonNE}) of C{LatLon} instances.
1270 @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2}
1271 or L{a_f2Tuple}), default C{WGS84}.
1272 @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for
1273 U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/
1274 reference/generated/scipy.interpolate.RectBivariateSpline.html>},
1275 -2 for U{interp2d linear<https://docs.SciPy.org/doc/scipy/
1276 reference/generated/scipy.interpolate.interp2d.html>}, -3
1277 for C{interp2d cubic} or -5 for C{interp2d quintic}.
1278 @kwarg name: Optional geoid name (C{str}).
1279 @kwarg smooth: Smoothing factor for U{RectBivariateSpline
1280 <https://docs.SciPy.org/doc/scipy/reference/generated/
1281 scipy.interpolate.RectBivariateSpline.html>}
1282 only (C{int}).
1284 @raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid B{C{crop}},
1285 B{C{kind}} or B{C{smooth}}.
1287 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.
1289 @raise LenError: EGM dataset B{C{egm_pgm}} axis mismatch.
1291 @raise SciPyError: A C{RectBivariateSpline} or C{inter2d} issue.
1293 @raise SciPyWarning: A C{RectBivariateSpline} or C{inter2d}
1294 warning as exception.
1296 @raise TypeError: Invalid B{C{datum}}.
1298 @note: The U{GeographicLib egm*.pgm<https://GeographicLib.SourceForge.io/
1299 html/geoid.html#geoidinst>} file sizes are based on a 2-byte
1300 C{int} height converted to 8-byte C{dtype float64} for C{scipy}
1301 interpolators. Therefore, internal memory usage is 4 times the
1302 C{egm*.pgm} file size and may exceed the available memory,
1303 especially with 32-bit Python. To reduce memory usage, set
1304 keyword argument B{C{crop}} to the region of interest. For example
1305 C{B{crop}=(20, -125, 50, -65)} covers the U{conterminous US<https://
1306 www.NGS.NOAA.gov/GEOID/GEOID12B/maps/GEOID12B_CONUS_grids.png>}
1307 (CONUS), less than 3% of the entire C{egm2008-1.pgm} dataset.
1309 @see: Class L{GeoidKarney} and function L{egmGeoidHeights}.
1310 '''
1311 np = self.numpy
1312 self._u2B = np.dtype(self.endian).itemsize
1314 g = self._open(egm_pgm, datum, kind, name, smooth)
1315 self._pgm = p = _PGM(g, pgm=egm_pgm, itemsize=self.u2B, sizeB=self.sizeB)
1316 if crop:
1317 g = p._cropped(g, abs(kind) + 1, *self._swne(crop))
1318 self._g2ll2 = self._g2ll2_cropped
1319 self._ll2g2 = self._ll2g2_cropped
1320 if map2(int, np.__version__.split(_DOT_)[:2]) < (1, 9):
1321 g = open(g.name, _rb_) # reopen tempfile for numpy 1.8.0-
1322 self._cropped = True
1323 else:
1324 self._cropped = False
1326 try:
1327 # U{numpy dtype formats are different from Python struct formats
1328 # <https://docs.SciPy.org/doc/numpy-1.15.0/reference/arrays.dtypes.html>}
1329 # read all heights, skipping the PGM header lines, converted to float
1330 hs = self._load(g, self.endian, p.knots, p.skip).reshape(p.nlat, p.nlon) * p.Scale
1331 if p.Offset: # offset
1332 hs = p.Offset + hs
1333 if p.dlat < 0: # flip the rows
1334 hs = np.flipud(hs)
1335 _GeoidBase.__init__(self, hs, p)
1336 except Exception as x:
1337 raise _SciPyIssue(x, _in_, repr(egm_pgm))
1338 finally:
1339 g.close()
1341 def _g2ll2(self, lat, lon):
1342 # convert grid (lat, lon) to earth (lat, lon), uncropped
1343 while lon > _180_0:
1344 lon -= _360_0
1345 return lat, lon
1347 def _g2ll2_cropped(self, lat, lon):
1348 # convert cropped grid (lat, lon) to earth (lat, lon)
1349 return lat, lon - self._lon_og
1351 def _ll2g2(self, lat, lon):
1352 # convert earth (lat, lon) to grid (lat, lon), uncropped
1353 while lon < 0:
1354 lon += _360_0
1355 return lat, lon
1357 def _ll2g2_cropped(self, lat, lon):
1358 # convert earth (lat, lon) to cropped grid (lat, lon)
1359 return lat, lon + self._lon_of
1361 if _FOR_DOCS:
1362 __call__ = _GeoidBase.__call__
1363 height = _GeoidBase.height
1366class _Gpars(_Named):
1367 '''(INTERNAL) Basic geoid parameters.
1368 '''
1369 # interpolator parameters
1370 dlat = 0 # +/- latitude resolution in C{degrees}
1371 dlon = 0 # longitude resolution in C{degrees}
1372 nlat = 1 # number of latitude knots (C{int})
1373 nlon = 0 # number of longitude knots (C{int})
1374 rlat = 0 # +/- latitude resolution in C{float}, 1 / .dlat
1375 rlon = 0 # longitude resolution in C{float}, 1 / .dlon
1376 slat = 0 # nothern- or southern most latitude (C{degrees90})
1377 wlon = 0 # western-most longitude in Eastern lon (C{degrees360})
1379 flon = 0 # forward, earth to grid longitude offset
1380 glon = 0 # reverse, grid to earth longitude offset
1382 knots = 0 # number of knots, nlat * nlon (C{int})
1383 skip = 0 # header bytes to skip (C{int})
1385 def __repr__(self):
1386 t = _COMMASPACE_.join(pairs((a, getattr(self, a)) for
1387 a in dir(self.__class__)
1388 if a[:1].isupper()))
1389 return _COLONSPACE_(self, t)
1391 def __str__(self):
1392 return Fmt.PAREN(self.classname, repr(self.name))
1395class _PGM(_Gpars):
1396 '''(INTERNAL) Parse an C{egm*.pgm} geoid dataset file.
1398 # Geoid file in PGM format for the GeographicLib::Geoid class
1399 # Description WGS84 EGM96, 5-minute grid
1400 # URL https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/egm96/egm96.html
1401 # DateTime 2009-08-29 18:45:03
1402 # MaxBilinearError 0.140
1403 # RMSBilinearError 0.005
1404 # MaxCubicError 0.003
1405 # RMSCubicError 0.001
1406 # Offset -108
1407 # Scale 0.003
1408 # Origin 90N 0E
1409 # AREA_OR_POINT Point
1410 # Vertical_Datum WGS84
1411 <width> <height>
1412 <pixel>
1413 ...
1414 '''
1415 crop4 = () # 4-tuple (C{south, west, north, east}).
1416 egm = None
1417 glon = 180 # reverse offset, uncropped
1418# pgm = NN # name
1419 sizeB = 0
1420 u2B = 2 # item size of grid height (C{int}).
1422 @staticmethod
1423 def _llstr2floats(latlon):
1424 # llstr to (lat, lon) floats
1425 lat, lon = latlon.split()
1426 return parseDMS2(lat, lon)
1428 # PGM file attributes, CamelCase but not .istitle()
1429 AREA_OR_POINT = str
1430 DateTime = str
1431 Description = str # 'WGS84 EGM96, 5-minute grid'
1432 Geoid = str # 'file in PGM format for the GeographicLib::Geoid class'
1433 MaxBilinearError = float
1434 MaxCubicError = float
1435 Offset = float
1436 Origin = _llstr2floats
1437 Pixel = 0
1438 RMSBilinearError = float
1439 RMSCubicError = float
1440 Scale = float
1441 URL = str # 'https://Earth-Info.NGA.mil/GandG/wgs84/...'
1442 Vertical_Datum = str
1444 def __init__(self, g, pgm=NN, itemsize=0, sizeB=0): # MCCABE 22
1445 '''(INTERNAL) New C{_PGM} parsed C{egm*.pgm} geoid dataset.
1446 '''
1447 self.name = pgm # geoid file name
1448 if itemsize:
1449 self._u2B = itemsize
1450 if sizeB:
1451 self.sizeB = sizeB
1453 t = g.readline() # make sure newline == '\n'
1454 if t != b'P5\n' and t.strip() != b'P5':
1455 raise self._Errorf(_format_, _header_, t)
1457 while True: # read all # Attr ... lines,
1458 try: # ignore empty ones or comments
1459 t = g.readline().strip()
1460 if t.startswith(_bHASH_):
1461 t = t.lstrip(_bHASH_).lstrip()
1462 a, v = map(_ub2str, t.split(None, 1))
1463 f = getattr(_PGM, a, None)
1464 if callable(f) and a[:1].isupper():
1465 setattr(self, a, f(v))
1466 elif t:
1467 break
1468 except (TypeError, ValueError):
1469 raise self._Errorf(_format_, 'Attr', t)
1470 else: # should never get here
1471 raise self._Errorf(_format_, _header_, g.tell())
1473 try: # must be (even) width and (odd) height
1474 nlon, nlat = map(int, t.split())
1475 if nlon < 2 or nlon > (360 * 60) or isodd(nlon) or \
1476 nlat < 2 or nlat > (181 * 60) or not isodd(nlat):
1477 raise ValueError
1478 except (TypeError, ValueError):
1479 raise self._Errorf(_format_, _SPACE_(_width_, _height_), t)
1481 try: # must be 16 bit pixel height
1482 t = g.readline().strip()
1483 self.Pixel = int(t)
1484 if not 255 < self.Pixel < 65536: # >u2 or >H only
1485 raise ValueError
1486 except (TypeError, ValueError):
1487 raise self._Errorf(_format_, 'pixel', t)
1489 for a in dir(_PGM): # set undefined # Attr ... to None
1490 if a[:1].isupper() and callable(getattr(self, a)):
1491 setattr(self, a, None)
1493 if self.Origin is None:
1494 raise self._Errorf(_format_, 'Origin', self.Origin)
1495 if self.Offset is None or self.Offset > 0:
1496 raise self._Errorf(_format_, 'Offset', self.Offset)
1497 if self.Scale is None or self.Scale < EPS:
1498 raise self._Errorf(_format_, 'Scale', self.Scale)
1500 self.skip = g.tell()
1501 self.knots = nlat * nlon
1503 self.nlat, self.nlon = nlat, nlon
1504 self.slat, self.wlon = self.Origin
1505 # note, negative .dlat and .rlat since rows
1506 # are from .slat 90N down in decreasing lat
1507 self.dlat, self.dlon = _180_0 / (1 - nlat), _360_0 / nlon
1508 self.rlat, self.rlon = (1 - nlat) / _180_0, nlon / _360_0
1510 # grid corners in earth (lat, lon), .slat = 90, .dlat < 0
1511 n = float(self.slat)
1512 s = n + self.dlat * (nlat - 1)
1513 w = self.wlon - self.glon
1514 e = w + self.dlon * nlon
1515 self.crop4 = s, w, n, e
1517 n = self.sizeB - self.skip
1518 if n > 0 and n != (self.knots * self.u2B):
1519 raise self._Errorf('%s(%s x %s != %s)', _assert_, nlat, nlon, n)
1521 def _cropped(self, g, k1, south, west, north, east): # MCCABE 15
1522 '''Crop the geoid to (south, west, north, east) box.
1523 '''
1524 # flon offset for both west and east
1525 f = 360 if west < 0 else 0
1526 # earth (lat, lon) to grid indices (y, x),
1527 # note y is decreasing, i.e. n < s
1528 s, w = self._lle2yx2(south, west, f)
1529 n, e = self._lle2yx2(north, east, f)
1530 s += 1 # s > n
1531 e += 1 # e > w
1533 hi, wi = self.nlat, self.nlon
1534 # handle special cases
1535 if (s - n) > hi:
1536 n, s = 0, hi # entire lat range
1537 if (e - w) > wi:
1538 w, e, f = 0, wi, 180 # entire lon range
1539 if s == hi and w == n == 0 and e == wi:
1540 return g # use entire geoid as-is
1542 if (e - w) < k1 or (s - n) < (k1 + 1):
1543 raise self._Errorf(_format_, 'swne', (north - south, east - west))
1545 if e > wi > w: # wrap around
1546 # read w..wi and 0..e
1547 r, p = (wi - w), (e - wi)
1548 elif e > w:
1549 r, p = (e - w), 0
1550 else:
1551 raise self._Errorf('%s(%s < %s)', _assert_, w, e)
1553 # convert to bytes
1554 r *= self.u2B
1555 p *= self.u2B
1556 q = wi * self.u2B # stride
1557 # number of rows and cols to skip from
1558 # the original (.slat, .wlon) origin
1559 z = self.skip + (n * wi + w) * self.u2B
1560 # sanity check
1561 if r < 2 or p < 0 or q < 2 or z < self.skip \
1562 or z > self.sizeB:
1563 raise self._Errorf(_format_, _assert_, (r, p, q, z))
1565 # can't use _BytesIO since numpy
1566 # needs .fileno attr in .fromfile
1567 t, c = 0, self._tmpfile()
1568 # reading (s - n) rows, forward
1569 for y in range(n, s): # PYCHOK y unused
1570 g.seek(z, _SEEK_SET)
1571 # Python 2 tmpfile.write returns None
1572 t += c.write(g.read(r)) or r
1573 if p: # wrap around to start of row
1574 g.seek(-q, _SEEK_CUR)
1575 # assert(g.tell() == (z - w * self.u2B))
1576 # Python 2 tmpfile.write returns None
1577 t += c.write(g.read(p)) or p
1578 z += q
1579 c.flush()
1580 g.close()
1582 s -= n # nlat
1583 e -= w # nlon
1584 k = s * e # knots
1585 z = k * self.u2B
1586 if t != z:
1587 raise self._Errorf('%s(%s != %s) %s', _assert_, t, z, self)
1589 # update the _Gpars accordingly, note attributes
1590 # .dlat, .dlon, .rlat and .rlon remain unchanged
1591 self.slat += n * self.dlat
1592 self.wlon += w * self.dlon
1593 self.nlat = s
1594 self.nlon = e
1595 self.flon = self.glon = f
1597 self.crop4 = south, west, north, east
1598 self.knots = k
1599 self.skip = 0 # no header lines in c
1601 c.seek(0, _SEEK_SET)
1602 # c = open(c.name, _rb_) # reopen for numpy 1.8.0-
1603 return c
1605 def _Errorf(self, fmt, *args): # PYCHOK no cover
1606 t = fmt % args
1607 e = self.pgm or NN
1608 if e:
1609 t = _SPACE_(t, _in_, repr(e))
1610 return PGMError(t)
1612 def _lle2yx2(self, lat, lon, flon):
1613 # earth (lat, lon) to grid indices (y, x)
1614 # with .dlat decreasing from 90N .slat
1615 lat -= self.slat
1616 lon += flon - self.wlon
1617 return (min(self.nlat - 1, max(0, int(lat * self.rlat))),
1618 max(0, int(lon * self.rlon)))
1620 def _tmpfile(self):
1621 # create a tmpfile to hold the cropped geoid grid
1622 try:
1623 from tempfile import NamedTemporaryFile as tmpfile
1624 except ImportError: # Python 2.7.16-
1625 from os import tmpfile
1626 t = _os_path.splitext(_os_path.basename(self.pgm))[0]
1627 f = tmpfile(mode='w+b', prefix=t or 'egm')
1628 f.seek(0, _SEEK_SET) # force overwrite
1629 return f
1631 @Property_RO
1632 def pgm(self):
1633 '''Get the geoid file name (C{str}).
1634 '''
1635 return self.name
1638class PGMError(GeoidError):
1639 '''Issue parsing or cropping an C{egm*.pgm} geoid dataset.
1640 '''
1641 pass
1644def egmGeoidHeights(GeoidHeights_dat):
1645 '''Generate geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/
1646 html/geoid.html#geoidinst>} height tests from U{GeoidHeights.dat
1647 <https://SourceForge.net/projects/geographiclib/files/testdata/>}
1648 U{Test data for Geoids<https://GeographicLib.SourceForge.io/C++/doc/
1649 geoid.html#testgeoid>}.
1651 @arg GeoidHeights_dat: The un-gz-ed C{GeoidHeights.dat} file
1652 (C{str} or C{file} handle).
1654 @return: For each test, yield a L{GeoidHeight5Tuple}C{(lat, lon,
1655 egm84, egm96, egm2008)}.
1657 @raise GeoidError: Invalid B{C{GeoidHeights_dat}}.
1659 @note: Function L{egmGeoidHeights} is used to test the geoids
1660 L{GeoidKarney} and L{GeoidPGM}, see PyGeodesy module
1661 C{test/testGeoids.py}.
1662 '''
1663 dat = GeoidHeights_dat
1664 if isinstance(dat, bytes):
1665 dat = _BytesIO(dat)
1667 try:
1668 dat.seek(0, _SEEK_SET) # reset
1669 except AttributeError as x:
1670 raise GeoidError(GeoidHeights_dat=type(dat), cause=x)
1672 for t in dat.readlines():
1673 t = t.strip()
1674 if t and not t.startswith(_bHASH_):
1675 lat, lon, egm84, egm96, egm2008 = map(float, t.split())
1676 while lon > _180_0: # EasternLon to earth lon
1677 lon -= _360_0
1678 yield GeoidHeight5Tuple(lat, lon, egm84, egm96, egm2008)
1681__all__ += _ALL_DOCS(_GeoidBase)
1683if __name__ == '__main__':
1685 from pygeodesy.lazily import printf, _sys
1687 _crop = ()
1688 _GeoidEGM = GeoidKarney
1689 _kind = 3
1691 geoids = _sys.argv[1:]
1692 while geoids:
1693 geoid = geoids.pop(0)
1695 if '-crop'.startswith(geoid.lower()):
1696 _crop = 20, -125, 50, -65 # CONUS
1698 elif '-karney'.startswith(geoid.lower()):
1699 _GeoidEGM = GeoidKarney
1701 elif '-kind'.startswith(geoid.lower()):
1702 _kind = int(geoids.pop(0))
1704 elif '-pgm'.startswith(geoid.lower()):
1705 _GeoidEGM = GeoidPGM
1707 elif geoid[-4:].lower() in ('.pgm',):
1708 g = _GeoidEGM(geoid, crop=_crop, kind=_kind)
1709 printf(g.toStr(), nt=1, nl=1)
1710 printf(repr(g.pgm), nt=1)
1711 # <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>:
1712 # The height of the EGM96 geoid at Timbuktu
1713 # echo 16:46:33N 3:00:34W | GeoidEval
1714 # => 28.7068 -0.02e-6 -1.73e-6
1715 # The 1st number is the height of the geoid, the 2nd and
1716 # 3rd are its slopes in northerly and easterly direction
1717 t = 'Timbuktu %s' % (g,)
1718 k = {'egm84-15.pgm': '31.2979',
1719 'egm96-5.pgm': '28.7067',
1720 'egm2008-1.pgm': '28.7880'}.get(g.name.lower(), '28.7880')
1721 ll = parseDMS2('16:46:33N', '3:00:34W', sep=':')
1722 for ll in (ll, (16.776, -3.009),):
1723 try:
1724 h, ll = g.height(*ll), fstr(ll, prec=6)
1725 printf('%s.height(%s): %.4F vs %s', t, ll, h, k)
1726 except (GeoidError, RangeError) as x:
1727 printf(_COLONSPACE_(t, str(x)))
1729 elif geoid[-4:].lower() in ('.bin',):
1730 g = GeoidG2012B(geoid, kind=_kind)
1731 printf(g.toStr())
1733 else:
1734 raise GeoidError(grid=repr(geoid))
1736_I = int # PYCHOK unused _I
1737del _intCs # trash ints cache
1739# **) MIT License
1740#
1741# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1742#
1743# Permission is hereby granted, free of charge, to any person obtaining a
1744# copy of this software and associated documentation files (the "Software"),
1745# to deal in the Software without restriction, including without limitation
1746# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1747# and/or sell copies of the Software, and to permit persons to whom the
1748# Software is furnished to do so, subject to the following conditions:
1749#
1750# The above copyright notice and this permission notice shall be included
1751# in all copies or substantial portions of the Software.
1752#
1753# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1754# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1755# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1756# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1757# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1758# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1759# OTHER DEALINGS IN THE SOFTWARE.
1761# <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>
1762# _lowerleft = -90, -179, -30.1500 # egm2008-1.pgm
1763# _lowerleft = -90, -179, -29.5350 # egm96-5.pgm
1764# _lowerleft = -90, -179, -29.7120 # egm84-15.pgm
1766# _center = 0, 0, 17.2260 # egm2008-1.pgm
1767# _center = 0, 0, 17.1630 # egm96-5.pgm
1768# _center = 0, 0, 18.3296 # egm84-15.pgm
1770# _upperright = 90, 180, 14.8980 # egm2008-1.pgm
1771# _upperright = 90, 180, 13.6050 # egm96-5.pgm
1772# _upperright = 90, 180, 13.0980 # egm84-15.pgm
1775# % python3 -m pygeodesy.geoids [-Karney] ../testGeoids/egm*.pgm
1776#
1777# 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)
1778#
1779# _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'
1780#
1781# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880
1782# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880
1783#
1784# 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)
1785#
1786# _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'
1787#
1788# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979
1789# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979
1790#
1791# 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)
1792#
1793# _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'
1794#
1795# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067
1796# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067
1799# % python3 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm
1800#
1801# 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)
1802#
1803# _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'
1804#
1805# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880
1806# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880
1807#
1808# 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)
1809#
1810# _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'
1811#
1812# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979
1813# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979
1814#
1815# 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)
1816#
1817# _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'
1818#
1819# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067
1820# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067
1823# % python2 -m pygeodesy.geoids -PGM ../testGeoids/egm*.pgm
1824#
1825# 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)
1826#
1827# _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'
1828#
1829# Timbuktu GeoidPGM('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880
1830# Timbuktu GeoidPGM('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880
1831#
1832# 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)
1833#
1834# _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'
1835#
1836# Timbuktu GeoidPGM('egm84-15.pgm').height(16.775833, -3.009444): 31.2979 vs 31.2979
1837# Timbuktu GeoidPGM('egm84-15.pgm').height(16.776, -3.009): 31.2975 vs 31.2979
1838#
1839# 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)
1840#
1841# _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'
1842#
1843# Timbuktu GeoidPGM('egm96-5.pgm').height(16.775833, -3.009444): 28.7065 vs 28.7067
1844# Timbuktu GeoidPGM('egm96-5.pgm').height(16.776, -3.009): 28.7064 vs 28.7067