1

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

3

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

5

6 Classes L{GeoidG2012B}, L{GeoidKarney} and L{GeoidPGM} to interpolate the

7 height of various U{geoid<https://WikiPedia.org/wiki/Geoid>}s at C{LatLon}

8 locations or separate lat-/longitudes using different interpolation methods

9 and C{geoid} model files.

10

11 L{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.

13 The L{GeoidG2012B} and L{GeoidPGM} interpolators both depend on U{scipy

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

15 require those packages to be installed.

16

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

18 a C{geoid} model file, I{specific to the interpolator}, more details below

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

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

21

22 Typical usage

23=============

24

25 1. Choose one of the interpolator classes L{GeoidG2012B}, L{GeoidKarney}

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

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

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

29

30 C{>>> from pygeodesy import GeoidG2012B # or -Karney or -PGM as GeoidXyz}

31

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

33 arguments to select different interpolation options

34

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

36

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

38

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

40

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

42

43 or

44

45 C{>>> h0, h1, h2, ... = ginterpolator(ll0, ll1, ll2, ...)}

46

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

48

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

50

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

52

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

54

55 or as 2 lists, 2 tuples, etc.

56

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

58

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

60 courtesy of SBFRF.

61

62@note: Classes L{GeoidG2012B} and L{GeoidPGM} require both U{numpy

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

64 to be installed.

65

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

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

68 accordingly, see L{SciPyWarning}.

69

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

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

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

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

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

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

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

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

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

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

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

81'''

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

83 from __future__ import division as _; del _ # PYCHOK semicolon

84

85 from pygeodesy.basics import len2, map1, isodd, ub2str as _ub2str

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

87 _N_180_0, _360_0

88# from pygeodesy.datums import _ellipsoidal_datum # from .heights

89# from pygeodesy.dms import parseDMS2 # _MODS

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

91 _xkwds_pop2

92 from pygeodesy.fmath import favg, Fdot, fdot, Fhorner, frange

93# from pygoedesy.formy import heightOrthometric # _MODS

94 from pygeodesy.heights import _as_llis2, _ascalar, _HeightBase, HeightError, \

95 _ellipsoidal_datum, _Wrap

96# from pygeodesy.internals import _version2 # _MODS

97 from pygeodesy.interns import NN, _COLONSPACE_, _COMMASPACE_, _E_, _height_, \

98 _in_, _kind_, _lat_, _lon_, _mean_, _N_, _n_a_, \

99 _numpy_, _on_, _outside_, _S_, _s_, _scipy_, \

100 _SPACE_, _stdev_, _tbd_, _W_, _width_, _4_

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

102 from pygeodesy.named import _name__, _Named, _NamedTuple

103# from pygeodesy.namedTuples import LatLon3Tuple # _MODS

104 from pygeodesy.props import Property_RO, property_RO, property_ROver

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

106 from pygeodesy.units import Height, Int_, Lat, Lon

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

108

109 from math import floor as _floor

110# from os import SEEK_CUR, SEEK_SET # _MODS

111# import os.path # _MODS

112 from struct import calcsize as _calcsize, unpack as _unpack

113 try:

114 from StringIO import StringIO as _BytesIO # reads bytes

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

116

117 except ImportError: # Python 3+

118 from io import BytesIO as _BytesIO # PYCHOK expected

119

120__all__ = _ALL_LAZY.geoids

121__version__ = '24.10.11'

122

123_assert_ = 'assert'

124_bHASH_ = b'#'

125_endian_ = 'endian'

126_format_ = '%s %r'

127_header_ = 'header'

128_intCs = {} # cache int value

129_lli_ = 'lli'

130_non_increasing_ = 'non-increasing'

131_rb_ = 'rb'

132_supported_ = 'supported'

133

134

135 class GeoidError(HeightError):

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

137 '''

138 pass

139

140

141 class _GeoidBase(_HeightBase):

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

143 '''

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

145 _cropped = None

146# _datum = _WGS84 # from _HeightBase

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

148 _endian = _tbd_

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

150 _geoid = _n_a_

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

152 _kind = 3 # order for interp2d, RectBivariateSpline

153# _kmin = 2 # min number of knots

154 _knots = 0 # nlat * nlon

155 _mean = None # fixed in GeoidKarney

156 _nBytes = 0 # numpy size in bytes, float64

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

158 _sizeB = 0 # geoid file size in bytes

159 _smooth = 0 # used only for RectBivariateSpline

160 _stdev = None # fixed in GeoidKarney

161 _u2B = 0 # np.itemsize or undefined

162 _yx_hits = None # cache hits, ala Karney

163

164# _lat_d = _0_0 # increment, +tive

165# _lat_lo = _0_0 # lower lat, south

166# _lat_hi = _0_0 # upper lat, noth

167# _lon_d = _0_0 # increment, +tive

168# _lon_lo = _0_0 # left lon, west

169# _lon_hi = _0_0 # right lon, east

170# _lon_of = _0_0 # forward lon offset

171# _lon_og = _0_0 # reverse lon offset

172

173 def __init__(self, hs, p):

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

175 several internal geoid attributes.

176

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

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

179 geoid parameters (C{INTERNAL}).

180 '''

181 spi = self.scipy_interpolate

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

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

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

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

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

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

188

189 # both axes and bounding box

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

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

192

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

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

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

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

197 k = self.kind

198 if k in self._k2interp2d: # see _HeightBase

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

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

201 k = self._kxky(k)

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

203 s=self._smooth).ev

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

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

206 self._knots = p.knots # grid knots

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

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

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

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

211 self._lat_lo, \

212 self._lat_hi, \

213 self._lon_lo, \

214 self._lon_hi = map(float, bb)

215 self._lon_lo -= g

216 self._lon_hi -= g

217

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

219 '''Interpolate the geoid height for one or several locations.

220

221 @arg llis: One or more locations (C{LatLon}s), all positional.

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

223 C{B{H}=False} (C{bool}). If C{B{wrap} is True},

224 wrap or I{normalize} all B{C{llis}} locations. If

225 C{B{H} is True}, return the I{orthometric} height

226 instead of the I{geoid} height at each location.

227

228 @return: A single interpolated geoid (or orthometric) height

229 (C{float}) or a list or tuple of interpolated geoid

230 (or orthometric) heights (C{float}s).

231

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

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

234

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

236 longitude range.

237

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

239

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

241

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

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

244 C{height=0} is used.

245

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

247 '''

248 return self._called(llis, True, **wrap_H)

249

250 def __enter__(self):

251 '''Open context.

252 '''

253 return self

254

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

256 '''Close context.

257 '''

258 self.close()

259 # return None # XXX False

260

261 def __repr__(self):

262 return self.toStr()

263

264 def __str__(self):

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

266

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

268 # handle __call__

269 _H = self._heightOrthometric if H else None

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

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

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

273 try:

274 for i, lli in enumerate(llis):

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

276 # orthometric or geoid height

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

278 return _as(hs)

279 except (GeoidError, RangeError) as x:

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

281 t = _lli_ if _as is _ascalar else Fmt.INDEX(llis=i)

282 lli = fstr((lli.lat, lli.lon), strepr=repr)

283 raise type(x)(t, lli, wrap=wrap, H=H, cause=x)

284 except Exception as x:

285 if iscipy and self.scipy:

286 raise _SciPyIssue(x, self._ev_name)

287 else:

288 raise

289

290 @Property_RO

291 def _center(self):

292 ''' Cache for method L{center}.

293 '''

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

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

296

297 def center(self, LatLon=None):

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

299

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

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

302

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

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

305 longitude and geoid height of the center grid location.

306 '''

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

308

309 def close(self):

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

311 '''

312 if not self.closed:

313 self._egm.close()

314 self._egm = None

315

316 @property_RO

317 def closed(self):

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

319 '''

320 return self._egm is None

321

322 @Property_RO

323 def cropped(self):

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

325 '''

326 return self._cropped

327

328 @Property_RO

329 def dtype(self):

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

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

332 '''

333 return self._hs_y_x.dtype

334

335 @Property_RO

336 def endian(self):

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

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

339 '''

340 return self._endian

341

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

343 # see methods _HeightBase._ev and -._interp2d

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

345

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

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

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

349 if m != n:

350 raise LenError(self.__class__, grid=m, **{name: n})

351 if d < 0:

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

353 for i in range(1, m):

354 e = a[i] - a[i-1]

355 if e < EPS: # non-increasing axis

356 i = Fmt.INDEX(name, i)

357 raise GeoidError(i, e, txt=_non_increasing_)

358 return self.numpy.array(a), d

359

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

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

362 self._notOverloaded(lat, lon)

363

364 def _gyx2g2(self, y, x):

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

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

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

368

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

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

371

372 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).

373 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).

374 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}

375 and B{C{lons}} locations (C{bool}).

376

377 @return: A single interpolated geoid height (C{float}) or a

378 list of interpolated geoid heights (C{float}s).

379

380 @raise GeoidError: Insufficient or non-matching number of

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

382

383 @raise RangeError: A B{C{lat}} or B{C{lon}} is outside this

384 geoid's lat- or longitude range.

385

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

387

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

389 '''

390 lls = self._as_lls(lats, lons)

391 return self(lls, **wrap) # __call__(ll) or __call__(lls)

392

393 @property_ROver

394 def _heightOrthometric(self):

395 return _MODS.formy.heightOrthometric # overwrite property_ROver

396

397 def _hGeoid(self, lat, lon):

398 out = self.outside(lat, lon)

399 if out:

400 lli = fstr((lat, lon), strepr=repr)

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

402 return float(self._ev(*self._ll2g2(lat, lon)))

403

404 @Property_RO

405 def _highest(self):

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

407 '''

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

409

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

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

412

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

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

415

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

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

418 longitude and geoid height of the highest grid location.

419 '''

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

421

422 @Property_RO

423 def hits(self):

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

425 '''

426 return self._yx_hits

427

428 @Property_RO

429 def kind(self):

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

431 '''

432 return self._kind

433

434 @Property_RO

435 def knots(self):

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

437 '''

438 return self._knots

439

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

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

442 self._notOverloaded(lat, lon)

443

444 @property_ROver

445 def _LL3T(self):

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

447 '''

448 return _MODS.namedTuples.LatLon3Tuple # overwrite property_ROver

449

450 def _llh3(self, lat, lon):

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

452

453 def _llh3LL(self, llh, LatLon):

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

455

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

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

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

459 # numpy.argmin.html#numpy.argmin>

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

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

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

463

464 def _load(self, g, dtype, n, offset=0):

465 # numpy.fromfile, like .frombuffer

466 g.seek(offset, _MODS.os.SEEK_SET)

467 return self.numpy.fromfile(g, dtype, n)

468

469 @Property_RO

470 def _lowerleft(self):

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

472 '''

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

474

475 def lowerleft(self, LatLon=None):

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

477

478 @kwarg LatLon: Optional class to return the location

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

480

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

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

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

484 '''

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

486

487 @Property_RO

488 def _loweright(self):

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

490 '''

491 return self._llh3(self._lat_lo, self._lon_hi)

492

493 def loweright(self, LatLon=None):

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

495

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

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

498

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

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

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

502 '''

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

504

505 lowerright = loweright # synonymous

506

507 @Property_RO

508 def _lowest(self):

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

510 '''

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

512

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

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

515

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

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

518

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

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

521 longitude and geoid height of the lowest grid location.

522 '''

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

524

525 @Property_RO

526 def mean(self):

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

528 '''

529 if self._mean is None: # see GeoidKarney

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

531 return self._mean

532

533 @property_RO

534 def name(self):

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

536 '''

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

538

539 @Property_RO

540 def nBytes(self):

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

542 '''

543 return self._nBytes

544

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

546 # open the geoid file

547 try:

548 self._geoid = _MODS.os.path.basename(geoid)

549 self._sizeB = _MODS.os.path.getsize(geoid)

550 g = open(geoid, _rb_)

551 except (IOError, OSError) as x:

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

553

554 if datum not in (None, self._datum):

555 self._datum = _ellipsoidal_datum(datum, name=name)

556 self._kind = int(kind)

557 if name:

558 _HeightBase.name.fset(self, name) # rename

559 if smooth:

560 self._smooth = Int_(smooth=smooth, Error=GeoidError, low=0)

561

562 return g

563

564 def outside(self, lat, lon):

565 '''Check whether a location is outside this geoid's

566 lat-/longitude or crop range.

567

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

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

570

571 @return: A 1- or 2-character C{str} if outside or an

572 empty C{str} if inside.

573 '''

574 return (_S_ if lat < self._lat_lo else

575 (_N_ if lat > self._lat_hi else NN)) + \

576 (_W_ if lon < self._lon_lo else

577 (_E_ if lon > self._lon_hi else NN))

578

579 @Property_RO

580 def pgm(self):

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

582 '''

583 return self._pgm

584

585 @Property_RO

586 def sizeB(self):

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

588 '''

589 return self._sizeB

590

591 @Property_RO

592 def smooth(self):

593 '''Get the C{RectBivariateSpline} smoothing (C{int}).

594 '''

595 return self._smooth

596

597 @Property_RO

598 def stdev(self):

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

600 '''

601 if self._stdev is None: # see GeoidKarney

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

603 return self._stdev

604

605 def _swne(self, crop):

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

607 try:

608 if len(crop) == 2:

609 try: # sw, ne LatLons

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

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

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

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

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

615 swne = crop

616 if len(swne) == 4:

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

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

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

620 return s, w, n, e

621 except (IndexError, TypeError, ValueError):

622 pass

623 raise GeoidError(crop=crop)

624

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

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

627

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

629 default). Trailing zero decimals are stripped

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

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

632 @kwarg sep: Separator to join (C{str}).

633

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

635 '''

636 s = 1 if self.kind < 0 else 2

637 t = tuple(Fmt.PAREN(m.__name__, fstr(m(), prec=prec)) for m in

638 (self.lowerleft, self.upperright,

639 self.center,

640 self.highest, self.lowest)) + \

641 attrs( _mean_, _stdev_, prec=prec, Nones=False) + \

642 attrs((_kind_, 'smooth')[:s], prec=prec, Nones=False) + \

643 attrs( 'cropped', 'dtype', _endian_, 'hits', 'knots', 'nBytes',

644 'sizeB', _scipy_, _numpy_, prec=prec, Nones=False)

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

646

647 @Property_RO

648 def u2B(self):

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

650 '''

651 return self._u2B

652

653 @Property_RO

654 def _upperleft(self):

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

656 '''

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

658

659 def upperleft(self, LatLon=None):

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

661

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

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

664

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

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

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

668 '''

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

670

671 @Property_RO

672 def _upperright(self):

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

674 '''

675 return self._llh3(self._lat_hi, self._lon_hi)

676

677 def upperright(self, LatLon=None):

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

679

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

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

682

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

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

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

686 '''

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

688

689

690 class GeoidG2012B(_GeoidBase):

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

692 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/>} grids U{CONUS

693 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_CONUS.shtml>},

694 U{Alaska<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_AK.shtml>},

695 U{Hawaii<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_HI.shtml>},

696 U{Guam and Northern Mariana Islands

697 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_GMNI.shtml>},

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

699 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_PRVI.shtml>} and

700 U{American Samoa<https://www.NGS.NOAA.gov/GEOID/GEOID12B/GEOID12B_AS.shtml>}

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

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

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

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

705 generated/scipy.interpolate.bisplrep.html>} interpolation.

706

707 Use any of the binary C{le} (little endian) or C{be} (big endian)

708 C{g2012b*.bin} grid files.

709 '''

710 def __init__(self, g2012b_bin, datum=None, # NAD 83 Ellipsoid

711 kind=3, smooth=0, **name_crop):

712 '''New L{GeoidG2012B} interpolator.

713

714 @arg g2012b_bin: A C{GEOID12B} grid file name (C{.bin}).

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

716 L{a_f2Tuple}), default C{WGS84}.

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

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

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

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

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

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

723 see note for more details.

724 @kwarg smooth: Smoothing factor for C{B{kind}=1..5} only (C{int}).

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

726 C{B{crop}=None}.

727

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

729 B{C{g2012b_bin}} issue.

730

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

732

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

734

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

736

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

738

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

740

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

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

743 '''

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

745 if crop is not None:

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

747

748 g = self._open(g2012b_bin, datum, kind, _name__(**name), smooth)

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

750

751 try:

752 p = _Gpars()

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

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

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

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

757 # skip 4xf8, get 3xi4

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

759 if ien == 1: # correct endian

760 p.knots = p.nlat * p.nlon

761 if p.knots == n and 1 < p.nlat < n \

762 and 1 < p.nlon < n:

763 self._endian = en_+'f4'

764 break

765 else: # couldn't validate endian

766 raise GeoidError(_endian_)

767

768 # get the first 4xf8

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

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

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

772 p.wlon -= _360_0 # western-most East longitude to earth (..., lon)

773 _GeoidBase.__init__(self, hs, p)

774

775 except Exception as x:

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

777 finally:

778 g.close()

779

780 def _g2ll2(self, lat, lon):

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

782 return lat, lon

783

784 def _ll2g2(self, lat, lon):

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

786 return lat, lon

787

788 if _FOR_DOCS:

789 __call__ = _GeoidBase.__call__

790 height = _GeoidBase.height

791

792

793 class GeoidHeight5Tuple(_NamedTuple): # .geoids.py

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

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

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

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

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

799 '''

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

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

802

803

804 def _I(i):

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

806 '''

807 return _intCs.setdefault(i, i) # PYCHOK undefined due to del _intCs

808

809

810 def _T(*cs):

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

812 '''

813 return map1(_I, *cs)

814

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

816

817

818 class GeoidKarney(_GeoidBase):

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

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

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

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

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

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

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

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

827

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

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

830 datasets.

831 '''

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

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

834 _C3 = ((_T(0, 0, 62, 124, 124, 62, 0, 0, 0, 0, 0, 0),

835 _T0s12,

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

837 _T0s12,

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

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

840 _T0s12,

841 _T(0, 0, 0, 0, 0, 0, 93, -93, -93, 93, 0, 0),

842 _T(-102, 102, 0, 0, 0, 0, 18, 12, -12, -18, -84, 84),

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

844

845 (_T(9, -9, 9, 186, 54, -9, -9, 54, -54, 9, -9, 9),

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

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

848 _T(0, 0, 90, -150, 30, 30, 30, -90, 90, -30, 0, 0),

849 _T(96, -96, 96, -96, -24, 24, -96, -24, 144, -24, 24, -24), # PYCHOK indent

850 _T(90, 30, 0, -150, -90, 0, 0, 30, 90, 0, 30, -30),

851 _T(0, 0, -20, 60, -60, 20, -20, 60, -60, 20, 0, 0),

852 _T(0, 0, -60, 60, 60, -60, 60, -60, -60, 60, 0, 0),

853 _T(-60, 60, 0, 60, -60, 0, 0, 60, -60, 0, -60, 60),

854 _T(-20, -20, 0, 60, 60, 0, 0, -60, -60, 0, 20, 20)),

855

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

857 _T(-36, 36, -165, 45, 141, -21, 0, 0, 0, 0, 36, -36),

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

859 _T(0, 0, 93, -93, -93, 93, 0, 0, 0, 0, 0, 0),

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

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

862 _T0s12,

863 _T(0, 0, -93, 93, 93, -93, 0, 0, 0, 0, 0, 0),

864 _T(-84, 84, 18, 12, -12, -18, 0, 0, 0, 0, -102, 102),

865 _T(-31, -31, 0, 93, 93, 0, -31, -62, -62, -31, 31, 31)))

866

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

868 _T(1, 0),

869 _T(0, 1),

870 _T(1, 1))

