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

67 C{scipy} can be thrown as L{SciPyWarning} exceptions, provided Python

68 C{warnings} are filtered 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.

75 interpolate.RectBivariateSpline.html>} and U{interp2d<https://docs.SciPy.org/doc/scipy/

76 reference/generated/scipy.interpolate.interp2d.html>}, functions L{elevations.elevation2}

77 and L{elevations.geoidHeight2}, U{I{Ellispoid vs Orthometric Elevations}<https://

78 www.YouTube.com/watch?v=dX6a6kCk3Po>} and U{I{Pitfalls Related to Ellipsoid Height

79 and Height Above Mean Sea Level (AMSL)}<https://Wiki.ROS.org/mavros#mavros.2FPlugins.

80 Avoiding_Pitfalls_Related_to_Ellipsoid_Height_and_Height_Above_Mean_Sea_Level>}.

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, _0_0, _1_0, _180_0, _360_0

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

88# from pygeodesy.dms import parseDMS2 # _MODS

89 from pygeodesy.errors import _incompatible, LenError, RangeError, SciPyError, \

90 _SciPyIssue, _xkwds_pop2

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

92# from pygoedesy.formy import heightOrthometric # _MODS

93 from pygeodesy.heights import _as_llis2, _ascalar, _height_called, HeightError, \

94 _HeightsBase, _ellipsoidal_datum, _Wrap

95# from pygeodesy.internals import _version2 # _MODS

96 from pygeodesy.interns import MISSING, NN, _4_, _COLONSPACE_, _COMMASPACE_, \

97 _cubic_, _E_, _height_, _in_, _kind_, _lat_, \

98 _linear_, _lon_, _mean_, _N_, _n_a_, _numpy_, \

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

100 _stdev_, _tbd_, _W_, _width_

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 deprecated_method, Property_RO, property_RO

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

110 import os.path as _os_path

111 from os import SEEK_CUR as _SEEK_CUR, SEEK_SET as _SEEK_SET

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.06.11'

122

123_assert_ = 'assert'

124_bHASH_ = b'#'

125_endian_ = 'endian'

126_format_ = '%s %r'

127_header_ = 'header'

128_intCs = {} # cache int value

129_interp2d_ks = {-2: _linear_,

130 -3: _cubic_,

131 -5: 'quintic'}

132_lli_ = 'lli'

133_non_increasing_ = 'non-increasing'

134_rb_ = 'rb'

135_supported_ = 'supported'

136

137

138 class _GeoidBase(_HeightsBase):

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

140 '''

141 _cropped = None

142# _datum = _WGS84 # from _HeightsBase

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

144 _endian = _tbd_

145 _geoid = _n_a_

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

147 _interp2d = None # interp2d interpolation

148 _kind = 3 # order for interp2d, RectBivariateSpline

149# _kmin = 2 # min number of knots

150 _knots = 0 # nlat * nlon

151 _mean = None # fixed in GeoidKarney

152# _name = NN # _Named

153 _nBytes = 0 # numpy size in bytes, float64

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

155 _sizeB = 0 # geoid file size in bytes

156 _smooth = 0 # used only for RectBivariateSpline

157 _stdev = None # fixed in GeoidKarney

158 _u2B = 0 # np.itemsize or undefined

159

160 _lat_d = _0_0 # increment, +tive

161 _lat_lo = _0_0 # lower lat, south

162 _lat_hi = _0_0 # upper lat, noth

163 _lon_d = _0_0 # increment, +tive

164 _lon_lo = _0_0 # left lon, west

165 _lon_hi = _0_0 # right lon, east

166 _lon_of = _0_0 # forward lon offset

167 _lon_og = _0_0 # reverse lon offset

168

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

170 _yx_hits = None # cache hits, ala Karney

171

172 def __init__(self, hs, p):

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

174 and several internal geoid attributes.

175

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

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

178 other geoid parameters (C{INTERNAL}).

179

180 @raise GeoidError: Incompatible grid B{C{hs}} shape or

181 invalid B{C{kind}}.

182

183 @raise LenError: Mismatch grid B{C{hs}} axis.

184

185 @raise SciPyError: A C{scipy.interpolate.inter2d} or

186 C{-.RectBivariateSpline} issue.

187

188 @raise SciPyWarning: A C{scipy.interpolate.inter2d} or

189 C{-.RectBivariateSpline} warning as

190 exception.

191

192 @note: C{scipy.interpolate.interp2d} has been C{DEPRECATED},

193 specify keyword argument C{B{kind}=1..5} to use

194 C{scipy.interpolate.RectBivariateSpline}.

195 '''

196 spi = self.scipy_interpolate

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

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

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

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

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

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

203

204 # both axes and bounding box

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

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

207

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

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

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

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

212 k = self.kind

213 if k in _interp2d_ks: # .interp2d DEPRECATED since scipy 1.10

214 if self._scipy_version() < (1, 10):

215 self._interp2d = spi.interp2d(xs, ys, hs, kind=_interp2d_ks[k])

