Coverage for pygeodesy/heights.py: 95%

1 

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

3 

4u'''Height interpolations at C{LatLon} points from known C{knots}. 

5 

6Classes L{HeightCubic}, L{HeightIDWcosineLaw}, L{HeightIDWdistanceTo}, 

7L{HeightIDWequirectangular}, L{HeightIDWeuclidean}, L{HeightIDWflatLocal}, 

8L{HeightIDWflatPolar}, L{HeightIDWhaversine}, L{HeightIDWhubeny}, 

9L{HeightIDWkarney}, L{HeightIDWthomas}, L{HeightIDWvincentys}, L{HeightLinear}, 

10L{HeightLSQBiSpline} and L{HeightSmoothBiSpline} to interpolate the height of 

11C{LatLon} locations or separate lat-/longitudes from a set of C{LatLon} points 

12with I{known heights}. 

13 

14Typical usage 

15============= 

16 

171. Get or create a set of C{LatLon} points with I{known heights}, called 

18C{knots}. The C{knots} do not need to be ordered in any particular way. 

19 

20C{>>> ...} 

21 

222. Select one of the C{Height} classes for height interpolation 

23 

24C{>>> from pygeodesy import HeightCubic as HeightXyz # or an other Height... class} 

25 

263. Instantiate a height interpolator with the C{knots} and use keyword 

27arguments to select different interpolation options 

28 

29C{>>> hinterpolator = HeightXyz(knots, **options)} 

30 

314. Get the interpolated height of C{LatLon} location(s) with 

32 

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

34 

35C{>>> h = hinterpolator(ll)} 

36 

37or 

38 

39C{>>> h0, h1, h2, ... = hinterpolator(ll0, ll1, ll2, ...)} 

40 

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

42 

43C{>>> hs = hinterpolator(lls)} 

44 

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

46 

47C{>>> h = hinterpolator.height(lat, lon)} 

48 

49or as 2 lists, 2 tuples, etc. 

50 

51C{>>> hs = hinterpolator.height(lats, lons)} 

52 

53or for several positionals use the C{height_} method 

54 

55C{>>> h1, h2, ... = hinterpolator.height_(lat1, lon1, lat2, lon2, ...)} 

56 

57@note: Classes L{HeightCubic} and L{HeightLinear} require package U{numpy 

58 <https://PyPI.org/project/numpy>}, classes L{HeightLSQBiSpline} and 

59 L{HeightSmoothBiSpline} require package U{scipy<https://SciPy.org>}. 

60 Classes L{HeightIDWkarney} and L{HeightIDWdistanceTo} -if used with 

61 L{ellipsoidalKarney.LatLon} points- require I{Karney}'s U{geographiclib 

62 <https://PyPI.org/project/geographiclib>} package to be installed. 

63 

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

65 by C{scipy} can be thrown as L{SciPyWarning} exceptions, provided 

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

67 

68@see: U{SciPy<https://docs.SciPy.org/doc/scipy/reference/interpolate.html>} 

69 Interpolation. 

70''' 

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

72from __future__ import division as _; del _ # noqa: E702 ; 

73 

74from pygeodesy.basics import isscalar, len2, map1, min2, _xnumpy, _xscipy 

75from pygeodesy.constants import EPS, PI, PI_2, PI2, _0_0, _90_0, _180_0 

76from pygeodesy.datums import _ellipsoidal_datum, _WGS84 

77from pygeodesy.errors import _AssertionError, LenError, PointsError, \ 

78 _SciPyIssue, _xattr, _xkwds, _xkwds_get, _xkwds_item2 

79# from pygeodesy.fmath import fidw # _MODS 

80# from pygeodesy import formy as _formy # _MODS.into 

81# from pygeodesy.internals import _version2 # _MODS 

82from pygeodesy.interns import NN, _COMMASPACE_, _insufficient_, _NOTEQUAL_, \ 

83 _PLUS_, _scipy_, _SPACE_, _STAR_ 

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

85from pygeodesy.named import _name2__, _Named 

86from pygeodesy.points import _distanceTo, LatLon_, Fmt, radians, _Wrap 

87from pygeodesy.props import Property_RO, property_RO, property_ROver 

88# from pygeodesy.streprs import Fmt # from .points 

89from pygeodesy.units import _isDegrees, Float_, Int_ 

90# from pygeodesy.utily import _Wrap # from .points 

91 

92# from math import radians # from .points 

93 

94__all__ = _ALL_LAZY.heights 

95__version__ = '26.02.02' 

96 