871

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

873 _T( 1, -1),

874 _T(-1, 0),

875 _T( 0, 0),

876 _T( 1, 0),

877 _T( 2, 0),

878 _T(-1, 1),

879 _T( 0, 1),

880 _T( 1, 1),

881 _T( 2, 1),

882 _T( 0, 2),

883 _T( 1, 2))

884

885# _cropped = None

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

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

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

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

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

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

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

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

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

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

896 _nBytes = None # not applicable

897 _nterms = len(_C3[0]) # columns length, number of row

898 _smooth = None # not applicable

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

900 _u2B = _calcsize(_endian) # pixelsize_ in bytes

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

902 _Ru2B = 0 # row of pixelsize_s in bytes

903 _yx_hits = 0 # cache hits

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

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

906

907 def __init__(self, egm_pgm, crop=None, datum=None, # WGS84

908 kind=3, **name_smooth):

909 '''New L{GeoidKarney} interpolator.

910

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

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

913 note below.

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

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

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

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

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

919 L{a_f2Tuple}), default C{WGS84}.

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

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

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

923

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

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

926

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

928

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

930

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

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

933 '''

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

935 if smooth is not None:

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

937

938 if kind in (2,):

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

940 elif kind not in (3,):

941 raise GeoidError(kind=kind)

942