216 else: # call and overwrite the DEPRECATED .interp2d

217 self._interp2d = self._interp2d(xs, ys, hs, k)

218 elif 1 <= k <= 5:

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

220 s=self._smooth).ev

221 else:

222 raise GeoidError(kind=k)

223

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

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

226 self._knots = p.knots # grid knots

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

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

229 # shrink the box by 1 unit on every side

230 # bb += self._lat_d, -self._lat_d, self._lon_d, -self._lon_d

231 self._lat_lo = float(bb[0])

232 self._lat_hi = float(bb[1])

233 self._lon_lo = float(bb[2] - p.glon)

234 self._lon_hi = float(bb[3] - p.glon)

235

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

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

238

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

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

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

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

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

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

245

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

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

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

249

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

251 invalid B{C{lli}} or the C{egm*.pgm}

252 geoid file is closed.

253

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

255 or longitude range.

256

257 @raise SciPyError: A C{scipy.interpolate.inter2d} or

258 C{-.RectBivariateSpline} issue.

259

260 @raise SciPyWarning: A C{scipy.interpolate.inter2d} or

261 C{-.RectBivariateSpline} warning as

262 exception.

263

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

265 location must have an ellipsoid C{height} or C{h}

266 attribute, otherwise C{height=0} is used.

267

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

269 '''

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

271

272 def __enter__(self):

273 '''Open context.

274 '''

275 return self

276

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

278 '''Close context.

279 '''

280 self.close()

281 # return None # XXX False

282

283 def __repr__(self):

284 return self.toStr()

285

286 def __str__(self):

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

288

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

290 # handle __call__

291 _H = _MODS.formy.heightOrthometric if H else None

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

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

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

295 try:

296 for i, lli in enumerate(llis):

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

298 # orthometric or geoid height

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

300 return _as(hs)

301

302 except (GeoidError, RangeError) as x:

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

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

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

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

307 except Exception as x:

308 if scipy and self.scipy:

309 raise _SciPyIssue(x)

310 else:

311 raise

312

313 @Property_RO

314 def _center(self):

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

316 '''

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

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

319

320 def center(self, LatLon=None):

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

322

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

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

325

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

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

328 longitude and geoid height of the center grid location.

329 '''

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

331

332 def close(self):

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

334 '''

335 if not self.closed:

336 self._egm.close()

337 self._egm = None

338

339 @property_RO

340 def closed(self):

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

342 '''

343 return self._egm is None

344

345 @Property_RO

346 def cropped(self):

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

348 '''

349 return self._cropped

350

351 @Property_RO

352 def dtype(self):

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

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

355 '''

356 return self._hs_y_x.dtype

357

358 @Property_RO

359 def endian(self):

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

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

362 '''

363 return self._endian

364

365 def _ev(self, y, x): # PYCHOK expected

366 # only used for .interpolate.interp2d, but

367 # overwritten for .RectBivariateSpline,

368 # note (y, x) must be flipped!

369 return self._interp2d(x, y)

370

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

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

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

374 if m != n:

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

376 if d < 0:

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

378 for i in range(1, m):

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

380 if e < EPS: # non-increasing axis

381 i = Fmt.INDEX(name, i)

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

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

384

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

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

387 self._notOverloaded(lat, lon)

388

389 def _gyx2g2(self, y, x):

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

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

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

393

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

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

396

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

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

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

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

401

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

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

404

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

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

407

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

409 geoid's lat- or longitude range.

410

411 @raise SciPyError: A C{scipy.interpolate.inter2d} or

412 C{-.RectBivariateSpline} issue.

413

414 @raise SciPyWarning: A C{scipy.interpolate.inter2d} or

415 C{-.RectBivariateSpline} warning as

416 exception.

417 '''

418 return _height_called(self, lats, lons, Error=GeoidError, **wrap)

419

420 def _hGeoid(self, lat, lon):

421 out = self.outside(lat, lon)

422 if out:

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

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

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

426

427 @Property_RO

428 def _highest(self):

429 '''(INTERNAL) Cache for L{highest} method.

430 '''

431 return self._llh3minmax(True)

432

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

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

435

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

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

438

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

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

441 longitude and geoid height of the highest grid location.

442 '''

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

444

445 @Property_RO

446 def hits(self):

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

448 '''

449 return self._yx_hits

450

451 @deprecated_method

452 def _interp2d(self, xs, ys, hs=(), k=0): # overwritten in .__init__ above

453 '''DEPRECATED on 23.01.06, use keyword argument C{B{kind}=1..5}.'''

454 # assert k in _interp2d_ks # and len(hs) == len(xs) == len(ys)

455 try:

456 return self.scipy_interpolate.interp2d(xs, ys, hs, kind=_interp2d_ks[k])

457 except AttributeError as x:

458 raise SciPyError(interp2d=MISSING, kind=k, cause=x)

459

460 @Property_RO

461 def kind(self):

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

463 '''

464 return self._kind

465

466 @Property_RO