97_error_ = 'error' 

98_formy = _MODS.into(formy=__name__) 

99_linear_ = 'linear' 

100_llis_ = 'llis' 

101 

102 

103class HeightError(PointsError): 

104 '''Height interpolator C{Height...} or interpolation issue. 

105 ''' 

106 pass 

107 

108 

109def _alist(ais): 

110 # return list of floats, not numpy.float64s 

111 return list(map(float, ais)) 

112 

113 

114def _ascalar(ais): # in .geoids 

115 # return single float, not numpy.float64 

116 ais = list(ais) # np.array, etc. to list 

117 if len(ais) != 1: 

118 n = Fmt.PAREN(len=repr(ais)) 

119 t = _SPACE_(len(ais), _NOTEQUAL_, 1) 

120 raise _AssertionError(n, txt=t) 

121 return float(ais[0]) # remove np.<type> 

122 

123 

124def _atuple(ais): 

125 # return tuple of floats, not numpy.float64s 

126 return tuple(map(float, ais)) 

127 

128 

129def _as_llis2(llis, m=1, Error=HeightError): # in .geoids 

130 # determine return type and convert lli C{LatLon}s to list 

131 if not isinstance(llis, tuple): # llis are *args 

132 n = Fmt.PAREN(type_=_STAR_(NN, _llis_)) 

133 raise _AssertionError(n, txt=repr(llis)) 

134 

135 n = len(llis) 

136 if n == 1: # convert single lli to 1-item list 

137 llis = llis[0] 

138 try: 

139 n, llis = len2(llis) 

140 _as = _alist # return list of interpolated heights 

141 except TypeError: # single lli 

142 n, llis = 1, [llis] 

143 _as = _ascalar # return single interpolated heights 

144 else: # of 0, 2 or more llis 

145 _as = _atuple # return tuple of interpolated heights 

146 

147 if n < m: 

148 raise _InsufficientError(m, Error=Error, llis=n) 

149 return _as, llis 

150 

151 

152def _InsufficientError(need, Error=HeightError, **name_value): # PYCHOK no cover 

153 # create an insufficient Error instance 

154 t = _COMMASPACE_(_insufficient_, str(need) + _PLUS_) 

155 return Error(txt=t, **name_value) 

156 

157 

158def _orderedup(ts, lo=EPS, hi=PI2-EPS): 

159 # clip, order and remove duplicates 

160 return sorted(set(max(lo, min(hi, t)) for t in ts)) # list 

161 

162 

163def _xyhs(wrap=False, _lat=_90_0, _lon=_180_0, height=True, **name_lls): 

164 # map (lat, lon, h) to (x, y, h) in radians, offset 

165 # x as 0 <= lon <= PI2 and y as 0 <= lat <= PI 

166 name, lls = _xkwds_item2(name_lls) 

167 _w, _r = _Wrap._latlonop(wrap), radians 

168 try: 

169 for i, ll in enumerate(lls): 

170 y, x = _w(ll.lat, ll.lon) 

171 h = ll.height if height else 0 

172 yield (max(_0_0, _r(x + _lon)), 

173 max(_0_0, _r(y + _lat)), h) 

174 except Exception as x: 

175 i = Fmt.INDEX(name, i) 

176 raise HeightError(i, ll, cause=x) 

177 

178 

179class _HeightNamed(_Named): # in .geoids 

180 '''(INTERNAL) Interpolator base class. 

181 ''' 

182 _datum = _WGS84 # default 

183 _Error = HeightError 

184 _kmin = 2 # min number of knots 

185 

186 _LLiC = LatLon_ # ._height class 

187 _np_sp = None # (numpy, scipy) 

188 _wrap = None # wrap knots and llis 

189 

190 def __call__(self, *llis, **wrap): # PYCHOK no cover 

191 '''Interpolate the height for one or several locations. I{Must be overloaded}. 

192 

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

194 @kwarg wrap: If C{B{wrap}=True} to wrap or I{normalize} all B{C{llis}} 

195 locations (C{bool}), overriding the B{C{knots}}' setting. 

196 

197 @return: A single interpolated height (C{float}) or a list or tuple of 

198 interpolated heights (each C{float}). 

199 

200 @raise HeightError: Insufficient number of B{C{llis}} or an invalid B{C{lli}}. 

201 

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

203 

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

205 ''' 

206 self._notOverloaded(callername='__call__', *llis, **wrap) 

207 

208 def _as_lls(self, lats, lons): # in .geoids 

209 LLiC, d = self._LLiC, self.datum 