943 self._egm = g = self._open(egm_pgm, datum, kind, _name__(**name), None)

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

945

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

947 self._Ru2B = _calcsize(self._Rendian)

948

949 self._knots = p.knots # grid knots

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

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

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

953 self._lat_lo, self._lon_lo, \

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

955 self._cropped = bool(crop)

956

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

958 '''Interpolate the geoid height for one or several locations.

959

960 @arg llis: One or more locations (C{LatLon}s), all positional.

961 @kwarg wrap_H: Keyword arguments C{B{wrap}=False, B{H}=False}.

962 If C{B{wrap} is True}, wrap or I{normalize} all

963 B{C{llis}} locations (C{bool}). If C{B{H} is True},

964 return the I{orthometric} height instead of the

965 I{geoid} height at each location (C{bool}).

966

967 @return: A single interpolated geoid (or orthometric) height

968 (C{float}) or a list or tuple of interpolated geoid

969 (or orthometric) heights (C{float}s).

970

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

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

973

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

975 longitude range.

976

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

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

979 C{height=0} is used.

980

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

982 '''

983 return self._called(llis, False, **wrap_H)

984

985 def _c0c3v(self, y, x):

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

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

988 p = self._pgm

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

990 # read 4x4 ushorts, drop the 4 corners

991 S = _MODS.os.SEEK_SET