467 def knots(self):

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

469 '''

470 return self._knots

471

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

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

474 self._notOverloaded(lat, lon)

475

476 @property_RO

477 def _LL3T(self):

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

479 '''

480 _GeoidBase._LL3T = T = _MODS.namedTuples.LatLon3Tuple # overwrite poperty_RO

481 return T

482

483 def _llh3(self, lat, lon):

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

485

486 def _llh3LL(self, llh, LatLon):

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

488

489 def _llh3minmax(self, highest=True, *unused):

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

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

492 # numpy.argmin.html#numpy.argmin>

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

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

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

496

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

498 # numpy.fromfile, like .frombuffer

499 g.seek(offset, _SEEK_SET)

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

501

502 @Property_RO

503 def _lowerleft(self):

504 '''(INTERNAL) Cache for L{lowerleft}.

505 '''

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

507

508 def lowerleft(self, LatLon=None):

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

510

511 @kwarg LatLon: Optional class to return the location

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

513

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

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

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

517 '''

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

519

520 @Property_RO

521 def _loweright(self):

522 '''(INTERNAL) Cache for L{loweright}.

523 '''

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

525

526 def loweright(self, LatLon=None):

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

528

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

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

531

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

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

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

535 '''

536

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

538

539 lowerright = loweright # synonymous

540

541 @Property_RO

542 def _lowest(self):

543 '''(INTERNAL) Cache for L{lowest}.

544 '''

545 return self._llh3minmax(False)

546

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

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

549

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

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

552

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

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

555 longitude and geoid height of the lowest grid location.

556 '''

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

558

559 @Property_RO

560 def mean(self):

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

562 '''

563 if self._mean is None: # see GeoidKarney

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

565 return self._mean

566

567 @property_RO

568 def name(self):

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

570 '''

571 return _HeightsBase.name.fget(self) or self._geoid # recursion

572

573 @Property_RO

574 def nBytes(self):

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

576 '''

577 return self._nBytes

578

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

580 # open the geoid file

581 try:

582 self._geoid = _os_path.basename(geoid)

583 self._sizeB = _os_path.getsize(geoid)

584 g = open(geoid, _rb_)

585 except (IOError, OSError) as x:

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

587

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

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

590 self._kind = int(kind)

591 if name:

592 _HeightsBase.name.fset(self, name) # rename

593 if smooth:

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

595

596 return g

597

598 def outside(self, lat, lon):

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

600 lat-/longitude or crop range.

601

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

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

604

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

606 empty C{str} if inside.

607 '''

608 return (_S_ if lat < self._lat_lo else

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

610 (_W_ if lon < self._lon_lo else

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

612

613 @Property_RO

614 def pgm(self):

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

616 '''

617 return self._pgm

618

619 @Property_RO

620 def sizeB(self):

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

622 '''

623 return self._sizeB

624

625 @Property_RO

626 def smooth(self):

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

628 '''

629 return self._smooth

630

631 @Property_RO

632 def stdev(self):

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

634 '''

635 if self._stdev is None: # see GeoidKarney

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

637 return self._stdev

638

639 def _swne(self, crop):

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

641 try:

642 if len(crop) == 2:

643 try: # sw, ne LatLons

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

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

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

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

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

649 swne = crop

650 if len(swne) == 4:

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

652 if -90 <= s <= (n - _1_0) <= 89 and \

653 -180 <= w <= (e - _1_0) <= 179:

654 return s, w, n, e

655 except (IndexError, TypeError, ValueError):

656 pass

657 raise GeoidError(crop=crop)

658

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

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

661

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

663 default). Trailing zero decimals are stripped

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

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

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

667

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

669 '''

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

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

672 (self.lowerleft, self.upperright,

673 self.center,

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

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

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

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

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

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

680

681 @Property_RO

682 def u2B(self):

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

684 '''

685 return self._u2B

686

687 @Property_RO

688 def _upperleft(self):

689 '''(INTERNAL) Cache for method L{upperleft}.

690 '''

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

692

693 def upperleft(self, LatLon=None):

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

695

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

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

698

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

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

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

702 '''

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

704

705 @Property_RO

706 def _upperright(self):

707 '''(INTERNAL) Cache for method L{upperright}.

708 '''

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

710

711 def upperright(self, LatLon=None):

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

713

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

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

716

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

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

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

720 '''

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

722

723

724 class GeoidError(HeightError):

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

726 '''

727 pass

728

729

730 class GeoidG2012B(_GeoidBase):

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

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

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

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

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

736 U{Guam and Northern Mariana Islands

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

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

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

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

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

742 scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>}

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

744 scipy.interpolate.interp2d.html>} interpolation.

745

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

747 C{g2012b*.bin} grid files.

748 '''

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

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

751 '''New L{GeoidG2012B} interpolator.

752

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

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

755 or L{a_f2Tuple}), default C{WGS84}.

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

757 U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/

758 reference/generated/scipy.interpolate.RectBivariateSpline.html>},