210 if _isDegrees(lats) and _isDegrees(lons): 

211 llis = LLiC(lats, lons, datum=d) 

212 else: 

213 n, lats = len2(lats) 

214 m, lons = len2(lons) 

215 if n != m: # format a LenError, but raise self._Error 

216 e = LenError(type(self), lats=n, lons=m, txt=None) 

217 raise self._Error(str(e)) 

218 llis = [LLiC(*t, datum=d) for t in zip(lats, lons)] 

219 return llis 

220 

221 @property_RO 

222 def datum(self): 

223 '''Get the C{datum} setting or the default (L{Datum}). 

224 ''' 

225 return self._datum 

226 

227 def height(self, lats, lons, **wrap): # PYCHOK no cover 

228 '''I{Must be overloaded}.''' 

229 self._notOverloaded(lats, lons, **wrap) 

230 

231 def height_(self, *latlons, **wrap): 

232 '''Interpolate the height for each M{(latlons[i], latlons[i+1]) pair 

233 for i in range(0, len(latlons), B{2})}. 

234 

235 @arg latlons: Alternating lat-/longitude pairs (each C{degrees}), 

236 all positional. 

237 

238 @see: Method C{height} for further details. 

239 

240 @return: A tuple of interpolated heights (each C{float}). 

241 ''' 

242 lls = self._as_lls(latlons[0::2], latlons[1::2]) 

243 return tuple(self(lls, **wrap)) 

244 

245 @property_RO 

246 def kmin(self): 

247 '''Get the minimum number of knots (C{int}). 

248 ''' 

249 return self._kmin 

250 

251 @property_RO 

252 def wrap(self): 

253 '''Get the C{wrap} setting (C{bool}) or C{None}. 

254 ''' 

255 return self._wrap 

256 

257 

258class _HeightBase(_HeightNamed): # in .geoids 

259 '''(INTERNAL) Interpolator base class. 

260 ''' 

261 _k2interp2d = {-1: _linear_, # in .geoids._GeoidBase.__init__ 

262 -2: _linear_, # for backward compatibility 

263 -3: 'cubic', 

264 -5: 'quintic'} 

265 

266 def _as_xyllis4(self, llis, **wrap): 

267 # convert lli C{LatLon}s to tuples or C{NumPy} arrays of 

268 # C{SciPy} sphericals and determine the return type 

269 atype = self.numpy.array 

270 kwds = _xkwds(wrap, wrap=self._wrap, height=False) 

271 _as, llis = _as_llis2(llis) 

272 xis, yis, _ = zip(*_xyhs(llis=llis, **kwds)) # PYCHOK yield 

273 return _as, atype(xis), atype(yis), llis 

274 

275 def _ev(self, *args): # PYCHOK no cover 

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

277 self._notOverloaded(*args) 

278 

279 def _evalls(self, llis, **wrap): # XXX single arg, not *args 

280 _as, xis, yis, _ = self._as_xyllis4(llis, **wrap) 

281 try: # SciPy .ev signature: y first, then x! 

282 return _as(self._ev(yis, xis)) 

283 except Exception as x: 

284 raise _SciPyIssue(x, self._ev_name) 

285 

286 def _ev2d(self, x, y): # PYCHOK no cover 

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

288 self._notOverloaded(x, y) 

289 

290 @property_RO 

291 def _ev_name(self): 

292 '''(INTERNAL) Get the name of the C{.ev} method. 

293 ''' 

294 _ev = str(self._ev) 

295 if _scipy_ not in _ev: 

296 _ev = str(self._ev2d) 

297 # '<scipy.interpolate._interpolate.interp2d object at ...> 

298 # '<function _HeightBase._interp2d.<locals>._bisplev at ...> 

299 # '<bound method BivariateSpline.ev of ... object at ...> 

300 _ev = _ev[1:].split(None, 4) 

301 return Fmt.PAREN(_ev['sfb'.index(_ev[0][0])]) 

302 

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

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

305 

306 @arg lats: Latitude or latitudes (each C{degrees}). 

307 @arg lons: Longitude or longitudes (each C{degrees}). 

308 @kwarg wrap: Kewyord argument C{B{wrap}=False} (C{bool}). Use C{True} to 

309 wrap or I{normalize} all B{C{lats}} and B{C{lons}} locationts, 

310 overriding the B{C{knots}}' setting. 

311 

312 @return: A single interpolated height (C{float}) or a list of interpolated 

313 heights (each C{float}). 

314 

315 @raise HeightError: Insufficient or unequal number of B{C{lats}} and B{C{lons}}. 

