Coverage for pygeodesy/geoids.py: 96%

664 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-07-12 13:40 -0400

1 

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

3 

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

5 

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

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

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

9and C{geoid} model files. 

10 

11L{GeoidKarney} is a transcoding of I{Charles Karney}'s C++ class U{Geoid 

12<https://GeographicLib.SourceForge.io/C++/doc/geoid.html>} to pure Python. 

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

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

15require those packages to be installed. 

16 

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

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

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

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

21 

22Typical usage 

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

24 

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

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

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

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

29 

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

31 

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

33arguments to select different interpolation options 

34 

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

36 

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

38 

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

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

41 

42or 

43 

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

45 

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

47 

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

49 

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

51 

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

53 

54or as 2 lists, 2 tuples, etc. 

55 

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

57 

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

59courtesy of SBFRF. 

60 

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

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

63 to be installed. 

64 

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

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

67 C{warnings} are filtered accordingly, see L{SciPyWarning}. 

68 

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

70 U{Geoid height<https://GeographicLib.SourceForge.io/html/geoid.html>} and U{Installing 

71 the Geoid datasets<https://GeographicLib.SourceForge.io/html/geoid.html#geoidinst>}, 

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

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

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

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

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

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

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

79 Avoiding_Pitfalls_Related_to_Ellipsoid_Height_and_Height_Above_Mean_Sea_Level>}. 

80''' 

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

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

83 

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

85from pygeodesy.constants import EPS, _float as _F, _0_0, _1_0, _180_0, _360_0 

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

87from pygeodesy.dms import parseDMS2 

88from pygeodesy.errors import _incompatible, LenError, RangeError, _SciPyIssue 

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

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

91 _HeightsBase, _ellipsoidal_datum, _Wrap 

92from pygeodesy.interns import NN, _4_, _COLONSPACE_, _COMMASPACE_, _cubic_, \ 

93 _DOT_, _E_, _height_, _in_, _kind_, _knots_, \ 

94 _lat_, _linear_, _lon_, _mean_, _N_, _n_a_, _not_, \ 

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

96 _SPACE_, _stdev_, _supported_, _tbd_, _W_, _width_ 

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

98from pygeodesy.named import _Named, _NamedTuple, notOverloaded 

99# from pygeodesy.namedTuples import LatLon3Tuple # _MODS 

100from pygeodesy.props import Property_RO, property_RO 

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

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

103# from pygeodesy.utils import _Wrap # from .heights 

104 

105from math import floor 

106import os.path as _os_path 

107from os import SEEK_CUR as _SEEK_CUR, SEEK_SET as _SEEK_SET 

108from struct import calcsize as _calcsize, unpack as _unpack 

109 

110try: 

111 from StringIO import StringIO as _BytesIO # reads bytes 

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

113 

114except ImportError: # Python 3+ 

115 from io import BytesIO as _BytesIO # PYCHOK expected 

116 

117__all__ = _ALL_LAZY.geoids 

118__version__ = '23.05.26' 

119 

120_assert_ = 'assert' 

121_bHASH_ = b'#' 

122_endian_ = 'endian' 

123_format_ = '%s %r' 

124_header_ = 'header' 

125# temporarily hold a single instance for each int value 

126_intCs = {} 

127_interp2d_ks = {-2: _linear_, 

128 -3: _cubic_, 

129 -5: 'quintic'} 

130_lli_ = 'lli' 

131_non_increasing_ = 'non-increasing' 

132_rb_ = 'rb' 

133 

134 

135class _GeoidBase(_HeightsBase): 

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

137 ''' 

138 _cropped = None 

139# _datum = _WGS84 # from _HeightsBase 

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

141 _endian = _tbd_ 

142 _geoid = _n_a_ 

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

144 _interp2d = None # interp2d interpolation 

145 _kind = 3 # order for interp2d, RectBivariateSpline 

146# _kmin = 2 # min number of knots 

147 _knots = 0 # nlat * nlon 

148 _mean = None # fixed in GeoidKarney 

149# _name = NN # _Named 

150 _nBytes = 0 # numpy size in bytes, float64 

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

152 _sizeB = 0 # geoid file size in bytes 

153 _smooth = 0 # used only for RectBivariateSpline 

154 _stdev = None # fixed in GeoidKarney 

155 _u2B = 0 # np.itemsize or undefined 

156 

157 _lat_d = _0_0 # increment, +tive 

158 _lat_lo = _0_0 # lower lat, south 

159 _lat_hi = _0_0 # upper lat, noth 

160 _lon_d = _0_0 # increment, +tive 

161 _lon_lo = _0_0 # left lon, west 

162 _lon_hi = _0_0 # right lon, east 

163 _lon_of = _0_0 # forward lon offset 

164 _lon_og = _0_0 # reverse lon offset 

165 

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

167 _yx_hits = None # cache hits, ala Karney 

168 

169 def __init__(self, hs, p): 

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

171 and several internal geoid attributes. 

172 

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

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

175 other geoid parameters (C{INTERNAL}). 

176 

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

178 invalid B{C{kind}}. 

179 

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

181 

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