759 -2 for U{interp2d linear<https://docs.SciPy.org/doc/scipy/

760 reference/generated/scipy.interpolate.interp2d.html>}, -3

761 for C{interp2d cubic} or -5 for C{interp2d quintic}.

762 @kwarg smooth: Smoothing factor for U{RectBivariateSpline

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

764 scipy.interpolate.RectBivariateSpline.html>}

765 only (C{int}).

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

767 keyword argument C{B{crop}}, use C{B{crop}=None} to ignore.

768

769 @raise GeoidError: G2012B grid file B{C{g2012b_bin}} issue or invalid

770 B{C{crop}}, B{C{kind}} or B{C{smooth}}.

771

772 @raise ImportError: Package C{numpy} or C{scipy} not found or not

773 installed.

774

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

776

777 @raise SciPyError: A C{RectBivariateSpline} or C{inter2d} issue.

778

779 @raise SciPyWarning: A C{RectBivariateSpline} or C{inter2d}

780 warning as exception.

781

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

783

784 @note: C{scipy.interpolate.interp2d} has been C{DEPRECATED}, specify

785 C{B{kind}=1..5} for C{scipy.interpolate.RectBivariateSpline}.

786 '''

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

788 if crop is not None:

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

790

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

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

793

794 try:

795 p = _Gpars()

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

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

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

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

800 # skip 4xf8, get 3xi4

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

802 if ien == 1: # correct endian

803 p.knots = p.nlat * p.nlon

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

805 and 1 < p.nlon < n:

806 self._endian = en_+'f4'

807 break

808 else: # couldn't validate endian

809 raise GeoidError(_endian_)

810

811 # get the first 4xf8

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

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

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

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

816 _GeoidBase.__init__(self, hs, p)

817

818 except Exception as x:

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

820 finally:

821 g.close()

822

823 def _g2ll2(self, lat, lon):

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

825 return lat, lon

826

827 def _ll2g2(self, lat, lon):

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

829 return lat, lon

830

831 if _FOR_DOCS:

832 __call__ = _GeoidBase.__call__

833 height = _GeoidBase.height

834

835

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

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

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

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

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

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

842 '''

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

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

845

846

847 def _I(i):

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

849 '''

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

851

852

853 def _T(*cs):

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

855 '''

856 return map1(_I, *cs)

857

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

859

860

861 class GeoidKarney(_GeoidBase):

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

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

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

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

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

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

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

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

870

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

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

873 datasets.

874 '''

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

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

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

878 _T0s12,

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

880 _T0s12,

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

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

883 _T0s12,

884 _T(0, 0, 0, 0, 0, 0, 93, -93, -93, 93, 0, 0),

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

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

887

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

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

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

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

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

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

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

895 _T(0, 0, -60, 60, 60, -60, 60, -60, -60, 60, 0, 0),

896 _T(-60, 60, 0, 60, -60, 0, 0, 60, -60, 0, -60, 60),

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

898

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

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

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

902 _T(0, 0, 93, -93, -93, 93, 0, 0, 0, 0, 0, 0),

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

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

905 _T0s12,

906 _T(0, 0, -93, 93, 93, -93, 0, 0, 0, 0, 0, 0),

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

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

909

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

911 _T(1, 0),

912 _T(0, 1),

913 _T(1, 1))

914

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

916 _T( 1, -1),

917 _T(-1, 0),

918 _T( 0, 0),

919 _T( 1, 0),

920 _T( 2, 0),

921 _T(-1, 1),

922 _T( 0, 1),

923 _T( 1, 1),

924 _T( 2, 1),

925 _T( 0, 2),

926 _T( 1, 2))

927

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

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

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

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

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

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

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

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

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

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

938 _nBytes = None # not applicable

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

940 _smooth = None # not applicable

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

942 _u2B = _calcsize(_endian) # pixelsize_ in bytes

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

944 _Ru2B = 0 # row of pixelsize_s in bytes

945 _yxH = () # cache (y, x) indices

946 _yxHt = () # cached 4- or 10-tuple for _ev2H resp. _ev3H

947 _yx_hits = 0 # cache hits

948

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

950 kind=3, **name_smooth):

951 '''New L{GeoidKarney} interpolator.

952

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

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

955 note below.

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

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

958 east)}) or 2, in C{degrees90} lat- and C{degrees180}

959 longitudes or a 2-tuple (C{LatLonSW, LatLonNE}) of

960 C{LatLon} instances.

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

962 or L{a_f2Tuple}), default C{WGS84}.

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

964 for C{cubic}.

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

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

967

968 @raise GeoidError: EGM dataset B{C{egm_pgm}} issue or invalid

969 B{C{crop}}, B{C{kind}} or B{C{smooth}}.

970

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

972

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

974

975 @note: Geoid file B{C{egm_pgm}} remains open and must be closed

976 by calling the C{close} method or by using this instance

977 in a C{with B{GeoidKarney}(...) as ...} context.

