Coverage for pygeodesy/geoids.py : 95%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
Classes L{GeoidG2012B}, L{GeoidKarney} and L{GeoidPGM} to interpolate the height of various U{geoid<https://WikiPedia.org/wiki/Geoid>}s at C{LatLon} locations or separate lat-/longitudes using different interpolation methods and C{geoid} model files.
L{GeoidKarney} is a transcoding of I{Charles Karney}'s C++ class U{Geoid <https://GeographicLib.SourceForge.io/html/geoid.html>} to pure Python. The L{GeoidG2012B} and L{GeoidPGM} interpolators both depend on U{scipy <https://SciPy.org>} and U{numpy<https://PyPI.org/project/numpy>} and require those packages to be installed.
In addition, each geoid interpolator needs C{grid knots} (down)loaded from a C{geoid} model file, specific to the interpolator, more details below. For each interpolator, there are several interpolation choices, for example I{linear}, I{cubic}, etc.
B{Typical usage} is as follows. First, create an interpolator from a C{geoid} model file, containing locations with known heights also referred to as the C{grid knots}.
C{>>> ginterpolator = GeoidXyz(geoid_file, **options)}
Then, get the interpolated geoid height of one or several C{LatLon} location(s) with
C{>>> h = ginterpolator(ll)}
or
C{>>> h0, h1, h2, ... = ginterpolator(ll0, ll1, ll2, ...)}
or
C{>>> hs = ginterpolator(lls) # list, tuple, generator, ...}
For separate lat- and longitudes invoke the C{.height} method as
C{>>> h = ginterpolator.height(lat, lon)}
or as
C{>>> h0, h1, h2, ... = ginterpolator.height(lats, lons) # lists, tuples, ...}
Errors from C{scipy} are raised as L{SciPyError}s. Warnings issued by C{scipy} can be thrown as L{SciPyWarning} exceptions, provided Python C{warnings} are filtered accordingly, see L{SciPyWarning}.
@see: I{Karney}'s U{GeographicLib<https://GeographicLib.SourceForge.io/ html/index.html>}, U{Geoid height<https://GeographicLib.SourceForge.io/ html/geoid.html>} and U{Installing the Geoid datasets<https:// GeographicLib.SourceForge.io/html/geoid.html#geoidinst>}, U{SciPy <https://docs.SciPy.org/doc/scipy/reference/interpolate.html>} interpolation U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/ reference/generated/scipy.interpolate.RectBivariateSpline.html>} and U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/ scipy.interpolate.interp2d.html>} and the functions L{elevations.elevation2} and L{elevations.geoidHeight2}. ''' # make sure int/int division yields float quotient, see .basics
_E_, _float as _F, _height_, _in_, _kind_, \ _knots_, _lat_, _linear_, _lon_, _mean_, _N_, \ _n_a_, _not_, _numpy_, _on_, _outside_, _S_, _s_, \ _scipy_, _SPACE_, _stdev_, _supported_, _tbd_, \ _W_, _width_, _4_, _0_0, _1_0, _180_0, _360_0
# temporarily hold a single instance for each int value -3: _cubic_, -5: 'quintic'}
'''(INTERNAL) Base class for C{Geoid...}s. ''' # _name = NN # _Named
'''(INTERNAL) Set up the grid axes, the C{SciPy} interpolator and several internal geoid attributes.
@arg hs: Grid knots with known height (C{numpy 2darray}). @arg p: The C{slat, wlon, nlat, nlon, dlat, dlon} and other geoid parameters (C{INTERNAL}).
@raise GeoidError: Incompatible grid B{C{hs}} shape or invalid B{C{kind}}.
@raise LenError: Mismatch grid B{C{hs}} axis.
@raise SciPyError: A C{scipy.interpolate.inter2d} or C{-.RectBivariateSpline} issue.
@raise SciPyWarning: A C{scipy.interpolate.inter2d} or C{-.RectBivariateSpline} warning as exception. ''' # for 2d scipy.interpolate.interp2d(xs, ys, hs, ...) and # scipy.interpolate.RectBivariateSpline(ys, xs, hs, ...) # require the shape of hs to be (len(ys), len(xs)), note # the different (xs, ys, ...) and (ys, xs, ...) orders raise GeoidError(shape=hs.shape, txt=_incompatible((p.nlat, p.nlon)))
# both axes and bounding box
# geoid grids are typically stored in row-major order, some # with rows (90..-90) reversed and columns (0..360) wrapped # to Easten longitude, 0 <= east < 180 and 180 <= west < 0 self._interp2d = self._spi.interp2d(xs, ys, hs, kind=_interp2d_ks[k]) ky=k, kx=k, s=self._smooth).ev else: raise GeoidError(kind=k)
# shrink the box by 1 unit on every side # bb += self._lat_d, -self._lat_d, self._lon_d, -self._lon_d
'''Interpolate the geoid height for one or several locations.
@arg llis: The location or locations (C{LatLon}, ... or C{LatLon}s).
@return: A single interpolated geoid height (C{float}) or a list or tuple of interpolated geoid heights (C{float}s).
@raise GeoidError: Insufficient number of B{C{llis}}, an invalid B{C{lli}} or the C{egm*.pgm} geoid file is closed.
@raise RangeError: An B{C{lli}} is outside this geoid's lat- or longitude range.
@raise SciPyError: A C{scipy.interpolate.inter2d} or C{-.RectBivariateSpline} issue.
@raise SciPyWarning: A C{scipy.interpolate.inter2d} or C{-.RectBivariateSpline} warning as exception. '''
'''Open context. '''
'''Close context. ''' # return None # XXX False
return self.toStr()
# handle __call__
# XXX avoid str(LatLon()) degree symbols raise type(x)(t, lli, txt=str(x)) except Exception as x: # PYCHOK no cover if scipy and self.scipy: raise _SciPyIssue(x) else: raise
'''Close the C{egm*.pgm} geoid file if open (and applicable). '''
'''Get the C{egm*.pgm} geoid file status. '''
# only used for .interpolate.interp2d, but # overwritten for .RectBivariateSpline, # note (y, x) must be flipped! return self._interp2d(x, y)
# build grid axis, hi = lo + (n - 1) * d raise LenError(self.__class__, grid=m, **{name: n}) i = Fmt.SQUARE(name, i) raise GeoidError(i, e, txt=_non_increasing_)
def _g2ll2(self, lat, lon): # PYCHOK no cover notOverloaded(self, lat, lon)
# convert grid (y, x) indices to grid (lat, lon) (self._lon_lo + self._lon_of + self._lon_d * x))
raise RangeError(lli=lli, txt=_SPACE_(_outside_, _on_, out))
def _ll2g2(self, lat, lon): # PYCHOK no cover notOverloaded(self, lat, lon)
# <https://docs.SciPy.org/doc/numpy/reference/generated/ # numpy.argmin.html#numpy.argmin>
# numpy.fromfile, like .frombuffer
# open the geoid file except (IOError, OSError) as x: raise GeoidError(geoid=geoid, txt=str(x))
self._datum = _ellipsoidal_datum(datum, name=name) _HeightBase.name.fset(self, name) # rename self._smooth = Int_(smooth=smooth, Error=GeoidError, low=0)
# crop box to 4-tuple (s, w, n, e) crop[1].lat, crop[1].lon) else: # (s, w, n, e) -180 <= w <= (e - _1_0) <= 179: except (IndexError, TypeError, ValueError): # PYCHOK no cover pass raise GeoidError(crop=crop)
''' Cache for L{center}. ''' favg(self._lon_lo, self._lon_hi))
'''Return the center location and height of this geoid.
@kwarg LatLon: Optional class to return the location and height (C{LatLon}) or C{None}.
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the grid center location. '''
'''Is geoid cropped (C{bool} or C{None} if crop not supported). '''
'''Get the grid C{scipy} U{dtype<https://docs.SciPy.org/doc/numpy/ reference/generated/numpy.ndarray.dtype.html>} (C{numpy.dtype}). '''
'''Get the geoid endianess and U{dtype<https://docs.SciPy.org/ doc/numpy/reference/generated/numpy.dtype.html>} (C{str}). '''
'''Interpolate the geoid height for one or several lat-/longitudes.
@arg lats: Latitude or latitudes (C{degrees} or C{degrees}s). @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
@return: A single interpolated geoid height (C{float}) or a list of interpolated geoid heights (C{float}s).
@raise GeoidError: Insufficient or non-matching number of B{C{lats}} and B{C{lons}}.
@raise RangeError: A B{C{lat}} or B{C{lon}} is outside this geoid's lat- or longitude range.
@raise SciPyError: A C{scipy.interpolate.inter2d} or C{-.RectBivariateSpline} issue.
@raise SciPyWarning: A C{scipy.interpolate.inter2d} or C{-.RectBivariateSpline} warning as exception. '''
'''(INTERNAL) Cache for L{highest}. '''
'''Return the location and largest height of this geoid.
@kwarg LatLon: Optional class to return the location and height (C{LatLon}) or C{None}.
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the highest grid location. '''
'''Get the number of cache hits (C{int} or C{None}). '''
'''Get the interpolator kind and order (C{int}). '''
'''Get the number of grid knots (C{int}). '''
'''(INTERNAL) Cache for L{lowerleft}. '''
'''Return the lower-left location and height of this geoid.
@kwarg LatLon: Optional class to return the location (C{LatLon}) and height or C{None}.
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the lower-left, SW grid corner. '''
'''(INTERNAL) Cache for L{loweright}. '''
'''Return the lower-right location and height of this geoid.
@kwarg LatLon: Optional class to return the location and height (C{LatLon}) or C{None}.
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the lower-right, SE grid corner. '''
'''(INTERNAL) Cache for L{lowest}. '''
'''Return the location and lowest height of this geoid.
@kwarg LatLon: Optional class to return the location and height (C{LatLon}) or C{None}.
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the lowest grid location. '''
'''Get the mean of this geoid's heights (C{float}). '''
'''Get the name of this geoid (C{str}). '''
'''Get the grid in-memory size in bytes (C{int}). '''
'''Get the imported C{numpy} version (C{str}). '''
'''Check whether a location is outside this geoid's lat-/longitude or crop range.
@arg lat: The latitude (C{degrees}). @arg lon: The longitude (C{degrees}).
@return: A 1- or 2-character C{str} if outside or an empty C{str} if inside. ''' (_N_ if lat > self._lat_hi else NN)) + \ (_W_ if lon < self._lon_lo else (_E_ if lon > self._lon_hi else NN))
'''Get the PGM attributes (C{_PGM} or C{None} if not available/applicable). ''' return self._pgm
'''Get the imported C{scipy} version (C{str}). '''
'''Get the geoid grid file size in bytes (C{int}). '''
'''Get the C{RectBivariateSpline} smoothing (C{int}). '''
'''Get the standard deviation of this geoid's heights (C{float}) or C{None}. '''
'''This geoid and all geoid attributes as a string.
@kwarg prec: Optional number of decimal digits (0..9 or C{None} for default). Trailing zero decimals are stripped for B{C{prec}} values of 1 and above, but kept for negative B{C{prec}} values. @kwarg sep: Optional separator (C{str}).
@return: Geoid name and attributes (C{str}). ''' (self.lowerleft, self.upperright, self.center, self.highest, self.lowest)) + \ attrs( _mean_, _stdev_, prec=prec, Nones=False) + \ attrs((_kind_, 'smooth')[:s], prec=prec, Nones=False) + \ attrs( 'cropped', 'dtype', _endian_, 'hits', _knots_, 'nBytes', 'sizeB', _scipy_, _numpy_, prec=prec, Nones=False)
'''(INTERNAL) Cache for L{upperleft}. '''
'''Return the upper-left location and height of this geoid.
@kwarg LatLon: Optional class to return the location and height (C{LatLon}) or C{None}.
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the upper-left, NW grid corner. '''
'''(INTERNAL) Cache for L{upperright}. '''
'''Return the upper-right location and height of this geoid.
@kwarg LatLon: Optional class to return the location and height (C{LatLon}) or C{None}.
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the upper-right, NE grid corner. '''
'''Geoid interpolator C{Geoid...} or interpolation issue. '''
'''Geoid height interpolator for U{GEOID12B Model <https://www.NGS.NOAA.gov/GEOID/GEOID12B/>} grids U{CONUS <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_CONUS.shtml>}, U{Alaska<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_AK.shtml>}, U{Hawaii<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_HI.shtml>}, U{Guam and Northern Mariana Islands <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_GMNI.shtml>}, U{Puerto Rico and U.S. Virgin Islands <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_PRVI.shtml>} and U{American Samoa<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_AS.shtml>} based on C{SciPy} U{RectBivariateSpline<https://docs.SciPy.org/doc/ scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>} or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/ scipy.interpolate.interp2d.html>} interpolation.
Use any of the binary C{le} (little endian) or C{be} (big endian) C{g2012b*.bin} grid files. ''' kind=3, name=NN, smooth=0): '''New L{GeoidG2012B} interpolator.
@arg g2012b_bin: A C{GEOID12B} grid file name (C{.bin}). @kwarg crop: Optional crop box, not supported (C{None}). @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}), default C{WGS84}. @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/ reference/generated/scipy.interpolate.RectBivariateSpline.html>}, -2 for U{interp2d linear<https://docs.SciPy.org/doc/scipy/ reference/generated/scipy.interpolate.interp2d.html>}, -3 for C{interp2d cubic} or -5 for C{interp2d quintic}. @kwarg name: Optional geoid name (C{str}). @kwarg smooth: Smoothing factor for U{RectBivariateSpline <https://docs.SciPy.org/doc/scipy/reference/generated/ scipy.interpolate.RectBivariateSpline.html>} only (C{int}).
@raise GeoidError: G2012B grid file B{C{g2012b_bin}} issue or invalid B{C{crop}}, B{C{kind}} or B{C{smooth}}.
@raise ImportError: Package C{numpy} or C{scipy} not found or not installed.
@raise LenError: Grid file B{C{g2012b_bin}} axis mismatch.
@raise SciPyError: A C{RectBivariateSpline} or C{inter2d} issue.
@raise SciPyWarning: A C{RectBivariateSpline} or C{inter2d} warning as exception.
@raise TypeError: Invalid B{C{datum}}. ''' raise GeoidError(crop=crop, txt=_not_(_supported_))
# U{numpy dtype formats are different from Python struct formats # <https://docs.SciPy.org/doc/numpy-1.15.0/reference/arrays.dtypes.html>} # skip 4xf8, get 3xi4 and 1 < p.nlon < n: else: # couldn't validate endian raise GeoidError(_endian_)
# get the first 4xf8 # read all f4 heights, ignoring the first 4xf8 and 3xi4
except Exception as x: # PYCHOK no cover raise _SciPyIssue(x, _in_, repr(g2012b_bin)) finally:
# convert grid (lat, lon) to earth (lat, lon)
# convert earth (lat, lon) to grid (lat, lon)
if _FOR_DOCS: __call__ = _GeoidBase.__call__ height = _GeoidBase.height
'''5-Tuple C{(lat, lon, egm84, egm96, egm2008)} for U{GeoidHeights.dat <https://SourceForge.net/projects/geographiclib/files/testdata/>} tests with the heights for 3 different EGM grids at C{degrees90} and C{degrees180} degrees (after converting C{lon} from original C{0 <= EasterLon <= 360}). '''
'''(INTERNAL) Cache a single C{int} constant. '''
'''(INTERNAL) Cache a tuple of single C{int} constants. '''
'''Geoid height interpolator for I{Karney}'s U{GeographicLib Earth Gravitational Model (EGM)<https://GeographicLib.SourceForge.io/html/ geoid.html>} geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/ html/geoid.html#geoidinst>} datasets using bilinear or U{cubic <https://dl.ACM.org/citation.cfm?id=368443>} interpolation and U{caching <https://GeographicLib.SourceForge.io/html/geoid.html#geoidcache>} in pure Python, transcoded from I{Karney}'s U{C++ class Geoid <https://GeographicLib.SourceForge.io/html/geoid.html#geoidinterp>}.
Use any of the geoid U{egm84-, egm96- or egm2008-*.pgm <https://GeographicLib.SourceForge.io/html/geoid.html#geoidinst>} datasets. ''' # matrices c3n_, c3, c3s_, transposed from GeographicLib/Geoid.cpp _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), _T(-131, 7, -31, -62, -62, -31, 45, 216, 156, -45, -55, -7), _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), _T(138, -138, 0, 0, 0, 0, -183, 33, 153, -3, 48, -48), # PYCHOK indent _T(144, 42, -62, -124, -124, -62, -9, 87, 99, 9, 42, -42), _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), _T(0, 0, 0, 0, 0, 0, 93, -93, -93, 93, 0, 0), _T(-102, 102, 0, 0, 0, 0, 18, 12, -12, -18, -84, 84), _T(-31, -31, 31, 62, 62, 31, 0, -93, -93, 0, 31, 31)), # PYCHOK indent
(_T(9, -9, 9, 186, 54, -9, -9, 54, -54, 9, -9, 9), _T(-18, 18, -88, -42, 162, -32, 8, -78, 78, -8, 18, -18), _T(-88, 8, -18, -42, -78, 18, 18, 162, 78, -18, -32, -8), _T(0, 0, 90, -150, 30, 30, 30, -90, 90, -30, 0, 0), _T(96, -96, 96, -96, -24, 24, -96, -24, 144, -24, 24, -24), # PYCHOK indent _T(90, 30, 0, -150, -90, 0, 0, 30, 90, 0, 30, -30), _T(0, 0, -20, 60, -60, 20, -20, 60, -60, 20, 0, 0), _T(0, 0, -60, 60, 60, -60, 60, -60, -60, 60, 0, 0), _T(-60, 60, 0, 60, -60, 0, 0, 60, -60, 0, -60, 60), _T(-20, -20, 0, 60, 60, 0, 0, -60, -60, 0, 20, 20)),
(_T(18, -18, 36, 210, 162, -36, 0, 0, 0, 0, -18, 18), # PYCHOK indent _T(-36, 36, -165, 45, 141, -21, 0, 0, 0, 0, 36, -36), _T(-122, -2, -27, -111, -75, 27, 62, 124, 124, 62, -64, 2), _T(0, 0, 93, -93, -93, 93, 0, 0, 0, 0, 0, 0), _T(120, -120, 147, -57, -129, 39, 0, 0, 0, 0, 66, -66), # PYCHOK indent _T(135, 51, -9, -192, -180, 9, 31, 62, 62, 31, 51, -51), _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), _T(0, 0, -93, 93, 93, -93, 0, 0, 0, 0, 0, 0), _T(-84, 84, 18, 12, -12, -18, 0, 0, 0, 0, -102, 102), _T(-31, -31, 0, 93, 93, 0, -31, -62, -62, -31, 31, 31)))
_T(1, 0), _T(0, 1), _T(1, 1))
_T( 1, -1), _T(-1, 0), _T( 0, 0), _T( 1, 0), _T( 2, 0), _T(-1, 1), _T( 0, 1), _T( 1, 1), _T( 2, 1), _T( 0, 2), _T( 1, 2))
# _highest = (-8.4, 147.367, 85.839) if egm2008-1.pgm else ( # (-8.167, 147.25, 85.422) if egm96-5.pgm else # (-4.5, 148.75, 81.33)) # egm84-15.pgm # _lowest = (4.7, 78.767, -106.911) if egm2008-1.pgm else ( # (4.667, 78.833, -107.043) if egm96-5.pgm else # (4.75, 79.25, -107.34)) # egm84-15.pgm
kind=3, name=NN, smooth=None): '''New L{GeoidKarney} interpolator.
@arg egm_pgm: An U{EGM geoid dataset<https://GeographicLib.SourceForge.io/ html/geoid.html#geoidinst>} file name (C{egm*.pgm}), see note below. @kwarg crop: Optional box to limit geoid locations, a 4-tuple (C{south, west, north, east}), 2-tuple (C{(south, west), (north, east)}) or 2, in C{degrees90} lat- and C{degrees180} longitudes or a 2-tuple (C{LatLonSW, LatLonNE}) of C{LatLon} instances. @kwarg datum: Optional grid datum (C{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}), default C{WGS84}. @kwarg kind: Interpolation order (C{int}), 2 for C{bilinear} or 3 for C{cubic}. @kwarg name: Optional geoid name (C{str}). @kwarg smooth: Smoothing factor, unsupported (C{None}).
@raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid B{C{crop}}, B{C{kind}} or B{C{smooth}}.
@raise TypeError: Invalid B{C{datum}}.
@see: Class L{GeoidPGM} and function L{egmGeoidHeights}.
@note: Geoid file B{C{egm_pgm}} remains open and must be closed by calling the C{close} method or by using this instance in a C{with B{GeoidKarney}(...) as ...} context. ''' raise GeoidError(smooth=smooth, txt=_not_(_supported_))
raise GeoidError(kind=kind)
# set earth (lat, lon) limits (s, w, n, e) self._lon_lo, \ self._lat_hi, \ self._lon_hi = self._swne(crop if crop else p.crop4)
'''Interpolate the geoid height for one or several locations.
@arg llis: The location or locations (C{LatLon}, ... or C{LatLon}s).
@return: A single interpolated geoid height (C{float}) or a list or tuple of interpolated geoid heights (C{float}s).
@raise GeoidError: Insufficient number of B{C{llis}}, an invalid B{C{lli}} or the C{egm*.pgm} geoid file is closed.
@raise RangeError: An B{C{lli}} is outside this geoid's lat- or longitude range. '''
# get the common denominator, the 10x12 cubic matrix and # the 12 cubic v-coefficients around geoid index (y, x) # read 4x4 ushorts, drop the 4 corners
else: # likely some wrapped y and/or x's
# interpolate the geoid height at grid (lat, lon)
# compute the bilinear 4-tuple and interpolate raw H else:
# compute the cubic 10-tuple and interpolate raw H else: # GeographicLib/Geoid.cpp Geoid::height(lat, lon) ... # real h = t[0] + fx * (t[1] + fx * (t[3] + fx * t[6])) + # fy * (t[2] + fx * (t[4] + fx * t[7]) + # fy * (t[5] + fx * t[8] + fy * t[9]));
# convert grid (lat, lon) to earth (lat, lon), uncropped
# convert grid (lat, lon) to grid (y, x) indices # note, slat = +90, rlat < 0 makes y >=0
# convert grid (y, x) indices to grid (lat, lon)
# convert earth lat(s) to min and max grid y indices
# convert earth (lat, lon) to grid (lat, lon), uncropped
# find highest or lowest, takes 10+ secs for egm2008-1.pgm geoid # (Python 2.7.16, macOS 10.13.6 High Sierra, iMac 3 GHz Core i3) else: # lowest
# get the ushort geoid height at geoid index (y, x), # like GeographicLib/Geoid.hpp real rawval(is, iy) x -= p.nlon else:
# get bilinear 4-tuple or 10x12 cubic matrix
# yield a 2-tuple (y, ushorts) for each row or for # the rows between two (or more) earth lat values # min(lat2) <= lat <= max(lat2) or 0 <= y < p.nlat
# position geoid to grid index (y, x) raise GeoidError('closed file: %r' % (p.egm,)) # IOError
'''Get the geoid's grid data type (C{str}). '''
'''Interpolate the geoid height for one or several lat-/longitudes.
@arg lats: Latitude or latitudes (C{degrees} or C{degrees}s). @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
@return: A single interpolated geoid height (C{float}) or a list of interpolated geoid heights (C{float}s).
@raise GeoidError: Insufficient or non-matching number of B{C{lats}} and B{C{lons}} or the C{egm*.pgm} geoid file is closed.
@raise RangeError: A B{C{lat}} or B{C{lon}} is outside this geoid's lat- or longitude range. '''
'''(INTERNAL) Cache for L{highest}. '''
'''Return the location and largest height of this geoid.
@kwarg LatLon: Optional class to return the location and height (C{LatLon}) or C{None}. @kwarg full: Search the full or limited latitude range (C{bool}).
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the highest grid location. '''
'''(INTERNAL) Cache for L{lowest}. '''
'''Return the location and lowest height of this geoid.
@kwarg LatLon: Optional class to return the location and height (C{LatLon}) or C{None}. @kwarg full: Search the full or limited latitude range (C{bool}).
@return: If B{C{LatLon}} is C{None}, a L{LatLon3Tuple}C{(lat, lon, height)} otherwise a B{C{LatLon}} instance with the lat-, longitude and height of the lowest grid location. '''
'''Get the PGM itemsize in bytes (C{int}). '''
'''Geoid height interpolator for I{Karney}'s U{GeographicLib Earth Gravitational Model (EGM)<https://GeographicLib.SourceForge.io/html/ geoid.html>} geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/ html/geoid.html#geoidinst>} datasets but based on C{SciPy} U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/reference/ generated/scipy.interpolate.RectBivariateSpline.html>} or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/ scipy.interpolate.interp2d.html>} interpolation.
Use any of the U{egm84-, egm96- or egm2008-*.pgm <https://GeographicLib.SourceForge.io/html/geoid.html#geoidinst>} datasets. However, unless cropped, an entire C{egm*.pgm} dataset is loaded into the C{SciPy} U{RectBivariateSpline<https://docs.SciPy.org/ doc/scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>} or U{interp2d<https://docs.SciPy.org/doc/scipy/reference/generated/ scipy.interpolate.interp2d.html>} interpolator and converted from 2-byte C{int} to 8-byte C{dtype float64}. Therefore, internal memory usage is 4x the U{egm*.pgm<https://GeographicLib.SourceForge.io/html/ geoid.html#geoidinst>} file size and may exceed the available memory, especially with 32-bit Python, see properties C{.nBytes} and C{.sizeB}. '''
kind=3, name=NN, smooth=0): '''New L{GeoidPGM} interpolator.
@arg egm_pgm: An U{EGM geoid dataset<https://GeographicLib.SourceForge.io/ html/geoid.html#geoidinst>} file name (C{egm*.pgm}). @kwarg crop: Optional box to crop B{C{egm_pgm}}, a 4-tuple (C{south, west, north, east}) or 2-tuple (C{(south, west), (north, east)}), in C{degrees90} lat- and C{degrees180} longitudes or a 2-tuple (C{LatLonSW, LatLonNE}) of C{LatLon} instances. @kwarg datum: Optional grid datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}), default C{WGS84}. @kwarg kind: C{scipy.interpolate} order (C{int}), use 1..5 for U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/ reference/generated/scipy.interpolate.RectBivariateSpline.html>}, -2 for U{interp2d linear<https://docs.SciPy.org/doc/scipy/ reference/generated/scipy.interpolate.interp2d.html>}, -3 for C{interp2d cubic} or -5 for C{interp2d quintic}. @kwarg name: Optional geoid name (C{str}). @kwarg smooth: Smoothing factor for U{RectBivariateSpline <https://docs.SciPy.org/doc/scipy/reference/generated/ scipy.interpolate.RectBivariateSpline.html>} only (C{int}).
@raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid B{C{crop}}, B{C{kind}} or B{C{smooth}}.
@raise ImportError: Package C{numpy} or C{scipy} not found or not installed.
@raise LenError: EGM dataset B{C{egm_pgm}} axis mismatch.
@raise SciPyError: A C{RectBivariateSpline} or C{inter2d} issue.
@raise SciPyWarning: A C{RectBivariateSpline} or C{inter2d} warning as exception.
@raise TypeError: Invalid B{C{datum}}.
@note: The U{GeographicLib egm*.pgm<https://GeographicLib.SourceForge.io/ html/geoid.html#geoidinst>} file sizes are based on a 2-byte C{int} height converted to 8-byte C{dtype float64} for C{scipy} interpolators. Therefore, internal memory usage is 4 times the C{egm*.pgm} file size and may exceed the available memory, especially with 32-bit Python. To reduce memory usage, set keyword argument B{C{crop}} to the region of interest. For example C{B{crop}=(20, -125, 50, -65)} covers the U{conterminous US<https:// www.NGS.NOAA.gov/GEOID/GEOID12B/maps/GEOID12B_CONUS_grids.png>} (CONUS), less than 3% of the entire C{egm2008-1.pgm} dataset.
@see: Class L{GeoidKarney} and function L{egmGeoidHeights}. '''
g = open(g.name, 'rb') # reopen tempfile for numpy 1.8.0- else: self._cropped = False
# U{numpy dtype formats are different from Python struct formats # <https://docs.SciPy.org/doc/numpy-1.15.0/reference/arrays.dtypes.html>} # read all heights, skipping the PGM header lines, converted to float except Exception as x: raise _SciPyIssue(x, _in_, repr(egm_pgm)) finally:
# convert grid (lat, lon) to earth (lat, lon), uncropped while lon > _180_0: lon -= _360_0 return lat, lon
# convert cropped grid (lat, lon) to earth (lat, lon)
# convert earth (lat, lon) to grid (lat, lon), uncropped while lon < 0: lon += _360_0 return lat, lon
# convert earth (lat, lon) to cropped grid (lat, lon)
if _FOR_DOCS: __call__ = _GeoidBase.__call__ height = _GeoidBase.height
'''Get the PGM itemsize in bytes (C{int}). '''
'''(INTERNAL) Basic geoid parameters. ''' # interpolator parameters
t = _COMMASPACE_.join(pairs((a, getattr(self, a)) for a in dir(self.__class__) if a[:1].isupper())) return _COLONSPACE_(self, t)
return Fmt.PAREN(self.classname, repr(self.name))
'''(INTERNAL) Parse an C{egm*.pgm} geoid dataset file.
# Geoid file in PGM format for the GeographicLib::Geoid class # Description WGS84 EGM96, 5-minute grid # URL https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/egm96/egm96.html # DateTime 2009-08-29 18:45:03 # MaxBilinearError 0.140 # RMSBilinearError 0.005 # MaxCubicError 0.003 # RMSCubicError 0.001 # Offset -108 # Scale 0.003 # Origin 90N 0E # AREA_OR_POINT Point # Vertical_Datum WGS84 <width> <height> <pixel> ... ''' # pgm = NN # name
# llstr to (lat, lon) floats
# PGM file attributes, CamelCase but not .istitle()
'''(INTERNAL) New C{_PGM} parsed C{egm*.pgm} geoid dataset. '''
raise self._Errorf(_format_, _header_, t)
while True: # read all # Attr ... lines, except (TypeError, ValueError): raise self._Errorf(_format_, 'Attr', t) else: # should never get here raise self._Errorf(_format_, _header_, g.tell())
nlat < 2 or nlat > (181 * 60) or not isodd(nlat): raise ValueError except (TypeError, ValueError): raise self._Errorf(_format_, _SPACE_(_width_, _height_), t)
raise ValueError except (TypeError, ValueError): raise self._Errorf(_format_, 'pixel', t)
setattr(self, a, None)
raise self._Errorf(_format_, 'Origin', self.Origin) raise self._Errorf(_format_, 'Offset', self.Offset) raise self._Errorf(_format_, 'Scale', self.Scale)
# note, negative .dlat and .rlat since rows # are from .slat 90N down in decreasing lat
# grid corners in earth (lat, lon), .slat = 90, .dlat < 0
raise self._Errorf('%s(%s x %s != %s)', _assert_, nlat, nlon, n)
'''Crop the geoid to (south, west, north, east) box. ''' # flon offset for both west and east # earth (lat, lon) to grid indices (y, x), # note y is decreasing, i.e. n < s
# handle special cases n, s = 0, hi # entire lat range w, e, f = 0, wi, 180 # entire lon range return g # use entire geoid as-is
raise self._Errorf(_format_, 'swne', (north - south, east - west))
# read w..wi and 0..e r, p = (wi - w), (e - wi) else: raise self._Errorf('%s(%s < %s)', _assert_, w, e)
# convert to bytes # number of rows and cols to skip from # the original (.slat, .wlon) origin # sanity check or z > self.sizeB: raise self._Errorf(_format_, _assert_, (r, p, q, z))
# can't use _BytesIO since numpy # needs .fileno attr in .fromfile # reading (s - n) rows, forward # Python 2 tmpfile.write returns None g.seek(-q, _SEEK_CUR) # assert(g.tell() == (z - w * self.u2B)) # Python 2 tmpfile.write returns None t += c.write(g.read(p)) or p
raise self._Errorf('%s(%s != %s) %s', _assert_, t, z, self)
# update the _Gpars accordingly, note attributes # .dlat, .dlon, .rlat and .rlon remain unchanged
# c = open(c.name, 'rb') # reopen for numpy 1.8.0-
def _Errorf(self, fmt, *args): # PYCHOK no cover t = fmt % args e = self.pgm or NN if e: t = _SPACE_(t, _in_, repr(e)) return PGMError(t)
# earth (lat, lon) to grid indices (y, x) # with .dlat decreasing from 90N .slat max(0, int(lon * self.rlon)))
# create a tmpfile to hold the cropped geoid grid except ImportError: # Python 2.7.16- from os import tmpfile
'''Get the geoid file name (C{str}). '''
'''Issue parsing or cropping an C{egm*.pgm} geoid dataset. '''
'''Generate geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/ html/geoid.html#geoidinst>} height tests from U{GeoidHeights.dat <https://SourceForge.net/projects/geographiclib/files/testdata/>} U{Test data for Geoids<https://GeographicLib.SourceForge.io/html/ geoid.html#testgeoid>}.
@arg GeoidHeights_dat: The un-gz-ed C{GeoidHeights.dat} file (C{str} or C{file} handle).
@return: For each test, yield a L{GeoidHeight5Tuple}C{(lat, lon, egm84, egm96, egm2008)}.
@raise GeoidError: Invalid B{C{GeoidHeights_dat}}.
@note: Function L{egmGeoidHeights} is used to test the geoids L{GeoidKarney} and L{GeoidPGM}, see PyGeodesy module C{test/testGeoids.py}. '''
except AttributeError as x: raise GeoidError(GeoidHeights_dat=type(dat), txt=str(x))
if __name__ == '__main__':
import sys
_crop = () _GeoidEGM = GeoidKarney _kind = 3
geoids = sys.argv[1:] while geoids: geoid = geoids.pop(0)
if '-crop'.startswith(geoid.lower()): _crop = 20, -125, 50, -65 # CONUS
elif '-karney'.startswith(geoid.lower()): _GeoidEGM = GeoidKarney
elif '-kind'.startswith(geoid.lower()): _kind = int(geoids.pop(0))
elif '-pgm'.startswith(geoid.lower()): _GeoidEGM = GeoidPGM
elif geoid[-4:].lower() in ('.pgm',): g = _GeoidEGM(geoid, crop=_crop, kind=_kind) print('\n%s\n' % (g.toStr(),)) print('%r\n' % (g.pgm,)) # <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>: # The height of the EGM96 geoid at Timbuktu # echo 16:46:33N 3:00:34W | GeoidEval # => 28.7068 -0.02e-6 -1.73e-6 # The 1st number is the height of the geoid, the 2nd and # 3rd are its slopes in northerly and easterly direction t = 'Timbuktu %s' % (g,) k = {'egm84-15.pgm': '31.2979', 'egm96-5.pgm': '28.7067', 'egm2008-1.pgm': '28.7880'}.get(g.name.lower(), '28.7880') ll = parseDMS2('16:46:33N', '3:00:34W', sep=':') for ll in (ll, (16.776, -3.009),): try: h, ll = g.height(*ll), fstr(ll, prec=6) print('%s.height(%s): %.4F vs %s' % (t, ll, h, k)) except (GeoidError, RangeError) as x: print(_COLONSPACE_(t, str(x)))
elif geoid[-4:].lower() in ('.bin',): g = GeoidG2012B(geoid, kind=_kind) print(g.toStr())
else: raise GeoidError(grid=repr(geoid))
# **) MIT License # # Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE.
# <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval> # _lowerleft = -90, -179, -30.1500 # egm2008-1.pgm # _lowerleft = -90, -179, -29.5350 # egm96-5.pgm # _lowerleft = -90, -179, -29.7120 # egm84-15.pgm
# _center = 0, 0, 17.2260 # egm2008-1.pgm # _center = 0, 0, 17.1630 # egm96-5.pgm # _center = 0, 0, 18.3296 # egm84-15.pgm
# _upperright = 90, 180, 14.8980 # egm2008-1.pgm # _upperright = 90, 180, 13.6050 # egm96-5.pgm # _upperright = 90, 180, 13.0980 # egm84-15.pgm
# % python pygeodesy/geoids.py [-Karney] ../geoids/egm*.pgm # # 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, -32.633, 85.839), lowest(4.683, -101.25, -106.911), mean=-1.317, stdev=29.244, # kind=3, dtype='ushort', endian='>H', hits=0, knots=233301600, sizeB=466603604 # _PGM('../geoids/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=(90, 0.0), Pixel=65535, RMSBilinearError=0.001, RMSCubicError=0.001, Scale=0.003, # URL='https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/egm2008', Vertical_Datum='WGS84', crop4=(-90.0, -180.0, 90.0, 180.0), # dlat=-0.016666666666666666, dlon=0.016666666666666666, egm=None, flon=0, glon=180, knots=233301600, nlat=10801, nlon=21600, # pgm='../geoids/egm2008-1.pgm', rlat=-60.0, rlon=60.0, sizeB=466603604, skip=404, slat=90, u2B=2, wlon=0.0 # Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880 # Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880 # # 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(-8.4, -32.633, 85.839), lowest(4.683, -101.25, -106.911), mean=-1.317, stdev=29.244, # kind=3, dtype='ushort', endian='>H', hits=0, knots=1038240, sizeB=2076896 # _PGM('../geoids/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=(90, 0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, # URL='https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84', crop4=(-90.0, -180.0, 90.0, 180.0), # dlat=-0.25, dlon=0.25, egm=None, flon=0, glon=180, knots=1038240, nlat=721, nlon=1440, # pgm='../geoids/egm84-15.pgm', rlat=-4.0, rlon=4.0, sizeB=2076896, skip=416, slat=90, u2B=2, wlon=0.0 # Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979 # Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979 # # 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.4, -32.633, 85.839), lowest(4.683, -101.25, -106.911), mean=-1.317, stdev=29.244, # kind=3, dtype='ushort', endian='>H', hits=0, knots=9335520, sizeB=18671448 # _PGM('../geoids/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=(90, 0.0), Pixel=65535, RMSBilinearError=0.005, RMSCubicError=0.001, Scale=0.003, # URL='https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/egm96/egm96.html', Vertical_Datum='WGS84', crop4=(), # dlat=-0.08333333333333333, dlon=0.08333333333333333, egm=None, flon=0, glon=180, knots=9335520, nlat=2161, nlon=4320, # pgm='../geoids/egm96-5.pgm', rlat=-12.0, rlon=12.0, sizeB=18671448, skip=408, slat=90, u2B=2, wlon=0.0 # Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067 # Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067
# % python pygeodesy/geoids.py -PGM ../geoids/egm*.pgm # # 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), mean=-1.317, stdev=29.244, # kind=3, smooth=0, dtype=dtype('float64'), endian='>u2', knots=233301600, nBytes=1866412800, sizeB=466603604, scipy='1.2.1', numpy='1.16.1' # _PGM('../geoids/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=(90, 0.0), Pixel=65535, RMSBilinearError=0.001, RMSCubicError=0.001, Scale=0.003, # URL='https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/egm2008', Vertical_Datum='WGS84', crop4=(-90.0, -180.0, 90.0, 180.0), # dlat=-0.016666666666666666, dlon=0.016666666666666666, egm=None, flon=0, glon=180, knots=233301600, nlat=10801, nlon=21600, # pgm='../geoids/egm2008-1.pgm', rlat=-60.0, rlon=60.0, sizeB=466603604, skip=404, slat=90, u2B=2, wlon=0.0 # Timbuktu GeoidPGM('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880 # Timbuktu GeoidPGM('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880 # # 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), mean=-0.855, stdev=29.183, # kind=3, smooth=0, dtype=dtype('float64'), endian='>u2', knots=1038240, nBytes=8305920, sizeB=2076896, scipy='1.2.1', numpy='1.16.1' # _PGM('../geoids/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=(90, 0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, # URL='https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84', crop4=(-90.0, -180.0, 90.0, 180.0), # dlat=-0.25, dlon=0.25, egm=None, flon=0, glon=180, knots=1038240, nlat=721, nlon=1440, # pgm='../geoids/egm84-15.pgm', rlat=-4.0, rlon=4.0, sizeB=2076896, skip=416, slat=90, u2B=2, wlon=0.0 # Timbuktu GeoidPGM('egm84-15.pgm').height(16.775833, -3.009444): 31.2979 vs 31.2979 # Timbuktu GeoidPGM('egm84-15.pgm').height(16.776, -3.009): 31.2975 vs 31.2979 # # 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), mean=-1.438, stdev=29.227, # kind=3, smooth=0, dtype=dtype('float64'), endian='>u2', knots=9335520, nBytes=74684160, sizeB=18671448, scipy='1.2.1', numpy='1.16.1' # _PGM('../geoids/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=(90, 0.0), Pixel=65535, RMSBilinearError=0.005, RMSCubicError=0.001, Scale=0.003, # URL='https://Earth-Info.NGA.mil/GandG/wgs84/gravitymod/egm96/egm96.html', Vertical_Datum='WGS84', crop4=(-90.0, -180.0, 90.0, 180.0), # dlat=-0.08333333333333333, dlon=0.08333333333333333, egm=None, flon=0, glon=180, knots=9335520, nlat=2161, nlon=4320, # pgm='../geoids/egm96-5.pgm', rlat=-12.0, rlon=12.0, sizeB=18671448, skip=408, slat=90, u2B=2, wlon=0.0 # Timbuktu GeoidPGM('egm96-5.pgm').height(16.775833, -3.009444): 28.7065 vs 28.7067 # Timbuktu GeoidPGM('egm96-5.pgm').height(16.776, -3.009): 28.7064 vs 28.7067 |