183 C{-.RectBivariateSpline} issue. 

184 

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

186 C{-.RectBivariateSpline} warning as 

187 exception. 

188 ''' 

189 spi = self.scipy_interpolate 

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

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

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

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

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

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

196 

197 # both axes and bounding box 

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

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

200 

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

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

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

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

205 k = self.kind 

206 if k in _interp2d_ks: 

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

208 elif 1 <= k <= 5: 

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

210 s=self._smooth).ev 

211 else: 

212 raise GeoidError(kind=k) 

213 

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

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

216 self._knots = p.knots # grid knots 

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

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

219 # shrink the box by 1 unit on every side 

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

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

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

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

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

225 

226 def __call__(self, *llis, **wrap): 

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

228 

229 @arg llis: The location or locations (C{LatLon}, ... or 

230 C{LatLon}s). 

231 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}} 

232 locations (C{bool}). 

233 

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

235 a list or tuple of interpolated geoid heights 

236 (C{float}s). 

237 

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

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

240 geoid file is closed. 

241 

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

243 or longitude range. 

244 

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

246 C{-.RectBivariateSpline} issue. 

247 

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

249 C{-.RectBivariateSpline} warning as 

250 exception. 

251 ''' 

252 return self._called(llis, True, **wrap) 

253 

254 def __enter__(self): 

255 '''Open context. 

256 ''' 

257 return self 

258 

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

260 '''Close context. 

261 ''' 

262 self.close() 

263 # return None # XXX False 

264 

265 def __repr__(self): 

266 return self.toStr() 

267 

268 def __str__(self): 

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

270 

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

272 # handle __call__ 

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

274 try: 

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

276 for i, lli in enumerate(llis): 

277 hs.append(self._hGeoid(*_w(lli.lat, lli.lon))) 

278 return _as(hs) 

279 

280 except (GeoidError, RangeError) as x: 

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

282 t = _lli_ if _as is _ascalar else Fmt.SQUARE(llis=i) 

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

284 raise type(x)(t, lli, wrap=wrap, cause=x) 

285 except Exception as x: 

286 if scipy and self.scipy: 

287 raise _SciPyIssue(x) 

288 else: 

289 raise 

290 

291 @Property_RO 

292 def _center(self): 

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

294 ''' 

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

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

297 

298 def center(self, LatLon=None): 

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

300 

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

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

303 

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

305 lon, height)} otherwise a B{C{LatLon}} instance 

306 with the lat-, longitude and height of the grid 

307 center location. 

308 ''' 

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

310 

311 def close(self): 

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

313 ''' 

314 if not self.closed: 

315 self._egm.close() 

316 self._egm = None 

317 

318 @property_RO 

319 def closed(self): 

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

321 ''' 

322 return self._egm is None 

323 

324 @Property_RO 

325 def cropped(self): 

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

327 ''' 

328 return self._cropped 

329 

330 @Property_RO 

331 def dtype(self): 

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

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

334 ''' 

335 return self._hs_y_x.dtype 

336 

337 @Property_RO 

338 def endian(self): 

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

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

341 ''' 

342 return self._endian 

343 

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

345 # only used for .interpolate.interp2d, but 

346 # overwritten for .RectBivariateSpline, 

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

348 return self._interp2d(x, y) 

349 

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

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

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

353 if m != n: 

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

355 if d < 0: 

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

357 for i in range(1, m): 

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

359 if e < EPS: # non-increasing axis 

360 i = Fmt.SQUARE(name, i) 

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

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

363 

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

365 notOverloaded(self, lat, lon) 

366 

367 def _gyx2g2(self, y, x): 

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

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

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

371 

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

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

374 

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

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

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

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

379 

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

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

382 

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

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

385 

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

387 geoid's lat- or longitude range. 

388 

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

390 C{-.RectBivariateSpline} issue. 

391 

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

393 C{-.RectBivariateSpline} warning as 

394 exception. 

395 ''' 

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

397 

398 def _hGeoid(self, lat, lon): 

399 out = self.outside(lat, lon) 

400 if out: 

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

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

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

404 

405 @Property_RO 

406 def _highest(self): 

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

408 ''' 

409 return self._llh3minmax(True) 

410 

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

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

413 

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

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

416 

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

418 lon, height)} otherwise a B{C{LatLon}} instance 

419 with the lat-, longitude and height of the highest 

420 grid location. 

421 ''' 

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

423 

424 @Property_RO 

425 def hits(self): 

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

427 ''' 

428 return self._yx_hits 

429 

430 @Property_RO 

431 def kind(self): 

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

433 ''' 

434 return self._kind 

435 

436 @Property_RO 

437 def knots(self): 

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

439 ''' 

440 return self._knots 

441 

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

443 notOverloaded(self, lat, lon) 

444 

445 @property_RO 

446 def _LL3T(self): 

447 '''(INTERNAL) Get L{LatLon3Tuple} once. 