978 '''

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

980 if smooth is not None:

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

982

983 if kind in (2,):

984 self._evH = self._ev2H

985 elif kind not in (3,):

986 raise GeoidError(kind=kind)

987

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

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

990

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

992 self._Ru2B = _calcsize(self._Rendian)

993

994 self._knots = p.knots # grid knots

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

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

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

998 self._lat_lo, self._lon_lo, \

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

1000 self._cropped = True if crop else False

1001

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

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

1004

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

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

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

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

1009 return the I{orthometric} height instead of the

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

1011

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

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

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

1015

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

1017 invalid B{C{lli}} or the C{egm*.pgm}

1018 geoid file is closed.

1019

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

1021 or longitude range.

1022

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

1024 location must have an ellipsoid C{height} or C{h}

1025 attribute, otherwise C{height=0} is used.

1026

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

1028 '''

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

1030

1031 def _c0c3v(self, y, x):

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

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

1034 p = self._pgm

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

1036 # read 4x4 ushorts, drop the 4 corners

1037 g = self._egm

1038 e = self._4endian

1039 n = self._4u2B

1040 R = self._Ru2B

1041

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

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

1044 b += R

1045 g.seek(b, _SEEK_SET)

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

1047 b += R

1048 g.seek(b, _SEEK_SET)

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

1050 b += R

1051 g.seek(b, _SEEK_SET)

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

1053 j = 1

1054

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

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

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

1058

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

1060

1061 @Property_RO

1062 def dtype(self):

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

1064 '''

1065 return 'ushort'

1066

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

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

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

1070 y, x = int(floor(fy)), int(floor(fx))

1071 fy -= y

1072 fx -= x

1073 H = self._evH(fy, fx, y, x) # ._ev3H or ._ev2H

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

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

1076 return H.fsum()

1077

1078 def _ev2H(self, fy, fx, *yx):

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

1080 if self._yxH == yx:

1081 t = self._yxHt

1082 self._yx_hits += 1

1083 else:

1084 y, x = self._yxH = yx

1085 self._yxHt = t = self._raws(y, x, GeoidKarney._BT)

1086 v = _1_0, -fx, fx

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

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

1089 return H

1090

1091 def _ev3H(self, fy, fx, *yx):

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

1093 if self._yxH == yx:

1094 t = self._yxHt

1095 self._yx_hits += 1

1096 else:

1097 self._yxH = yx

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

1099 t = [fdot(v, *c3[i]) / c0 for i in range(self._nterms)]

1100 self._yxHt = t = tuple(t)

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

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

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

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

1105 v = _1_0, fx, fy

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

1107 H *= fy

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

1109 H *= fy

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

1111 return H

1112

1113 _evH = _ev3H # overriden for kind == 2

1114

1115 def _g2ll2(self, lat, lon):

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

1117 while lon > _180_0:

1118 lon -= _360_0

1119 return lat, lon

1120

1121 def _g2yx2(self, lat, lon):

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

1123 p = self._pgm

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

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

1126

1127 def _gyx2g2(self, y, x):

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

1129 p = self._pgm

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

1131

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

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

1134

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

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

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

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

1139

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

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

1142

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

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

1145 geoid file is closed.

1146

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

1148 geoid's lat- or longitude range.

1149 '''

1150 return _height_called(self, lats, lons, Error=GeoidError, **wrap)

1151

1152 @Property_RO

1153 def _highest_ltd(self):

1154 '''(INTERNAL) Cache for L{highest} mesthod.

1155 '''

1156 return self._llh3minmax(True, -12, -4)

1157

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

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

1160

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

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

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

1164

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

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

1167 longitude and geoid height of the highest grid location.

1168 '''

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

1170 return self._llh3LL(llh, LatLon)

1171

1172 def _lat2y2(self, lat2):

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

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

1175 for lat in lat2:

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

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

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

1179

1180 def _ll2g2(self, lat, lon):

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

1182 while lon < 0:

1183 lon += _360_0

1184 return lat, lon

1185

1186 def _llh3minmax(self, highest=True, *lat2):

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

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

1189 y = x = 0

1190 h = self._raw(y, x)

1191 if highest:

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

1193 m = max(r)

1194 if m > h:

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

1196 else: # lowest

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

1198 m = min(r)

1199 if m < h:

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

1201 h *= self._pgm.Scale

1202 h += self._pgm.Offset

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

1204

1205 @Property_RO

1206 def _lowest_ltd(self):

1207 '''(INTERNAL) Cache for L{lowest}.

1208 '''

1209 return self._llh3minmax(False, 0, 8)

1210

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

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

1213

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

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

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

1217

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

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

1220 longitude and geoid height of the lowest grid location.

1221 '''

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

1223 return self._llh3LL(llh, LatLon)

1224

1225 def _raw(self, y, x):

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

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

1228 p = self._pgm

1229 if x < 0:

1230 x += p.nlon

1231 elif x >= p.nlon:

1232 x -= p.nlon

1233 h = p.nlon // 2

1234 if y < 0:

1235 y = -y