316 

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

318 

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

320 ''' 

321 lls = self._as_lls(lats, lons) # dup of _HeightIDW.height 

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

323 

324 def _interp2d(self, xs, ys, hs, kind=-3): 

325 '''Create a C{scipy.interpolate.interp2d} or C{-.bisplrep/-ev} 

326 interpolator before, respectively since C{SciPy} version 1.14. 

327 ''' 

328 try: 

329 spi = self.scipy_interpolate 

330 if self._scipy_version() < (1, 14) and kind in self._k2interp2d: 

331 # SciPy.interpolate.interp2d kind 'linear', 'cubic' or 'quintic' 

332 # DEPRECATED since scipy 1.10, removed altogether in 1.14 

333 self._ev2d = spi.interp2d(xs, ys, hs, kind=self._k2interp2d[kind]) 

334 

335 else: # <https://scipy.GitHub.io/devdocs/tutorial/interpolate/interp_transition_guide.html> 

336 k = self._kxky(abs(kind)) 

337 # spi.RectBivariateSpline needs strictly ordered xs and ys 

338 r = spi.bisplrep(xs, ys, hs.T, kx=k, ky=k) 

339 

340 def _bisplev(x, y): 

341 return spi.bisplev(x, y, r) # .T 

342 

343 self._ev2d = _bisplev 

344 

345 except Exception as x: 

346 raise _SciPyIssue(x, self._ev_name) 

347 

348 def _kxky(self, kind): 

349 return Int_(kind=kind, low=1, high=5, Error=self._Error) 

350 

351 def _np_sp2(self, throwarnings=False): # PYCHOK no cover 

352 '''(INTERNAL) Import C{numpy} and C{scipy}, once. 

353 ''' 

354 # raise SciPyWarnings, but not if 

355 # scipy has already been imported 

356 if throwarnings: # PYCHOK no cover 

357 import sys 

358 if _scipy_ not in sys.modules: 

359 import warnings 

360 warnings.filterwarnings(_error_) 

361 return self.numpy, self.scipy 

362 

363 @property_ROver 

364 def numpy(self): 

365 '''Get the C{numpy} module or C{None}. 

366 ''' 

367 return _xnumpy(type(self), 1, 9) # overwrite property_ROver 

368 

369 @property_ROver 

370 def scipy(self): 

371 '''Get the C{scipy} module or C{None}. 

372 ''' 

373 return _xscipy(type(self), 1, 2) # overwrite property_ROver 

374 

375 @property_ROver 

376 def scipy_interpolate(self): 

377 '''Get the C{scipy.interpolate} module or C{None}. 

378 ''' 

379 _ = self.scipy 

380 import scipy.interpolate as spi # scipy 1.2.2 

381 return spi # overwrite property_ROver 

382 

383 def _scipy_version(self, **n): 

384 '''Get the C{scipy} version as 2- or 3-tuple C{(major, minor, micro)}. 

385 ''' 

386 return _MODS.internals._version2(self.scipy.version.version, **n) 

387 

388 def _xyhs3(self, knots, wrap=False, **name): 

389 # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals 

390 xs, ys, hs = zip(*_xyhs(knots=knots, wrap=wrap)) # PYCHOK yield 

391 n = len(hs) 

392 if n < self.kmin: 

393 raise _InsufficientError(self.kmin, knots=n) 

394 if name: 

395 self.name = name 

396 return map1(self.numpy.array, xs, ys, hs) 

397 

398 

399class HeightCubic(_HeightBase): 

400 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/ 

401 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>} 

402 C{kind='cubic'} or U{bisplrep/-ev<https://docs.SciPy.org/doc/scipy/ 

403 reference/generated/scipy.interpolate.interp2d.html>} C{kx=ky=3}. 

404 ''' 

405 _kind = -3 

406 _kmin = 16 

407 

408 def __init__(self, knots, **name_wrap): 

409 '''New L{HeightCubic} interpolator. 

410 

411 @arg knots: The points with known height (C{LatLon}s). 

412 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator (C{str}) 

413 and keyword argument C{b{wrap}=False} to wrap or I{normalize} all 

414 B{C{knots}} and B{C{llis}} locations iff C{True} (C{bool}). 

415 

416 @raise HeightError: Insufficient number of B{C{knots}} or invalid B{C{knot}}. 

417 

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

419 

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

421 

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

423 ''' 

424 xs_yx_hs = self._xyhs3(knots, **name_wrap) 

425 self._interp2d(*xs_yx_hs, kind=self._kind) 