448 ''' 

449 T = _MODS.namedTuples.LatLon3Tuple 

450 _GeoidBase._LL3T = T # overwrite poperty_RO 

451 return T 

452 

453 def _llh3(self, lat, lon): 

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

455 

456 def _llh3LL(self, llh, LatLon): 

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

458 

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

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

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

462 # numpy.argmin.html#numpy.argmin> 

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

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

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

466 

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

468 # numpy.fromfile, like .frombuffer 

469 g.seek(offset, _SEEK_SET) 

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

471 

472 @Property_RO 

473 def _lowerleft(self): 

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

475 ''' 

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

477 

478 def lowerleft(self, LatLon=None): 

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

480 

481 @kwarg LatLon: Optional class to return the location 

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

483 

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

485 lon, height)} otherwise a B{C{LatLon}} instance 

486 with the lat-, longitude and height of the lower-left, 

487 SW grid corner. 

488 ''' 

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

490 

491 @Property_RO 

492 def _loweright(self): 

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

494 ''' 

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

496 

497 def loweright(self, LatLon=None): 

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

499 

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

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

502 

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

504 lon, height)} otherwise a B{C{LatLon}} instance 

505 with the lat-, longitude and height of the lower-right, 

506 SE grid corner. 

507 ''' 

508 

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

510 

511 lowerright = loweright # synonymous 

512 

513 @Property_RO 

514 def _lowest(self): 

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

516 ''' 

517 return self._llh3minmax(False) 

518 

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

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

521 

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

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

524 

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

526 lon, height)} otherwise a B{C{LatLon}} instance 

527 with the lat-, longitude and height of the lowest 

528 grid location. 

529 ''' 

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

531 

532 @Property_RO 

533 def mean(self): 

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

535 ''' 

536 if self._mean is None: # see GeoidKarney 

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

538 return self._mean 

539 

540 @property_RO 

541 def name(self): 

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

543 ''' 

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

545 

546 @Property_RO 

547 def nBytes(self): 

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

549 ''' 

550 return self._nBytes 

551 

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

553 # open the geoid file 

554 try: 

555 self._geoid = _os_path.basename(geoid) 

556 self._sizeB = _os_path.getsize(geoid) 

557 g = open(geoid, _rb_) 

558 except (IOError, OSError) as x: 

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

560 

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

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

563 self._kind = int(kind) 

564 if name: 

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

566 if smooth: 

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

568 

569 return g 

570 

571 def outside(self, lat, lon): 

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

573 lat-/longitude or crop range. 

574 

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

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

577 

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

579 empty C{str} if inside. 

580 ''' 

581 return (_S_ if lat < self._lat_lo else 

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

583 (_W_ if lon < self._lon_lo else 

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

585 

586 @Property_RO 

587 def pgm(self): 

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

589 ''' 

590 return self._pgm 

591 

592 @Property_RO 

593 def sizeB(self): 

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

595 ''' 

596 return self._sizeB 

597 

598 @Property_RO 

599 def smooth(self): 

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

601 ''' 

602 return self._smooth 

603 

604 @Property_RO 

605 def stdev(self): 

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

607 ''' 

608 if self._stdev is None: # see GeoidKarney 

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

610 return self._stdev 

611 

612 def _swne(self, crop): 

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

614 try: 

615 if len(crop) == 2: 

616 try: # sw, ne LatLons 

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

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

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

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

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

622 swne = crop 

623 if len(swne) == 4: 

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

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

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

627 return s, w, n, e 

628 except (IndexError, TypeError, ValueError): 

629 pass 

630 raise GeoidError(crop=crop) 

631 

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

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

634 

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

636 default). Trailing zero decimals are stripped 

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

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

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

640 

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

642 ''' 

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

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

645 (self.lowerleft, self.upperright, 

646 self.center, 

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

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

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

650 attrs( 'cropped', 'dtype', _endian_, 'hits', _knots_, 'nBytes', 

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

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

653 

654 @Property_RO 

655 def u2B(self): 

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

657 ''' 

658 return self._u2B 

659 

660 @Property_RO 

661 def _upperleft(self): 

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

663 ''' 

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

665 

666 def upperleft(self, LatLon=None): 

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

668 

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

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

671 

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

673 lon, height)} otherwise a B{C{LatLon}} instance 

674 with the lat-, longitude and height of the upper-left, 

675 NW grid corner. 

676 ''' 

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

678 

679 @Property_RO 

680 def _upperright(self): 

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

682 ''' 

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

684 

685 def upperright(self, LatLon=None): 

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

687 

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

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

690 

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

692 lon, height)} otherwise a B{C{LatLon}} instance 

693 with the lat-, longitude and height of the upper-right, 

694 NE grid corner. 

695 ''' 

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

697 

698 

699class GeoidError(HeightError): 

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

701 ''' 

702 pass 

703 

704 

705class GeoidG2012B(_GeoidBase): 

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

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

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

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

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