992 e = self._4endian

993 g = self._egm

994 n = self._4u2B

995 R = self._Ru2B

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

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

998 b += R

999 g.seek(b, S)

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

1001 b += R

1002 g.seek(b, S)

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

1004 b += R

1005 g.seek(b, S)

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

1007 j = 1

1008

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

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

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

1012

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

1014

1015 @Property_RO

1016 def dtype(self):

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

1018 '''

1019 return 'ushort'

1020

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

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

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

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

1025 fy -= y

1026 fx -= x

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

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

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

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

1031

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

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

1034 if self._yx_i == yx:

1035 self._yx_hits += 1

1036 else:

1037 y, x = self._yx_i = yx

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

1039 t = self._yx_t

1040 v = _1_0, -fx, fx

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

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

1043 return H

1044

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

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

1047 if self._yx_i == yx:

1048 self._yx_hits += 1

1049 else:

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

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

1052 self._yx_t = tuple(fdot(v, *r3) / c0 for r3 in c3)

1053 self._yx_i = yx

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

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

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

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

1058 t = self._yx_t

1059 v = _1_0, fx, fy

1060 H = Fdot(v, t[5], t[8], t[9])

1061 H *= fy

1062 H += Fhorner(fx, t[2], t[4], t[7])

1063 H *= fy

1064 H += Fhorner(fx, t[0], t[1], t[3], t[6])

1065 return H

1066

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

1068

1069 def _g2ll2(self, lat, lon):

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

1071 while lon > _180_0:

1072 lon -= _360_0

1073 return lat, lon

1074

1075 def _g2yx2(self, lat, lon):

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

1077 p = self._pgm

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

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

1080

1081 def _gyx2g2(self, y, x):

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

1083 p = self._pgm

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

1085

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

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

1088

1089 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).

1090 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).

1091 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}

1092 and B{C{lons}} locations (C{bool}).

1093

1094 @return: A single interpolated geoid height (C{float}) or a

1095 list of interpolated geoid heights (C{float}s).

1096

1097 @raise GeoidError: Insufficient or non-matching number of

1098 B{C{lats}} and B{C{lons}} or the C{egm*.pgm}

1099 geoid file is closed.

1100

1101 @raise RangeError: A B{C{lat}} or B{C{lon}} is outside this

1102 geoid's lat- or longitude range.

1103 '''