1236 elif y >= p.nlat:

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

1238 else:

1239 h = 0

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

1241 self._seek(y, x)

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

1243 return h[0]

1244

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

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

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

1248

1249 def _raw2(self, *lat2):

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

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

1252 p = self._pgm

1253 g = self._egm

1254 e = self._Rendian

1255 n = self._Ru2B

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

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

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

1259 for y in range(s, t):

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

1261

1262 def _seek(self, y, x):

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

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

1265 if g:

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

1267 g.seek(b, _SEEK_SET)

1268 return b # position

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

1270

1271

1272 class GeoidPGM(_GeoidBase):

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

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

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

1276 C++/doc/geoid.html#geoidinst>} datasets but based on C{SciPy}

1277 U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/reference/

1278 generated/scipy.interpolate.RectBivariateSpline.html>} or

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

1280 scipy.interpolate.interp2d.html>} interpolation.

1281

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

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

1284 datasets. However, unless cropped, an entire C{egm*.pgm} dataset

1285 is loaded into the C{SciPy} U{RectBivariateSpline<https://docs.SciPy.org/

1286 doc/scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html>}

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

1288 scipy.interpolate.interp2d.html>} interpolator and converted from

1289 2-byte C{int} to 8-byte C{dtype float64}. Therefore, internal memory

1290 usage is 4x the U{egm*.pgm<https://GeographicLib.SourceForge.io/C++/doc/

1291 geoid.html#geoidinst>} file size and may exceed the available memory,

1292 especially with 32-bit Python, see properties C{.nBytes} and C{.sizeB}.

1293 '''

1294 _cropped = False

1295 _endian = '>u2'

1296

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

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

1299 '''New L{GeoidPGM} interpolator.

1300

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

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

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

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

1305 in C{degrees90} lat- and C{degrees180} longitudes or a

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

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

1308 or L{a_f2Tuple}), default C{WGS84}.

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

1310 U{RectBivariateSpline<https://docs.SciPy.org/doc/scipy/

1311 reference/generated/scipy.interpolate.RectBivariateSpline.html>},

1312 -2 for U{interp2d linear<https://docs.SciPy.org/doc/scipy/

1313 reference/generated/scipy.interpolate.interp2d.html>}, -3

1314 for C{interp2d cubic} or -5 for C{interp2d quintic}.

1315 @kwarg smooth: Smoothing factor for U{RectBivariateSpline

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

1317 scipy.interpolate.RectBivariateSpline.html>}

1318 only (C{int}).

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

1320

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

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

1323

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

1325

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

1327

1328 @raise SciPyError: A C{RectBivariateSpline} or C{inter2d} issue.

1329

1330 @raise SciPyWarning: A C{RectBivariateSpline} or C{inter2d}

1331 warning as exception.

1332

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

1334

1335 @note: C{scipy.interpolate.interp2d} has been C{DEPRECATED}, specify

1336 C{B{kind}=1..5} for C{scipy.interpolate.RectBivariateSpline}.

1337