426 

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

428 '''Interpolate the height for one or several locations. 

429 

430 @see: L{Here<_HeightBase.__call__>} for further details. 

431 ''' 

432 return self._evalls(llis, **wrap) 

433 

434 def _ev(self, yis, xis): # PYCHOK overwritten with .RectBivariateSpline.ev 

435 # to make SciPy .interp2d single (x, y) signature 

436 # match SciPy .ev signature(ys, xs), flipped multiples 

437 return map(self._ev2d, xis, yis) 

438 

439 

440class HeightLinear(HeightCubic): 

441 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/ 

442 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>} 

443 C{kind='linear'} or U{bisplrep/-ev<https://docs.SciPy.org/doc/scipy/ 

444 reference/generated/scipy.interpolate.interp2d.html>} C{kx=ky=1}. 

445 ''' 

446 _kind = -1 

447 _kmin = 2 

448 

449 def __init__(self, knots, **name_wrap): 

450 '''New L{HeightLinear} interpolator. 

451 

452 @see: L{Here<HeightCubic.__init__>} for all details. 

453 ''' 

454 HeightCubic.__init__(self, knots, **name_wrap) 

455 

456 if _FOR_DOCS: 

457 __call__ = HeightCubic.__call__ 

458 height = HeightCubic.height 

459 

460 

461class HeightLSQBiSpline(_HeightBase): 

462 '''Height interpolator using C{SciPy} U{LSQSphereBivariateSpline 

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

464 interpolate.LSQSphereBivariateSpline.html>}. 

465 ''' 

466 _kmin = 16 # k = 3, always 

467 

468 def __init__(self, knots, weight=None, low=1e-4, **name_wrap): 

469 '''New L{HeightLSQBiSpline} interpolator. 

470 

471 @arg knots: The points with known height (C{LatLon}s). 

472 @kwarg weight: Optional weight or weights for each B{C{knot}} 

473 (C{scalar} or C{scalar}s). 

474 @kwarg low: Optional lower bound for I{ordered knots} (C{radians}). 

475 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator 

476 (C{str}) and keyword argument C{b{wrap}=False} to wrap or 

477 I{normalize} all B{C{knots}} and B{C{llis}} locations iff 

478 C{True} (C{bool}). 

479 

480 @raise HeightError: Insufficient number of B{C{knots}} or an invalid 

481 B{C{knot}}, B{C{weight}} or B{C{eps}}. 

482 

483 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s. 

484 

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

486 installed. 

487 

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

489 

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

491 ''' 

492 np = self.numpy 

493 spi = self.scipy_interpolate 

494 

495 xs, ys, hs = self._xyhs3(knots, **name_wrap) 

496 n = len(hs) 

497 

498 w = weight 

499 if isscalar(w): 

500 w = float(w) 

501 if w <= 0: 

502 raise HeightError(weight=w) 

503 w = (w,) * n 

504 elif w is not None: 

505 m, w = len2(w) 

506 if m != n: 

507 raise LenError(HeightLSQBiSpline, weight=m, knots=n) 

508 m, i = min2(*map(float, w)) 

509 if m <= 0: # PYCHOK no cover 

510 raise HeightError(Fmt.INDEX(weight=i), m) 

511 try: 

512 if not EPS < low < (PI_2 - EPS): # 1e-4 like SciPy example 

513 raise HeightError(low=low) 

514 ps = np.array(_orderedup(xs, low, PI2 - low)) 

515 ts = np.array(_orderedup(ys, low, PI - low)) 

516 self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs, 

517 ts, ps, eps=EPS, w=w).ev 

518 except Exception as x: 

519 raise _SciPyIssue(x, self._ev_name) 

520 

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

522 '''Interpolate the height for one or several locations. 

523 

524 @see: L{Here<_HeightBase.__call__>} for further details. 

525 ''' 

526 return self._evalls(llis, **wrap) 

527 

528 

529class HeightSmoothBiSpline(_HeightBase): 

530 '''Height interpolator using C{SciPy} U{SmoothSphereBivariateSpline 

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

532 interpolate.SmoothSphereBivariateSpline.html>}. 

533 ''' 

534 _kmin = 16 # k = 3, always 

535 

536 def __init__(self, knots, s=4, **name_wrap): 

537 '''New L{HeightSmoothBiSpline} interpolator. 

538 

539 @arg knots: The points with known height (C{LatLon}s). 

540 @kwarg s: The spline smoothing factor (C{scalar}), default C{4}. 

541 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator 