1104 lls = self._as_lls(lats, lons)

1105 return self(lls, **wrap) # __call__(ll) or __call__(lls)

1106

1107 @Property_RO

1108 def _highest_ltd(self):

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

1110 '''

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

1112

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

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

1115

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

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

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

1119

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

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

1122 longitude and geoid height of the highest grid location.

1123 '''

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

1125 return self._llh3LL(llh, LatLon)

1126

1127 def _lat2y2(self, lat2):

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

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

1130 for lat in lat2:

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

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

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

1134

1135 def _ll2g2(self, lat, lon):

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

1137 while lon < 0:

1138 lon += _360_0

1139 return lat, lon

1140

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

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

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

1144 if highest:

1145 def _mt(r, h):

1146 m = max(r)

1147 return m, (m > h)

1148

1149 else: # lowest

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

1151 m = min(r)

1152 return m, (m < h)

1153

1154 y = x = 0

1155 h = self._raw(y, x)

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

1157 m, t = _mt(r, h)

1158 if t:

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

1160 h *= self._pgm.Scale

1161 h += self._pgm.Offset

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

1163

1164 @Property_RO

1165 def _lowest_ltd(self):

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

1167 '''

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

1169

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

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

1172

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

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

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

1176

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

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

1179 longitude and geoid height of the lowest grid location.

1180 '''

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

1182 return self._llh3LL(llh, LatLon)

1183

1184 def _raw(self, y, x):

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

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

1187 p = self._pgm

1188 if x < 0:

1189 x += p.nlon

1190 elif x >= p.nlon:

1191 x -= p.nlon

1192 h = p.nlon // 2

1193 if y < 0:

1194 y = -y

1195 elif y >= p.nlat:

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

1197 else:

1198 h = 0

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

1200 self._seek(y, x)

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

1202 return h[0]

1203

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

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

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

1207

1208 def _raw2(self, *lat2):

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

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

1211 p = self._pgm

1212 g = self._egm

1213 e = self._Rendian

1214 n = self._Ru2B

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

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

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

1218 for y in range(s, t):

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

1220

1221 def _seek(self, y, x):

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

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

1224 if g:

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

1226 g.seek(b, _MODS.os.SEEK_SET)

1227 return b # position

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

1229

1230

1231 class GeoidPGM(_GeoidBase):

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

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

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

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

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

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

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

1239 interp2d.html>} interpolation.

1240

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

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

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

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

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

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

1247 and C{.sizeB}.

1248 '''

1249 _cropped = False

1250 _endian = '>u2'

1251

1252 def __init__(self, egm_pgm, crop=None, datum=None, # WGS84

1253 kind=3, smooth=0, **name):

1254 '''New L{GeoidPGM} interpolator.

1255

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

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

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

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

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

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

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

1263 L{a_f2Tuple}), default C{WGS84}.

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

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

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

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

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

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

1270 see note for more details.

1271 @kwarg smooth: Smoothing factor for C{B{kind}=1..5} only (C{int}).

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

1273

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

1275 or B{C{smooth}}.

1276

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

1278

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

1280

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

1282

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

1284

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

1286

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

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

1289

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

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

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

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

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

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

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

1297 <https://www.NGS.NOAA.gov/GEOID/GEOID12B/maps/GEOID12B_CONUS_grids.png>}

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

1299

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

1301 '''

1302 np = self.numpy

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

1304

1305 g = self._open(egm_pgm, datum, kind, _name__(**name), smooth)

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

1307 if crop:

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

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

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

1311 self._cropped = True

1312 try:

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

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

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

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

1317 if p.Offset: # offset

1318 hs = p.Offset + hs

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

1320 hs = np.flipud(hs)

1321 _GeoidBase.__init__(self, hs, p)

1322 except Exception as x:

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

1324 finally:

1325 g.close()

1326

1327 def _g2ll2(self, lat, lon):

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

1329 if self._cropped:

1330 lon -= self._lon_of

1331 else:

1332 while lon > _180_0:

1333 lon -= _360_0

1334 return lat, lon

1335

1336 def _ll2g2(self, lat, lon):

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

1338 if self._cropped:

1339 lon += self._lon_of

1340 else:

1341 while lon < 0:

1342 lon += _360_0

1343 return lat, lon

1344

1345 if _FOR_DOCS:

1346 __call__ = _GeoidBase.__call__

1347 height = _GeoidBase.height

1348

1349

1350 class _Gpars(_Named):

1351 '''(INTERNAL) Basic geoid parameters.