1338 @note: The U{GeographicLib egm*.pgm<https://GeographicLib.SourceForge.io/

1339 C++/doc/geoid.html#geoidinst>} file sizes are based on a 2-byte

1340 C{int} height converted to 8-byte C{dtype float64} for C{scipy}

1341 interpolators. Therefore, internal memory usage is 4 times the

1342 C{egm*.pgm} file size and may exceed the available memory,

1343 especially with 32-bit Python. To reduce memory usage, set

1344 keyword argument B{C{crop}} to the region of interest. For example

1345 C{B{crop}=(20, -125, 50, -65)} covers the U{conterminous US<https://

1346 www.NGS.NOAA.gov/GEOID/GEOID12B/maps/GEOID12B_CONUS_grids.png>}

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

1348

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

1350 '''

1351 np = self.numpy

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

1353

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

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

1356 if crop:

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

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

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

1360 self._cropped = True

1361 try:

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

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

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

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

1366 if p.Offset: # offset

1367 hs = p.Offset + hs

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

1369 hs = np.flipud(hs)

1370 _GeoidBase.__init__(self, hs, p)

1371 except Exception as x:

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

1373 finally:

1374 g.close()

1375

1376 def _g2ll2(self, lat, lon):

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

1378 if self._cropped:

1379 lon -= self._lon_of

1380 else:

1381 while lon > _180_0:

1382 lon -= _360_0

1383 return lat, lon

1384

1385 def _ll2g2(self, lat, lon):

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

1387 if self._cropped:

1388 lon += self._lon_of

1389 else:

1390 while lon < 0:

1391 lon += _360_0

1392 return lat, lon

1393

1394 if _FOR_DOCS:

1395 __call__ = _GeoidBase.__call__

1396 height = _GeoidBase.height

1397

1398

1399 class _Gpars(_Named):

1400 '''(INTERNAL) Basic geoid parameters.

1401 '''

1402 # interpolator parameters

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

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

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

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

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

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

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

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

1411

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

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

1414

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

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

1417

1418 def __repr__(self):

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

1420 a in dir(self.__class__)

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

1422 return _COLONSPACE_(self, t)

1423

1424 def __str__(self):

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

1426

1427

1428 class _PGM(_Gpars):

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

1430

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

1432 # Description WGS84 EGM96, 5-minute grid

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

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

1435 # MaxBilinearError 0.140

1436 # RMSBilinearError 0.005

1437 # MaxCubicError 0.003

1438 # RMSCubicError 0.001

1439 # Offset -108

1440 # Scale 0.003

1441 # Origin 90N 0E

1442 # AREA_OR_POINT Point

1443 # Vertical_Datum WGS84

1444 <width> <height>

1445 <pixel>

1446 ...

1447 '''

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

1449 egm = None

1450 glon = 180 # reverse offset, uncropped

1451# pgm = NN # name

1452 sizeB = 0

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

1454

1455 @staticmethod

1456 def _llstr2floats(latlon):

1457 # llstr to (lat, lon) floats

1458 lat, lon = latlon.split()

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

1460

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

1462 AREA_OR_POINT = str

1463 DateTime = str

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

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

1466 MaxBilinearError = float

1467 MaxCubicError = float

1468 Offset = float

1469 Origin = _llstr2floats

1470 Pixel = 0

1471 RMSBilinearError = float

1472 RMSCubicError = float

1473 Scale = float

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

1475 Vertical_Datum = str

1476

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

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

1479 '''

1480 self.name = pgm # geoid file name

1481 if itemsize:

1482 self._u2B = itemsize

1483 if sizeB:

1484 self.sizeB = sizeB

1485

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

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

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

1489

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

1491 try: # ignore empty ones or comments

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

1493 if t.startswith(_bHASH_):

1494 t = t.lstrip(_bHASH_).lstrip()

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

1496 f = getattr(_PGM, a, None)

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

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

1499 elif t:

1500 break

1501 except (TypeError, ValueError):

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

1503 else: # should never get here

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

1505

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

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

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

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

1510 raise ValueError

1511 except (TypeError, ValueError):

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

1513

1514 try: # must be 16 bit pixel height

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

1516 self.Pixel = int(t)

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

1518 raise ValueError

1519 except (TypeError, ValueError):

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

1521

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

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

1524 setattr(self, a, None)

1525

1526 if self.Origin is None:

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

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

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

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

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

1532

1533 self.skip = g.tell()

1534 self.knots = nlat * nlon

1535

1536 self.nlat, self.nlon = nlat, nlon

1537 self.slat, self.wlon = self.Origin

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

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

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

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

1542

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

1544 n = float(self.slat)

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

1546 w = self.wlon - self.glon

1547 e = w + self.dlon * nlon

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

1549

1550 n = self.sizeB - self.skip

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

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

1553

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

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

1556 '''

1557 # flon offset for both west and east

1558 f = 360 if west < 0 else 0

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

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

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

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

1563 s += 1 # s > n

1564 e += 1 # e > w

1565

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

1567 # handle special cases

1568 if (s - n) > hi:

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

1570 if (e - w) > wi:

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

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

1573 return g # use entire geoid as-is

1574

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

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

1577

1578 if e > wi > w: # wrap around

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

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

1581 elif e > w:

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

1583 else:

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

1585

1586 # convert to bytes

1587 r *= self.u2B

1588 p *= self.u2B

1589 q = wi * self.u2B # stride

1590 # number of rows and cols to skip from

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

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

1593 # sanity check

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

1595 or z > self.sizeB:

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

1597

1598 # can't use _BytesIO since numpy

1599 # needs .fileno attr in .fromfile

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

1601 # reading (s - n) rows, forward

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

1603 g.seek(z, _SEEK_SET)

1604 # Python 2 tmpfile.write returns None

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

1606 if p: # wrap around to start of row

1607 g.seek(-q, _SEEK_CUR)

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

1609 # Python 2 tmpfile.write returns None

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

1611 z += q

1612 c.flush()

1613 g.close()

1614

1615 s -= n # nlat

1616 e -= w # nlon

1617 k = s * e # knots

1618 z = k * self.u2B

1619 if t != z:

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

1621

1622 # update the _Gpars accordingly, note attributes

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

1624 self.slat += n * self.dlat

1625 self.wlon += w * self.dlon

1626 self.nlat = s

1627 self.nlon = e

1628 self.flon = self.glon = f

1629

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

1631 self.knots = k

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

1633

1634 c.seek(0, _SEEK_SET)

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

1636 return c

1637

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

1639 t = fmt % args

1640 e = self.pgm or NN

1641 if e:

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

1643 return PGMError(t)

1644

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

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

1647 # with .dlat decreasing from 90N .slat

1648 lat -= self.slat

1649 lon += flon - self.wlon

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

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

1652

1653 def _tmpfile(self):

1654 # create a tmpfile to hold the cropped geoid grid

1655 try:

1656 from tempfile import NamedTemporaryFile as tmpfile

1657 except ImportError: # Python 2.7.16-

1658 from os import tmpfile

1659 t = _os_path.splitext(_os_path.basename(self.pgm))[0]

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

1661 f.seek(0, _SEEK_SET) # force overwrite

1662 return f

1663

1664 @Property_RO

1665 def pgm(self):

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

1667 '''

1668 return self.name

1669

1670

1671 class PGMError(GeoidError):

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

1673 '''

1674 pass

1675

1676

1677 def egmGeoidHeights(GeoidHeights_dat):

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

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

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

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

1682 geoid.html#testgeoid>}.

1683

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

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

1686

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

1688 egm84, egm96, egm2008)}.