542 (C{str}) and keyword argument C{b{wrap}=False} to wrap or 

543 I{normalize} all B{C{knots}} and B{C{llis}} locations iff 

544 C{True} (C{bool}). 

545 

546 @raise HeightError: Insufficient number of B{C{knots}} or an invalid 

547 B{C{knot}} or B{C{s}}. 

548 

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

550 installed. 

551 

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

553 

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

555 ''' 

556 spi = self.scipy_interpolate 

557 

558 s = Float_(smoothing=s, Error=HeightError, low=4) 

559 

560 xs, ys, hs = self._xyhs3(knots, **name_wrap) 

561 try: 

562 self._ev = spi.SmoothSphereBivariateSpline(ys, xs, hs, 

563 eps=EPS, s=s).ev 

564 except Exception as x: 

565 raise _SciPyIssue(x, self._ev_name) 

566 

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

568 '''Interpolate the height for one or several locations. 

569 

570 @see: L{Here<_HeightBase.__call__>} for further details. 

571 ''' 

572 return self._evalls(llis, **wrap) 

573 

574 

575class _HeightIDW(_HeightNamed): 

576 '''(INTERNAL) Base class for U{Inverse Distance Weighting 

577 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) height 

578 interpolators. 

579 

580 @see: U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/ 

581 geostatistics/Inverse-Distance-Weighting/index.html>}, 

582 U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/ 

583 shepard_interp_2d/shepard_interp_2d.html>} and other C{_HeightIDW*} 

584 classes. 

585 ''' 

586 _beta = 0 # fidw inverse power 

587 _func = None # formy function 

588 _knots = () # knots list or tuple 

589 _kwds = {} # func_ options 

590 

591 def __init__(self, knots, beta=2, **name__kwds): 

592 '''New C{_HeightIDW*} interpolator. 

593 

594 @arg knots: The points with known height (C{LatLon}s). 

595 @kwarg beta: Inverse distance power (C{int} 1, 2, or 3). 

596 @kwarg name__kwds: Optional C{B{name}=NN} for this height interpolator 

597 (C{str}) and any keyword arguments for the distance function, 

598 retrievable with property C{kwds}. 

599 

600 @raise HeightError: Insufficient number of B{C{knots}} or an invalid 

601 B{C{knot}} or B{C{beta}}. 

602 ''' 

603 name, kwds = _name2__(**name__kwds) 

604 if name: 

605 self.name = name 

606 

607 n, self._knots = len2(knots) 

608 if n < self.kmin: 

609 raise _InsufficientError(self.kmin, knots=n) 

610 self.beta = beta 

611 self._kwds = kwds or {} 

612 

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

614 '''Interpolate the height for one or several locations. 

615 

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

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

618 locations (C{bool}). 

619 

620 @return: A single interpolated height (C{float}) or a list 

621 or tuple of interpolated heights (C{float}s). 

622 

623 @raise HeightError: Insufficient number of B{C{llis}}, an 

624 invalid B{C{lli}} or L{pygeodesy.fidw} 

625 issue. 

626 ''' 

627 def _xy2(wrap=False): 

628 _w = _Wrap._latlonop(wrap) 

629 try: # like _xyhs above, but degrees 

630 for i, ll in enumerate(llis): 

631 yield _w(ll.lon, ll.lat) 

632 except Exception as x: 

633 i = Fmt.INDEX(llis=i) 

634 raise HeightError(i, ll, cause=x) 

635 

636 _as, llis = _as_llis2(llis) 

637 return _as(map(self._hIDW, *zip(*_xy2(**wrap)))) 

638 

639 @property_RO 

640 def adjust(self): 

641 '''Get the C{adjust} setting (C{bool}) or C{None}. 

642 ''' 

643 return _xkwds_get(self._kwds, adjust=None) 

644 

645 @property 

646 def beta(self): 

647 '''Get the inverse distance power (C{int}). 

648 ''' 

649 return self._beta 

650 

651 @beta.setter # PYCHOK setter! 

652 def beta(self, beta): 

653 '''Set the inverse distance power (C{int} 1, 2, or 3). 

654 

655 @raise HeightError: Invalid B{C{beta}}. 

656 ''' 

657 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3) 

658 

659 @property_RO 

660 def datum(self): 

661 '''Get the C{datum} setting or the default (L{Datum}). 

662 ''' 

663 return _xkwds_get(self._kwds, datum=self._datum) 

664 

665 def _datum_setter(self, datum): 

666 '''(INTERNAL) Set the default C{datum}. 