1352 '''

1353 # interpolator parameters

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

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

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

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

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

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

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

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

1362

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

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

1365

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

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

1368

1369 def __repr__(self):

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

1371 a in dir(self.__class__)

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

1373 return _COLONSPACE_(self, t)

1374

1375 def __str__(self):

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

1377

1378

1379 class _PGM(_Gpars):

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

1381

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

1383 # Description WGS84 EGM96, 5-minute grid

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

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

1386 # MaxBilinearError 0.140

1387 # RMSBilinearError 0.005

1388 # MaxCubicError 0.003

1389 # RMSCubicError 0.001

1390 # Offset -108

1391 # Scale 0.003

1392 # Origin 90N 0E

1393 # AREA_OR_POINT Point

1394 # Vertical_Datum WGS84

1395 <width> <height>

1396 <pixel>

1397 ...

1398 '''

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

1400 egm = None

1401 glon = 180 # reverse offset, uncropped

1402# pgm = NN # name

1403 sizeB = 0

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

1405

1406 @staticmethod

1407 def _llstr2floats(latlon):

1408 # llstr to (lat, lon) floats

1409 lat, lon = latlon.split()

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

1411

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

1413 AREA_OR_POINT = str

1414 DateTime = str

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

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

1417 MaxBilinearError = float

1418 MaxCubicError = float

1419 Offset = float

1420 Origin = _llstr2floats

1421 Pixel = 0

1422 RMSBilinearError = float

1423 RMSCubicError = float

1424 Scale = float

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

1426 Vertical_Datum = str

1427

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

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

1430 '''

1431 self.name = pgm # geoid file name

1432 if itemsize:

1433 self._u2B = itemsize

1434 if sizeB:

1435 self.sizeB = sizeB

1436

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

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

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

1440

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

1442 try: # ignore empty ones or comments

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

1444 if t.startswith(_bHASH_):

1445 t = t.lstrip(_bHASH_).lstrip()

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

1447 f = getattr(_PGM, a, None)

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

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

1450 elif t:

1451 break

1452 except (TypeError, ValueError):

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

1454 else: # should never get here

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

1456

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

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

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

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

1461 raise ValueError

1462 except (TypeError, ValueError):

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

1464

1465 try: # must be 16 bit pixel height

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

1467 self.Pixel = int(t)

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

1469 raise ValueError

1470 except (TypeError, ValueError):

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

1472

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

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

1475 setattr(self, a, None)

1476

1477 if self.Origin is None:

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

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

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

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

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

1483

1484 self.skip = g.tell()

1485 self.knots = nlat * nlon

1486

1487 self.nlat, self.nlon = nlat, nlon

1488 self.slat, self.wlon = self.Origin

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

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

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

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

1493

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

1495 n = float(self.slat)

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

1497 w = self.wlon - self.glon

1498 e = w + self.dlon * nlon

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

1500

1501 n = self.sizeB - self.skip

1502 if n > 0 and n != (self.knots * self.u2B):

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

1504

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

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

1507 '''

1508 # flon offset for both west and east

1509 f = 360 if west < 0 else 0

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

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

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

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

1514 s += 1 # s > n

1515 e += 1 # e > w

1516

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

1518 # handle special cases

1519 if (s - n) > hi:

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

1521 if (e - w) > wi:

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

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

1524 return g # use entire geoid as-is

1525

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

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

1528

1529 if e > wi > w: # wrap around

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

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

1532 elif e > w:

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

1534 else:

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

1536

1537 # convert to bytes

1538 r *= self.u2B

1539 p *= self.u2B

1540 q = wi * self.u2B # stride

1541 # number of rows and cols to skip from

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

1543 z = self.skip + (n * wi + w) * self.u2B

1544 # sanity check

1545 if r < 2 or p < 0 or q < 2 or z < self.skip \

1546 or z > self.sizeB:

1547 raise self._Errorf(_format_, _assert_, (r, p, q, z))

1548

1549 # can't use _BytesIO since numpy

1550 # needs .fileno attr in .fromfile

1551 t, c = 0, self._tmpfile()

1552 # reading (s - n) rows, forward

1553 for y in range(n, s): # PYCHOK y unused

1554 g.seek(z, _MODS.os.SEEK_SET)

1555 # Python 2 tmpfile.write returns None

1556 t += c.write(g.read(r)) or r

1557 if p: # wrap around to start of row

1558 g.seek(-q, _MODS.os.SEEK_CUR)

1559 # assert(g.tell() == (z - w * self.u2B))

1560 # Python 2 tmpfile.write returns None

1561 t += c.write(g.read(p)) or p

1562 z += q

1563 c.flush()

1564 g.close()

1565

1566 s -= n # nlat

1567 e -= w # nlon

1568 k = s * e # knots

1569 z = k * self.u2B

1570 if t != z:

1571 raise self._Errorf('%s(%s != %s) %s', _assert_, t, z, self)

1572

1573 # update the _Gpars accordingly, note attributes

1574 # .dlat, .dlon, .rlat and .rlon remain unchanged

1575 self.slat += n * self.dlat

1576 self.wlon += w * self.dlon

1577 self.nlat = s

1578 self.nlon = e

1579 self.flon = self.glon = f

1580

1581 self.crop4 = south, west, north, east

1582 self.knots = k

1583 self.skip = 0 # no header lines in c