711 U{Guam and Northern Mariana Islands 

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

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

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

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

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

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

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

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

720 

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

722 C{g2012b*.bin} grid files. 

723 ''' 

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

725 kind=3, name=NN, smooth=0): 

726 '''New L{GeoidG2012B} interpolator. 

727 

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

729 @kwarg crop: Optional crop box, not supported (C{None}). 

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

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

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

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

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

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

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

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

738 @kwarg name: Optional geoid name (C{str}). 

739 @kwarg smooth: Smoothing factor for U{RectBivariateSpline 

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

741 scipy.interpolate.RectBivariateSpline.html>} 

742 only (C{int}). 

743 

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

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

746 

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

748 not installed. 

749 

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

751 

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

753 

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

755 warning as exception. 

756 

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

758 ''' 

759 if crop is not None: 

760 raise GeoidError(crop=crop, txt=_not_(_supported_)) 

761 

762 g = self._open(g2012b_bin, datum, kind, name, smooth) 

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

764 

765 try: 

766 p = _Gpars() 

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

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

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

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

771 # skip 4xf8, get 3xi4 

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

773 if ien == 1: # correct endian 

774 p.knots = p.nlat * p.nlon 

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

776 and 1 < p.nlon < n: 

777 self._endian = en_+'f4' 

778 break 

779 else: # couldn't validate endian 

780 raise GeoidError(_endian_) 

781 

782 # get the first 4xf8 

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

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

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

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

787 _GeoidBase.__init__(self, hs, p) 

788 

789 except Exception as x: 

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

791 finally: 

792 g.close() 

793 

794 def _g2ll2(self, lat, lon): 

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

796 return lat, lon 

797 

798 def _ll2g2(self, lat, lon): 

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

800 return lat, lon 

801 

802 if _FOR_DOCS: 

803 __call__ = _GeoidBase.__call__ 

804 height = _GeoidBase.height 

805 

806 

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

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

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

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

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

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

813 ''' 

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

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

816 

817 

818def _I(i): 

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

820 ''' 

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

822 

823 

824def _T(*cs): 

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

826 ''' 

827 return map1(_I, *cs) 

828 

829 

830class GeoidKarney(_GeoidBase): 

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

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

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

834 html/geoid.html#geoidinst>} datasets using bilinear or U{cubic 

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

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

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

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

839 

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

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

842 datasets. 

843 ''' 

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

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

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

847 _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 

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

849 _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 

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

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

852 _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 

853 _T(0, 0, 0, 0, 0, 0, 93, -93, -93, 93, 0, 0), 

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

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

856 

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

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

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

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

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

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

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

864 _T(0, 0, -60, 60, 60, -60, 60, -60, -60, 60, 0, 0), 

865 _T(-60, 60, 0, 60, -60, 0, 0, 60, -60, 0, -60, 60), 

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

867 

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

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

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

871 _T(0, 0, 93, -93, -93, 93, 0, 0, 0, 0, 0, 0), 

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

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

874 _T(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 

875 _T(0, 0, -93, 93, 93, -93, 0, 0, 0, 0, 0, 0), 

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

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

878 

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

880 _T(1, 0), 

881 _T(0, 1), 

882 _T(1, 1)) 

883 

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

885 _T( 1, -1), 

886 _T(-1, 0), 

887 _T( 0, 0), 

888 _T( 1, 0), 

889 _T( 2, 0), 

890 _T(-1, 1), 

891 _T( 0, 1), 

892 _T( 1, 1), 

893 _T( 2, 1), 

894 _T( 0, 2), 

895 _T( 1, 2)) 

896 

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

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

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

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

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

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

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

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

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

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

907 _nBytes = None # not applicable 

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

909 _smooth = None # not applicable 

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

911 _u2B = _calcsize(_endian) # pixelsize_ in bytes 

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

913 _Ru2B = 0 # row of pixelsize_s in bytes 

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

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

916 _yx_hits = 0 # cache hits 

917 

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

919 kind=3, name=NN, smooth=None): 

920 '''New L{GeoidKarney} interpolator. 

921 

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

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

924 note below. 

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

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

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

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

929 C{LatLon} instances. 

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

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

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

933 for C{cubic}. 

934 @kwarg name: Optional geoid name (C{str}). 

935 @kwarg smooth: Smoothing factor, unsupported (C{None}). 

936 

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

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

939 

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

941 

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

943 

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

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

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

947 ''' 

948 if smooth is not None: 

949 raise GeoidError(smooth=smooth, txt=_not_(_supported_)) 

950 

951 if kind in (2,): 

952 self._evH = self._ev2H 

953 elif kind not in (3,): 

954 raise GeoidError(kind=kind) 

955 

956 self._egm = g = self._open(egm_pgm, datum, kind, name, smooth) 

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

958 

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

960 self._Ru2B = _calcsize(self._Rendian) 

961 

962 self._knots = p.knots # grid knots 

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

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

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

966 self._lat_lo, \ 

967 self._lon_lo, \ 

968 self._lat_hi, \ 

969 self._lon_hi = self._swne(crop if crop else p.crop4) 

970 self._cropped = True if crop else False 

971 

972 def __call__(self, *llis, **wrap): 

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

974 

975 @arg llis: The location or locations (C{LatLon}, ... or 

976 C{LatLon}s). 

977 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}} 

978 locations (C{bool}). 

979 

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

981 a list or tuple of interpolated geoid heights 

982 (C{float}s). 

983 

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

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

986 geoid file is closed. 

987 

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

989 or longitude range. 

990 ''' 

991 return self._called(llis, False, **wrap) 

992 

993 def _c0c3v(self, y, x): 

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

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

996 p = self._pgm 

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

998 # read 4x4 ushorts, drop the 4 corners 

999 g = self._egm 

1000 e = self._4endian 

1001 n = self._4u2B 

1002 R = self._Ru2B 

1003 

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

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

1006 b += R 

1007 g.seek(b, _SEEK_SET) 

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

1009 b += R 

1010 g.seek(b, _SEEK_SET) 

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

1012 b += R 

1013 g.seek(b, _SEEK_SET) 

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

1015 j = 1 

1016 

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

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

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

1020 

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

1022 

1023 @Property_RO 

1024 def dtype(self): 

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

1026 ''' 

1027 return 'ushort' 

1028 

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

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

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

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

1033 fy -= y 

1034 fx -= x 

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

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

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

1038 return H.fsum() 

1039 

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

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

1042 if self._yxH == yx: 

1043 t = self._yxHt 

1044 self._yx_hits += 1 

1045 else: 

1046 y, x = self._yxH = yx 

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

1048 v = _1_0, -fx, fx 

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

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

1051 return H 

1052 

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

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

1055 if self._yxH == yx: 

1056 t = self._yxHt 

1057 self._yx_hits += 1 

1058 else: 

1059 self._yxH = yx 

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

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

1062 self._yxHt = t = tuple(t) 

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

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

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

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

1067 v = _1_0, fx, fy 

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

1069 H *= fy 

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

1071 H *= fy 

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

1073 return H 

1074 

1075 _evH = _ev3H # overriden for kind == 2 

1076 

1077 def _g2ll2(self, lat, lon): 

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

1079 while lon > _180_0: 

1080 lon -= _360_0 

1081 return lat, lon 

1082 

1083 def _g2yx2(self, lat, lon): 

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

1085 p = self._pgm 

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

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

1088 

1089 def _gyx2g2(self, y, x): 

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

1091 p = self._pgm 

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

1093 

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

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

1096 

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

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

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

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

1101 

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

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

1104 

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

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

1107 geoid file is closed. 

1108 

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

1110 geoid's lat- or longitude range. 

1111 ''' 

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

1113 

1114 @Property_RO 

1115 def _highest_ltd(self): 

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

1117 ''' 

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

1119 

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

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

1122 

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

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

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

1126 

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

1128 lon, height)} otherwise a B{C{LatLon}} instance 

1129 with the lat-, longitude and height of the highest 

1130 grid location. 

1131 ''' 

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

1133 return self._llh3LL(llh, LatLon) 

1134 

1135 def _lat2y2(self, lat2): 

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

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

1138 for lat in lat2: 

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

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

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

1142 

1143 def _ll2g2(self, lat, lon): 

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

1145 while lon < 0: 

1146 lon += _360_0 

1147 return lat, lon 

1148 

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

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

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

1152 y = x = 0 

1153 h = self._raw(y, x) 

1154 if highest: 

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

1156 m = max(r) 

1157 if m > h: 

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

1159 else: # lowest 

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

1161 m = min(r) 

1162 if m < h: 

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

1164 h *= self._pgm.Scale 

1165 h += self._pgm.Offset 

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

1167 

1168 @Property_RO 

1169 def _lowest_ltd(self): 

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

1171 ''' 

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

1173 

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

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

1176 

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

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

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

1180 

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

1182 lon, height)} otherwise a B{C{LatLon}} instance 

1183 with the lat-, longitude and height of the lowest 

1184 grid location. 

1185 ''' 

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

1187 return self._llh3LL(llh, LatLon) 

1188 

1189 def _raw(self, y, x): 

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

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

1192 p = self._pgm 

1193 if x < 0: 

1194 x += p.nlon 

1195 elif x >= p.nlon: 

1196 x -= p.nlon 

1197 h = p.nlon // 2 

1198 if y < 0: 

1199 y = -y 

1200 elif y >= p.nlat: 

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

1202 else: 

1203 h = 0 

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

1205 self._seek(y, x) 

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

1207 return h[0] 

1208 

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

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

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

1212 

1213 def _raw2(self, *lat2): 

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

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

1216 p = self._pgm 

1217 g = self._egm 

1218 e = self._Rendian 

1219 n = self._Ru2B 

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

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

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

1223 for y in range(s, t): 

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

1225 

1226 def _seek(self, y, x): 

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

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

1229 if g: 

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

1231 g.seek(b, _SEEK_SET) 

1232 return b # position 

1233 raise GeoidError('closed file: %r' % (p.egm,)) # IOError 

1234 

1235 

1236class GeoidPGM(_GeoidBase): 

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

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

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

1240 html/geoid.html#geoidinst>} datasets but based on C{SciPy} 

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

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

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

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

1245 

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

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

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

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

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

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

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

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

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

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

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

1257 ''' 

1258 _endian = '>u2' 

1259 

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

1261 kind=3, name=NN, smooth=0): 

1262 '''New L{GeoidPGM} interpolator. 

1263 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1278 @kwarg name: Optional geoid name (C{str}). 

1279 @kwarg smooth: Smoothing factor for U{RectBivariateSpline 

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

1281 scipy.interpolate.RectBivariateSpline.html>} 

1282 only (C{int}). 

1283 

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

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

1286 

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

1288 

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

1290 

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

1292 

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

1294 warning as exception. 

1295 

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

1297 

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

1299 html/geoid.html#geoidinst>} file sizes are based on a 2-byte 

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

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

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

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

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

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

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

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

1308 

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

1310 ''' 

1311 np = self.numpy 

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

1313 

1314 g = self._open(egm_pgm, datum, kind, name, smooth) 

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

1316 if crop: 

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

1318 self._g2ll2 = self._g2ll2_cropped 

1319 self._ll2g2 = self._ll2g2_cropped 

1320 if map2(int, np.__version__.split(_DOT_)[:2]) < (1, 9): 

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

1322 self._cropped = True 

1323 else: 

1324 self._cropped = False 

1325 

1326 try: 

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

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

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

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

1331 if p.Offset: # offset 

1332 hs = p.Offset + hs 

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

1334 hs = np.flipud(hs) 

1335 _GeoidBase.__init__(self, hs, p) 

1336 except Exception as x: 

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

1338 finally: 

1339 g.close() 

1340 

1341 def _g2ll2(self, lat, lon): 

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

1343 while lon > _180_0: 

1344 lon -= _360_0 

1345 return lat, lon 

1346 

1347 def _g2ll2_cropped(self, lat, lon): 

1348 # convert cropped grid (lat, lon) to earth (lat, lon) 

1349 return lat, lon - self._lon_og 

1350 

1351 def _ll2g2(self, lat, lon): 

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

1353 while lon < 0: 

1354 lon += _360_0 

1355 return lat, lon 

1356 

1357 def _ll2g2_cropped(self, lat, lon): 

1358 # convert earth (lat, lon) to cropped grid (lat, lon) 

1359 return lat, lon + self._lon_of 

1360 

1361 if _FOR_DOCS: 

1362 __call__ = _GeoidBase.__call__ 

1363 height = _GeoidBase.height 

1364 

1365 

1366class _Gpars(_Named): 

1367 '''(INTERNAL) Basic geoid parameters. 

1368 ''' 

1369 # interpolator parameters 

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

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

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

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

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

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

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

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

1378 

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

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

1381 

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

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

1384 

1385 def __repr__(self): 

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

1387 a in dir(self.__class__) 

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

1389 return _COLONSPACE_(self, t) 

1390 

1391 def __str__(self): 

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

1393 

1394 

1395class _PGM(_Gpars): 

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

1397 

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

1399 # Description WGS84 EGM96, 5-minute grid 

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

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

1402 # MaxBilinearError 0.140 

1403 # RMSBilinearError 0.005 

1404 # MaxCubicError 0.003 

1405 # RMSCubicError 0.001 

1406 # Offset -108 

1407 # Scale 0.003 

1408 # Origin 90N 0E 

1409 # AREA_OR_POINT Point 

1410 # Vertical_Datum WGS84 

1411 <width> <height> 

1412 <pixel> 

1413 ... 

1414 ''' 

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

1416 egm = None 

1417 glon = 180 # reverse offset, uncropped 

1418# pgm = NN # name 

1419 sizeB = 0 

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

1421 

1422 @staticmethod 

1423 def _llstr2floats(latlon): 

1424 # llstr to (lat, lon) floats 

1425 lat, lon = latlon.split() 

1426 return parseDMS2(lat, lon) 

1427 

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

1429 AREA_OR_POINT = str 

1430 DateTime = str 

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

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

1433 MaxBilinearError = float 

1434 MaxCubicError = float 

1435 Offset = float 

1436 Origin = _llstr2floats 

1437 Pixel = 0 

1438 RMSBilinearError = float 

1439 RMSCubicError = float 

1440 Scale = float 

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

1442 Vertical_Datum = str 

1443 

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

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

1446 ''' 

1447 self.name = pgm # geoid file name 

1448 if itemsize: 

1449 self._u2B = itemsize 

1450 if sizeB: 

1451 self.sizeB = sizeB 

1452 

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

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

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

1456 

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

1458 try: # ignore empty ones or comments 

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

1460 if t.startswith(_bHASH_): 

1461 t = t.lstrip(_bHASH_).lstrip() 

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

1463 f = getattr(_PGM, a, None) 

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

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

1466 elif t: 

1467 break 

1468 except (TypeError, ValueError): 

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

1470 else: # should never get here 

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

1472 

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

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

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

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

1477 raise ValueError 

1478 except (TypeError, ValueError): 

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

1480 

1481 try: # must be 16 bit pixel height 

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

1483 self.Pixel = int(t) 

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

1485 raise ValueError 

1486 except (TypeError, ValueError): 

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

1488 

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

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

1491 setattr(self, a, None) 

1492 

1493 if self.Origin is None: 

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

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

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

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

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

1499 

1500 self.skip = g.tell() 

1501 self.knots = nlat * nlon 

1502 

1503 self.nlat, self.nlon = nlat, nlon 

1504 self.slat, self.wlon = self.Origin 

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

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

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

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

1509 

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

1511 n = float(self.slat) 

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

1513 w = self.wlon - self.glon 

1514 e = w + self.dlon * nlon 

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

1516 

1517 n = self.sizeB - self.skip 

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

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

1520 

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

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

1523 ''' 

1524 # flon offset for both west and east 

1525 f = 360 if west < 0 else 0 

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

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

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

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

1530 s += 1 # s > n 

1531 e += 1 # e > w 

1532 

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

1534 # handle special cases 

1535 if (s - n) > hi: 

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

1537 if (e - w) > wi: 

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

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

1540 return g # use entire geoid as-is 

1541 

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

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

1544 

1545 if e > wi > w: # wrap around 

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

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

1548 elif e > w: 

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

1550 else: 

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

1552 

1553 # convert to bytes 

1554 r *= self.u2B 

1555 p *= self.u2B 

1556 q = wi * self.u2B # stride 

1557 # number of rows and cols to skip from 

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

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

1560 # sanity check 

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

1562 or z > self.sizeB: 

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

1564 

1565 # can't use _BytesIO since numpy 

1566 # needs .fileno attr in .fromfile 

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

1568 # reading (s - n) rows, forward 

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

1570 g.seek(z, _SEEK_SET) 

1571 # Python 2 tmpfile.write returns None 

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

1573 if p: # wrap around to start of row 

1574 g.seek(-q, _SEEK_CUR) 

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

1576 # Python 2 tmpfile.write returns None 

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

1578 z += q 

1579 c.flush() 

1580 g.close() 

1581 

1582 s -= n # nlat 

1583 e -= w # nlon 

1584 k = s * e # knots 

1585 z = k * self.u2B 

1586 if t != z: 

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

1588 

1589 # update the _Gpars accordingly, note attributes 

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

1591 self.slat += n * self.dlat 

1592 self.wlon += w * self.dlon 

1593 self.nlat = s 

1594 self.nlon = e 

1595 self.flon = self.glon = f 

1596 

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

1598 self.knots = k 

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

1600 

1601 c.seek(0, _SEEK_SET) 

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

1603 return c 

1604 

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

1606 t = fmt % args 

1607 e = self.pgm or NN 

1608 if e: 

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

1610 return PGMError(t) 

1611 

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

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

1614 # with .dlat decreasing from 90N .slat 

1615 lat -= self.slat 

1616 lon += flon - self.wlon 

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

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

1619 

1620 def _tmpfile(self): 

1621 # create a tmpfile to hold the cropped geoid grid 

1622 try: 

1623 from tempfile import NamedTemporaryFile as tmpfile 

1624 except ImportError: # Python 2.7.16- 

1625 from os import tmpfile 

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

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

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

1629 return f 

1630 

1631 @Property_RO 

1632 def pgm(self): 

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

1634 ''' 

1635 return self.name 

1636 

1637 

1638class PGMError(GeoidError): 

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

1640 ''' 

1641 pass 

1642 

1643 

1644def egmGeoidHeights(GeoidHeights_dat): 

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

1646 html/geoid.html#geoidinst>} height tests from U{GeoidHeights.dat 

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

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

1649 geoid.html#testgeoid>}. 

1650 

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

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

1653 

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

1655 egm84, egm96, egm2008)}. 