667 ''' 

668 d = datum or _xattr(self._knots[0], datum=None) 

669 if d and d is not self._datum: 

670 self._datum = _ellipsoidal_datum(d, name=self.name) 

671 

672 def _distances(self, x, y): 

673 '''(INTERNAL) Yield distances to C{(x, y)}. 

674 ''' 

675 _f, kwds = self._func, self._kwds 

676 if not callable(_f): # PYCHOK no cover 

677 self._notOverloaded(distance_function=_f) 

678 try: 

679 for i, k in enumerate(self._knots): 

680 yield _f(y, x, k.lat, k.lon, **kwds) 

681 except Exception as x: 

682 i = Fmt.INDEX(knots=i) 

683 raise HeightError(i, k, cause=x) 

684 

685 def _distancesTo(self, _To): 

686 '''(INTERNAL) Yield distances C{_To}. 

687 ''' 

688 try: 

689 for i, k in enumerate(self._knots): 

690 yield _To(k) 

691 except Exception as x: 

692 i = Fmt.INDEX(knots=i) 

693 raise HeightError(i, k, cause=x) 

694 

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

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

697 

698 @arg lats: Latitude or latitudes (each C{degrees}). 

699 @arg lons: Longitude or longitudes (each C{degrees}). 

700 @kwarg wrap: Keyword argument C{B{wrap}=False} (C{bool}). Use 

701 C{B{wrap}=True} to wrap or I{normalize} all B{C{lats}} 

702 and B{C{lons}}. 

703 

704 @return: A single interpolated height (C{float}) or a list of 

705 interpolated heights (each C{float}). 

706 

707 @raise HeightError: Insufficient or unequal number of B{C{lats}} 

708 and B{C{lons}} or a L{pygeodesy.fidw} issue. 

709 ''' 

710 lls = self._as_lls(lats, lons) # dup of _HeightBase.height 

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

712 

713 @Property_RO 

714 def _heights(self): 

715 '''(INTERNAL) Get the knots' heights. 

716 ''' 

717 return tuple(_xattr(k, height=0) for k in self.knots) 

718 

719 def _hIDW(self, x, y): 

720 '''(INTERNAL) Return the IDW-interpolated height at 

721 location (x, y), both C{degrees} or C{radians}. 

722 ''' 

723 ds, hs = self._distances(x, y), self._heights 

724 try: 

725 return _MODS.fmath.fidw(hs, ds, beta=self.beta) 

726 except (TypeError, ValueError) as e: 

727 raise HeightError(x=x, y=y, cause=e) 

728 

729 @property_RO 

730 def hypot(self): 

731 '''Get the C{hypot} setting (C{callable}) or C{None}. 

732 ''' 

733 return _xkwds_get(self._kwds, hypot=None) 

734 

735 @property_RO 

736 def knots(self): 

737 '''Get the B{C{knots}} (C{list} or C{tuple}). 

738 ''' 

739 return self._knots 

740 

741 @property_RO 

742 def kwds(self): 

743 '''Get the optional keyword arguments (C{dict}). 

744 ''' 

745 return self._kwds 

746 

747 @property_RO 

748 def limit(self): 

749 '''Get the C{limit} setting (C{degrees}) or C{None}. 

750 ''' 

751 return _xkwds_get(self._kwds, limit=None) 

752 

753 @property_RO 

754 def radius(self): 

755 '''Get the C{radius} setting (C{bool}) or C{None}. 

756 ''' 

757 return _xkwds_get(self._kwds, radius=None) 

758 

759 @property_RO 

760 def scaled(self): 

761 '''Get the C{scaled} setting (C{bool}) or C{None}. 

762 ''' 

763 return _xkwds_get(self._kwds, scaled=None) 

764 

765 @property_RO 

766 def wrap(self): 

767 '''Get the C{wrap} setting or the default (C{bool}) or C{None}. 

768 ''' 

769 return _xkwds_get(self._kwds, wrap=self._wrap) 

770 

771 

772class HeightIDWcosineLaw(_HeightIDW): 

773 '''Height interpolator using U{Inverse Distance Weighting 

774 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) 

775 and function L{pygeodesy.cosineLaw}. 

776 

777 @note: See note at function L{pygeodesy.vincentys_}. 

778 ''' 

779 def __init__(self, knots, beta=2, **name__corr_earth_datum_radius_wrap): 