1584

1585 c.seek(0, _MODS.os.SEEK_SET)

1586 # c = open(c.name, _rb_) # reopen for numpy 1.8.0-

1587 return c

1588

1589 def _Errorf(self, fmt, *args): # PYCHOK no cover

1590 t = fmt % args

1591 e = self.pgm or NN

1592 if e:

1593 t = _SPACE_(t, _in_, repr(e))

1594 return PGMError(t)

1595

1596 def _lle2yx2(self, lat, lon, flon):

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

1598 # with .dlat decreasing from 90N .slat

1599 lat -= self.slat

1600 lon += flon - self.wlon

1601 return (min(self.nlat - 1, max(0, int(lat * self.rlat))),

1602 max(0, int(lon * self.rlon)))

1603

1604 def _tmpfile(self):

1605 # create a tmpfile to hold the cropped geoid grid

1606 try:

1607 from tempfile import NamedTemporaryFile as tmpfile

1608 except ImportError: # Python 2.7.16-

1609 from _MODS.os import tmpfile

1610 t = _MODS.os.path.basename(self.pgm)

1611 t = _MODS.os.path.splitext(t)[0]

1612 f = tmpfile(mode='w+b', prefix=t or 'egm')

1613 f.seek(0, _MODS.os.SEEK_SET) # force overwrite

1614 return f

1615

1616 @Property_RO

1617 def pgm(self):

1618 '''Get the geoid file name (C{str}).

1619 '''

1620 return self.name

1621

1622

1623 class PGMError(GeoidError):

1624 '''Issue parsing or cropping an C{egm*.pgm} geoid dataset.

1625 '''

1626 pass

1627

1628

1629 def egmGeoidHeights(GeoidHeights_dat):

1630 '''Generate geoid U{egm*.pgm<https://GeographicLib.SourceForge.io/

1631 C++/doc/geoid.html#geoidinst>} height tests from U{GeoidHeights.dat

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

1633 U{Test data for Geoids<https://GeographicLib.SourceForge.io/C++/doc/

1634 geoid.html#testgeoid>}.

1635

1636 @arg GeoidHeights_dat: The un-gz-ed C{GeoidHeights.dat} file

1637 (C{str} or C{file} handle).

1638

1639 @return: For each test, yield a L{GeoidHeight5Tuple}C{(lat, lon,

1640 egm84, egm96, egm2008)}.

1641

1642 @raise GeoidError: Invalid B{C{GeoidHeights_dat}}.

1643

1644 @note: Function L{egmGeoidHeights} is used to test the geoids

1645 L{GeoidKarney} and L{GeoidPGM}, see PyGeodesy module

1646 C{test/testGeoids.py}.

1647 '''

1648 dat = GeoidHeights_dat

1649 if isinstance(dat, bytes):

1650 dat = _BytesIO(dat)

1651

1652 try:

1653 dat.seek(0, _MODS.os.SEEK_SET) # reset

1654 except AttributeError as x:

1655 raise GeoidError(GeoidHeights_dat=type(dat), cause=x)

1656

1657 for t in dat.readlines():

1658 t = t.strip()

1659 if t and not t.startswith(_bHASH_):

1660 lat, lon, egm84, egm96, egm2008 = map(float, t.split())

1661 while lon > _180_0: # EasternLon to earth lon

1662 lon -= _360_0

1663 yield GeoidHeight5Tuple(lat, lon, egm84, egm96, egm2008)

1664

1665

1666__all__ += _ALL_DOCS(_GeoidBase)

1667

1668 if __name__ == '__main__': # MCCABE 14

1669

1670 from pygeodesy.internals import printf, _secs2str, _versions, _sys

1671 from time import time

1672

1673 _crop = ()

1674 _GeoidEGM = GeoidKarney

1675 _kind = 3

1676

1677 geoids = _sys.argv[1:]

1678 while geoids:

1679 geoid = geoids.pop(0)

1680

1681 if '-crop'.startswith(geoid.lower()):

1682 _crop = 20, -125, 50, -65 # CONUS

1683

1684 elif '-karney'.startswith(geoid.lower()):

1685 _GeoidEGM = GeoidKarney

1686

1687 elif '-kind'.startswith(geoid.lower()):

1688 _kind = int(geoids.pop(0))

1689

1690 elif '-pgm'.startswith(geoid.lower()):

1691 _GeoidEGM = GeoidPGM

1692

1693 elif geoid[-4:].lower() in ('.pgm',):

1694 g = _GeoidEGM(geoid, crop=_crop, kind=_kind)

1695 t = time()

1696 _ = g.highest()

1697 t = _secs2str(time() - t)

1698 printf('%s: %s (%s)', g.toStr(), t, _versions(), nl=1, nt=1)

1699 printf(repr(g.pgm), nt=1)

1700 # <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>:

1701 # The height of the EGM96 geoid at Timbuktu

1702 # echo 16:46:33N 3:00:34W | GeoidEval

1703 # => 28.7068 -0.02e-6 -1.73e-6

1704 # The 1st number is the height of the geoid, the 2nd and

1705 # 3rd are its slopes in northerly and easterly direction

1706 t = 'Timbuktu %s' % (g,)

1707 k = {'egm84-15.pgm': '31.2979',

1708 'egm96-5.pgm': '28.7067',

1709 'egm2008-1.pgm': '28.7880'}.get(g.name.lower(), '28.7880')

1710 ll = _MODS.dms.parseDMS2('16:46:33N', '3:00:34W', sep=':')

1711 for ll in (ll, (16.776, -3.009),):

1712 try:

1713 h, ll = g.height(*ll), fstr(ll, prec=6)

1714 printf('%s.height(%s): %.4F vs %s', t, ll, h, k)

1715 except (GeoidError, RangeError) as x:

1716 printf(_COLONSPACE_(t, str(x)))