1656 

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

1658 

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

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

1661 C{test/testGeoids.py}. 

1662 ''' 

1663 dat = GeoidHeights_dat 

1664 if isinstance(dat, bytes): 

1665 dat = _BytesIO(dat) 

1666 

1667 try: 

1668 dat.seek(0, _SEEK_SET) # reset 

1669 except AttributeError as x: 

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

1671 

1672 for t in dat.readlines(): 

1673 t = t.strip() 

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

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

1676 while lon > _180_0: # EasternLon to earth lon 

1677 lon -= _360_0 

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

1679 

1680 

1681__all__ += _ALL_DOCS(_GeoidBase) 

1682 

1683if __name__ == '__main__': 

1684 

1685 from pygeodesy.lazily import printf, _sys 

1686 

1687 _crop = () 

1688 _GeoidEGM = GeoidKarney 

1689 _kind = 3 

1690 

1691 geoids = _sys.argv[1:] 

1692 while geoids: 

1693 geoid = geoids.pop(0) 

1694 

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

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

1697 

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

1699 _GeoidEGM = GeoidKarney 

1700 

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

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

1703 

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

1705 _GeoidEGM = GeoidPGM 

1706 

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

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

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

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

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

1712 # The height of the EGM96 geoid at Timbuktu 

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

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

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

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

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

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

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

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

1721 ll = parseDMS2('16:46:33N', '3:00:34W', sep=':') 

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

1723 try: 

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

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

1726 except (GeoidError, RangeError) as x: 

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

1728 

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

1730 g = GeoidG2012B(geoid, kind=_kind) 

1731 printf(g.toStr()) 

1732 

1733 else: 

1734 raise GeoidError(grid=repr(geoid)) 

1735 

1736_I = int # PYCHOK unused _I 

1737del _intCs # trash ints cache 

1738 

1739# **) MIT License 

1740# 

1741# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

1742# 

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

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

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

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

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

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

1749# 

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

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

1752# 

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

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

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

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

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

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

1759# OTHER DEALINGS IN THE SOFTWARE. 

1760 

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

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

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

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

1765 

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

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

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

1769 

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

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

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

1773 

1774 

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

1776# 

1777# GeoidKarney('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, 147.367, 85.839), lowest(4.7, 78.767, -106.911) 

1778# 

1779# _PGM('../testGeoids/egm2008-1.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-31 06:54:00', Description='WGS84 EGM2008, 1-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.025, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.001, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008', Vertical_Datum='WGS84' 

1780# 

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

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

1783# 

1784# GeoidKarney('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, 148.75, 81.33), lowest(4.75, 79.25, -107.34) 

1785# 

1786# _PGM('../testGeoids/egm84-15.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 15-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.413, MaxCubicError=0.02, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84' 

1787# 

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

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

1790# 

1791# GeoidKarney('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, 0.0, 17.163), highest(-8.167, 147.25, 85.422), lowest(4.667, 78.833, -107.043) 

1792# 

1793# _PGM('../testGeoids/egm96-5.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:03', Description='WGS84 EGM96, 5-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.14, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.005, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html', Vertical_Datum='WGS84' 

1794# 

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

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

1797 

1798 

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

1800# 

1801# GeoidKarney('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, 147.367, 85.839), lowest(4.7, 78.767, -106.911) 

1802# 

1803# _PGM('../testGeoids/egm2008-1.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-31 06:54:00', Description='WGS84 EGM2008, 1-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.025, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.001, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008', Vertical_Datum='WGS84' 

1804# 

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

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

1807# 

1808# GeoidKarney('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, 148.75, 81.33), lowest(4.75, 79.25, -107.34) 

1809# 

1810# _PGM('../testGeoids/egm84-15.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 15-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.413, MaxCubicError=0.02, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84' 

1811# 

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

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

1814# 

1815# GeoidKarney('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, 0.0, 17.163), highest(-8.167, 147.25, 85.422), lowest(4.667, 78.833, -107.043) 

1816# 

1817# _PGM('../testGeoids/egm96-5.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:03', Description='WGS84 EGM96, 5-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.14, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.005, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html', Vertical_Datum='WGS84' 

1818# 

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

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

1821 

1822 

1823# % python2 -m pygeodesy.geoids -PGM ../testGeoids/egm*.pgm 

1824# 

1825# GeoidPGM('egm2008-1.pgm'): lowerleft(-90.0, -180.0, -30.15), upperright(90.0, 180.0, 14.898), center(0.0, 0.0, 17.226), highest(-8.4, -32.633, 85.839), lowest(4.683, -101.25, -106.911) 

1826# 

1827# _PGM('../testGeoids/egm2008-1.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-31 06:54:00', Description='WGS84 EGM2008, 1-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.025, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.001, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008', Vertical_Datum='WGS84' 

1828# 

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

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

1831# 

1832# GeoidPGM('egm84-15.pgm'): lowerleft(-90.0, -180.0, -29.712), upperright(90.0, 180.0, 13.098), center(0.0, 0.0, 18.33), highest(-4.5, -31.25, 81.33), lowest(4.75, -100.75, -107.34) 

1833# 

1834# _PGM('../testGeoids/egm84-15.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:02', Description='WGS84 EGM84, 15-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.413, MaxCubicError=0.02, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.018, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html', Vertical_Datum='WGS84' 

1835# 

1836# Timbuktu GeoidPGM('egm84-15.pgm').height(16.775833, -3.009444): 31.2979 vs 31.2979 

1837# Timbuktu GeoidPGM('egm84-15.pgm').height(16.776, -3.009): 31.2975 vs 31.2979 

1838# 

1839# GeoidPGM('egm96-5.pgm'): lowerleft(-90.0, -180.0, -29.535), upperright(90.0, 180.0, 13.605), center(0.0, -0.0, 17.179), highest(-8.167, -32.75, 85.422), lowest(4.667, -101.167, -107.043) 

1840# 

1841# _PGM('../testGeoids/egm96-5.pgm'): AREA_OR_POINT='Point', DateTime='2009-08-29 18:45:03', Description='WGS84 EGM96, 5-minute grid', Geoid='file in PGM format for the GeographicLib::Geoid class', MaxBilinearError=0.14, MaxCubicError=0.003, Offset=-108.0, Origin=LatLon2Tuple(lat=90.0, lon=0.0), Pixel=65535, RMSBilinearError=0.005, RMSCubicError=0.001, Scale=0.003, URL='http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html', Vertical_Datum='WGS84' 

1842# 

1843# Timbuktu GeoidPGM('egm96-5.pgm').height(16.775833, -3.009444): 28.7065 vs 28.7067 

1844# Timbuktu GeoidPGM('egm96-5.pgm').height(16.776, -3.009): 28.7064 vs 28.7067