1689

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

1691

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

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

1694 C{test/testGeoids.py}.

1695 '''

1696 dat = GeoidHeights_dat

1697 if isinstance(dat, bytes):

1698 dat = _BytesIO(dat)

1699

1700 try:

1701 dat.seek(0, _SEEK_SET) # reset

1702 except AttributeError as x:

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

1704

1705 for t in dat.readlines():

1706 t = t.strip()

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

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

1709 while lon > _180_0: # EasternLon to earth lon

1710 lon -= _360_0

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

1712

1713

1714__all__ += _ALL_DOCS(_GeoidBase)

1715

1716 if __name__ == '__main__':

1717

1718 from pygeodesy.internals import printf, _sys

1719

1720 _crop = ()

1721 _GeoidEGM = GeoidKarney

1722 _kind = 3

1723

1724 geoids = _sys.argv[1:]

1725 while geoids:

1726 geoid = geoids.pop(0)

1727

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

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

1730

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

1732 _GeoidEGM = GeoidKarney

1733

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

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

1736

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

1738 _GeoidEGM = GeoidPGM

1739

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

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

1742 printf(g.toStr(), nt=1, nl=1)

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

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

1745 # The height of the EGM96 geoid at Timbuktu

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

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

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

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

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

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

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

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

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

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

1756 try:

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

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

1759 except (GeoidError, RangeError) as x:

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

1761

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

1763 g = GeoidG2012B(geoid, kind=_kind)

1764 printf(g.toStr())

1765

1766 else:

1767 raise GeoidError(grid=repr(geoid))

1768

1769_I = int # PYCHOK unused _I

1770 del _intCs # trash ints cache

1771

1772# **) MIT License

1773#

1774# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1775#

1776# Permission is hereby granted, free of charge, to any person obtaining a

1777# copy of this software and associated documentation files (the "Software"),

1778# to deal in the Software without restriction, including without limitation

1779# the rights to use, copy, modify, merge, publish, distribute, sublicense,

1780# and/or sell copies of the Software, and to permit persons to whom the

1781# Software is furnished to do so, subject to the following conditions:

1782#

1783# The above copyright notice and this permission notice shall be included

1784# in all copies or substantial portions of the Software.

1785#

1786# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

1787# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

1788# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

1789# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

1790# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

1791# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

1792# OTHER DEALINGS IN THE SOFTWARE.

1793

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

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

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

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

1798

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

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

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

1802

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

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

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

1806

1807

1808# % python3 -m pygeodesy.geoids [-Karney] ../testGeoids/egm*.pgm

1809#

1810# 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)

1811#

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

1813#

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

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

1816#

1817# 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)

1818#

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

1820#

1821# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979

1822# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979

1823#

1824# 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)

1825#

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

1827#

1828# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067

1829# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067

1830

1831

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

1833#

1834# 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)

1835#

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

1837#

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

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

1840#

1841# 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)

1842#

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

1844#

1845# Timbuktu GeoidKarney('egm84-15.pgm').height(16.775833, -3.009444): 31.2983 vs 31.2979

1846# Timbuktu GeoidKarney('egm84-15.pgm').height(16.776, -3.009): 31.2979 vs 31.2979

1847#

1848# 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)

1849#

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

1851#

1852# Timbuktu GeoidKarney('egm96-5.pgm').height(16.775833, -3.009444): 28.7068 vs 28.7067

1853# Timbuktu GeoidKarney('egm96-5.pgm').height(16.776, -3.009): 28.7067 vs 28.7067

1854

1855

1856# % python2 -m pygeodesy.geoids -PGM ../testGeoids/egm*.pgm

1857#

1858# 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)

1859#

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

1861#

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

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

1864#

1865# 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)

1866#

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

1868#

1869# Timbuktu GeoidPGM('egm84-15.pgm').height(16.775833, -3.009444): 31.2979 vs 31.2979

1870# Timbuktu GeoidPGM('egm84-15.pgm').height(16.776, -3.009): 31.2975 vs 31.2979

1871#

1872# 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)

1873#

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

1875#

1876# Timbuktu GeoidPGM('egm96-5.pgm').height(16.775833, -3.009444): 28.7065 vs 28.7067

1877# Timbuktu GeoidPGM('egm96-5.pgm').height(16.776, -3.009): 28.7064 vs 28.7067