1717

1718 elif geoid[-4:].lower() in ('.bin',):

1719 g = GeoidG2012B(geoid, kind=_kind)

1720 printf(g.toStr())

1721

1722 else:

1723 raise GeoidError(grid=repr(geoid))

1724

1725_I = int # PYCHOK unused _I

1726 del _intCs # trash ints cache

1727

1728

1729# <https://GeographicLib.SourceForge.io/cgi-bin/GeoidEval>

1730# _lowerleft = -90, -179, -30.1500 # egm2008-1.pgm

1731# _lowerleft = -90, -179, -29.5350 # egm96-5.pgm

1732# _lowerleft = -90, -179, -29.7120 # egm84-15.pgm

1733

1734# _center = 0, 0, 17.2260 # egm2008-1.pgm

1735# _center = 0, 0, 17.1630 # egm96-5.pgm

1736# _center = 0, 0, 18.3296 # egm84-15.pgm

1737

1738# _upperright = 90, 180, 14.8980 # egm2008-1.pgm

1739# _upperright = 90, 180, 13.6050 # egm96-5.pgm

1740# _upperright = 90, 180, 13.0980 # egm84-15.pgm

1741

1742

1743# % python3.12 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm

1744#

1745# GeoidKarney('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, 147.367, 85.839), lowest(4.7, 78.767, -106.911): 204.334 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)

1746#

1747# _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'

1748#

1749# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880

1750# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880

1751#

1752# GeoidKarney('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, 148.75, 81.33), lowest(4.75, 79.25, -107.34): 1.007 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)

1753#

1754# _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'

1755#

1756# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979

1757# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979

1758#

1759# GeoidKarney('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, 0.0, 17.163), highest(-8.167, 147.25, 85.422), lowest(4.667, 78.833, -107.043): 8.509 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)

1760#

1761# _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'

1762#

1763# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067

1764# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067

1765

1766

1767# % python3.8 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm

1768#

1769# GeoidKarney('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, 147.367, 85.839), lowest(4.7, 78.767, -106.911): 353.050 ms (pygeodesy 24.8.24 Python 3.8.10 64bit arm64_x86_64 macOS 10.16)

1770#

1771# _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'

1772#

1773# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880

1774# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880

1775#

1776# GeoidKarney('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, 148.75, 81.33), lowest(4.75, 79.25, -107.34): 1.727 ms (pygeodesy 24.8.24 Python 3.8.10 64bit arm64_x86_64 macOS 10.16)

1777#

1778# _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'

1779#

1780# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979

1781# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979

1782#

1783# GeoidKarney('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, 0.0, 17.163), highest(-8.167, 147.25, 85.422), lowest(4.667, 78.833, -107.043): 14.807 ms (pygeodesy 24.8.24 Python 3.8.10 64bit arm64_x86_64 macOS 10.16)

1784#

1785# _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'

1786#

1787# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067

1788# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067

1789

1790

1791# % python2 -m pygeodesy.geoids -Karney ../testGeoids/egm*.pgm

1792#

1793# GeoidKarney('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, 147.367, 85.839), lowest(4.7, 78.767, -106.911): 283.362 ms (pygeodesy 24.8.24 Python 2.7.18 64bit arm64_x86_64 macOS 10.16)

1794#

1795# _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'

1796#

1797# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880

1798# Timbuktu GeoidKarney('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880

1799#

1800# GeoidKarney('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, 148.75, 81.33), lowest(4.75, 79.25, -107.34): 1.378 ms (pygeodesy 24.8.24 Python 2.7.18 64bit arm64_x86_64 macOS 10.16)

1801#

1802# _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'

1803#

1804# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979

1805# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979

1806#

1807# GeoidKarney('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, 0.0, 17.163), highest(-8.167, 147.25, 85.422), lowest(4.667, 78.833, -107.043): 11.612 ms (pygeodesy 24.8.24 Python 2.7.18 64bit arm64_x86_64 macOS 10.16)

1808#

1809# _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'

1810#

1811# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067

1812# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067

1813

1814

1815# % python3.12 -m pygeodesy.geoids -PGM ../testGeoids/egm*.pgm

1816#

1817# GeoidPGM('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, -32.633, 85.839), lowest(4.683, -101.25, -106.911): 543.148 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)

1818#

1819# _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'

1820#

1821# Timbuktu GeoidPGM('egm2008-1.pgm').height(16.775833, -3.009444): 28.7881 vs 28.7880

1822# Timbuktu GeoidPGM('egm2008-1.pgm').height(16.776, -3.009): 28.7880 vs 28.7880

1823#

1824# GeoidPGM('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, -31.25, 81.33), lowest(4.75, -100.75, -107.34): 1.762 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)

1825#

1826# _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'

1827#

1828# Timbuktu GeoidPGM('egm84-15.pgm').height(16.775833, -3.009444): 31.2979 vs 31.2979

1829# Timbuktu GeoidPGM('egm84-15.pgm').height(16.776, -3.009): 31.2975 vs 31.2979

1830#

1831# GeoidPGM('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, -0.0, 17.179), highest(-8.167, -32.75, 85.422), lowest(4.667, -101.167, -107.043): 12.594 ms (pygeodesy 24.8.24 Python 3.12.5 64bit arm64 macOS 14.6.1)

1832#

1833# _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'

1834#

1835# Timbuktu GeoidPGM('egm96-5.pgm').height(16.775833, -3.009444): 28.7065 vs 28.7067

1836# Timbuktu GeoidPGM('egm96-5.pgm').height(16.776, -3.009): 28.7064 vs 28.7067

1837

1838# **) MIT License

1839#

1840# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1841#

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

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

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

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

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

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

1848#

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

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

1851#

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

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

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

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

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

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

1858# OTHER DEALINGS IN THE SOFTWARE.