780 '''New L{HeightIDWcosineLaw} interpolator. 

781 

782 @kwarg name__corr_earth_datum_radius_wrap: Optional C{B{name}=NN} 

783 for this height interpolator (C{str}) and any keyword 

784 arguments for function L{pygeodesy.cosineLaw}. 

785 

786 @see: L{Here<_HeightIDW.__init__>} for further details. 

787 ''' 

788 _HeightIDW.__init__(self, knots, beta=beta, **name__corr_earth_datum_radius_wrap) 

789 self._func = _formy.cosineLaw 

790 

791 if _FOR_DOCS: 

792 __call__ = _HeightIDW.__call__ 

793 height = _HeightIDW.height 

794 

795 

796class HeightIDWdistanceTo(_HeightIDW): 

797 '''Height interpolator using U{Inverse Distance Weighting 

798 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) 

799 and the points' C{LatLon.distanceTo} method. 

800 ''' 

801 def __init__(self, knots, beta=2, **name__distanceTo_kwds): 

802 '''New L{HeightIDWdistanceTo} interpolator. 

803 

804 @kwarg name__distanceTo_kwds: Optional C{B{name}=NN} for this 

805 height interpolator (C{str}) and keyword arguments 

806 for B{C{knots}}' method C{LatLon.distanceTo}. 

807 

808 @see: L{Here<_HeightIDW.__init__>} for further details. 

809 

810 @note: All B{C{points}} I{must} be instances of the same 

811 ellipsoidal or spherical C{LatLon} class, I{not 

812 checked}. 

813 ''' 

814 _HeightIDW.__init__(self, knots, beta=beta, **name__distanceTo_kwds) 

815 ks0 = _distanceTo(HeightError, knots=self._knots)[0] 

816 # use knots[0] class and datum to create compatible points 

817 # in ._as_lls instead of class LatLon_ and datum None 

818 self._datum = ks0.datum 

819 self._LLiC = ks0.classof # type(ks0) 

820 

821 def _distances(self, x, y): 

822 '''(INTERNAL) Yield distances to C{(x, y)}. 

823 ''' 

824 kwds, ll = self._kwds, self._LLiC(y, x) 

825 

826 def _To(k): 

827 return k.distanceTo(ll, **kwds) 

828 

829 return self._distancesTo(_To) 

830 

831 if _FOR_DOCS: 

832 __call__ = _HeightIDW.__call__ 

833 height = _HeightIDW.height 

834 

835 

836class HeightIDWequirectangular(_HeightIDW): 

837 '''Height interpolator using U{Inverse Distance Weighting 

838 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) 

839 and function L{pygeodesy.equirectangular4}. 

840 ''' 

841 def __init__(self, knots, beta=2, **name__adjust_limit_wrap): # XXX beta=1 

842 '''New L{HeightIDWequirectangular} interpolator. 

843 

844 @kwarg name__adjust_limit_wrap: Optional C{B{name}=NN} for this 

845 height interpolator (C{str}) and keyword arguments 

846 for function L{pygeodesy.equirectangular4}. 

847 

848 @see: L{Here<_HeightIDW.__init__>} for further details. 

849 ''' 

850 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_limit_wrap) 

851 

852 def _distances(self, x, y): 

853 '''(INTERNAL) Yield distances to C{(x, y)}. 

854 ''' 

855 _f, kwds = _formy.equirectangular4, self._kwds 

856 

857 def _To(k): 

858 return _f(y, x, k.lat, k.lon, **kwds).distance2 

859 

860 return self._distancesTo(_To) 

861 

862 if _FOR_DOCS: 

863 __call__ = _HeightIDW.__call__ 

864 height = _HeightIDW.height 

865 

866 

867class HeightIDWeuclidean(_HeightIDW): 

868 '''Height interpolator using U{Inverse Distance Weighting 

869 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) 

870 and function L{pygeodesy.euclidean_}. 

871 ''' 

872 def __init__(self, knots, beta=2, **name__adjust_radius_wrap): 

873 '''New L{HeightIDWeuclidean} interpolator. 

874 

875 @kwarg name__adjust_radius_wrap: Optional C{B{name}=NN} for this 

876 height interpolator (C{str}) and keyword arguments 

877 for function function L{pygeodesy.euclidean}. 

878 

879 @see: L{Here<_HeightIDW.__init__>} for further details. 

880 ''' 

881 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_radius_wrap) 

882 self._func = _formy.euclidean 

883 

884 if _FOR_DOCS: 

885 __call__ = _HeightIDW.__call__ 

886 height = _HeightIDW.height 

887 

888 

889class HeightIDWexact(_HeightIDW): 

890 '''Height interpolator using U{Inverse Distance Weighting