Coverage for pygeodesy/heights.py: 96%

361 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-12 11:45 -0400

1 

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

3 

4u'''Height interpolations of C{LatLon} points. 

5 

6Classes L{HeightCubic}, L{HeightIDWcosineAndoyerLambert}, 

7L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWcosineLaw}, 

8L{HeightIDWdistanceTo}, L{HeightIDWequirectangular}, L{HeightIDWeuclidean}, 

9L{HeightIDWflatLocal}, L{HeightIDWflatPolar}, L{HeightIDWhaversine}, 

10L{HeightIDWhubeny}, L{HeightIDWkarney}, L{HeightIDWthomas}, L{HeightIDWvincentys}, 

11L{HeightLinear}, L{HeightLSQBiSpline} and L{HeightSmoothBiSpline} 

12to interpolate the height of C{LatLon} locations or separate 

13lat-/longitudes from a set of C{LatLon} points with I{known heights}. 

14 

15Typical usage 

16============= 

17 

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

19called C{knots}. The C{knots} do not need to be ordered in any 

20particular way. 

21 

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

23 

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

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 other C{LatLon} location(s) with 

32 

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

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

35 

36or 

37 

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

39 

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

41 

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

43 

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

45 

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

47 

48or as 2 lists, 2 tuples, etc. 

49 

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

51 

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

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

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

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

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

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

58 

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

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

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

62 

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

64 Interpolation. 

65''' 

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

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

68 

69from pygeodesy.basics import isscalar, len2, map1, map2, _xnumpy, _xscipy 

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

71from pygeodesy.datums import _ellipsoidal_datum, _WGS84 

72from pygeodesy.errors import _AssertionError, LenError, PointsError, _SciPyIssue 

73from pygeodesy.fmath import fidw, hypot2 

74from pygeodesy.formy import cosineAndoyerLambert_, cosineForsytheAndoyerLambert_, \ 

75 cosineLaw_, euclidean_, flatPolar_, haversine_, \ 

76 _scale_rad, thomas_, vincentys_ 

77from pygeodesy.interns import NN, _COMMASPACE_, _cubic_, _datum_, _distanceTo_, \ 

78 _knots_, _linear_, _NOTEQUAL_, _scipy_, _SPACE_, _STAR_ 

79from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS 

80from pygeodesy.named import _Named, notOverloaded 

81from pygeodesy.points import LatLon_, Property_RO 

82# from pygeodesy.props import Property_RO # from .points 

83from pygeodesy.streprs import _boolkwds, Fmt 

84from pygeodesy.units import Float_, Int_ 

85from pygeodesy.utily import radiansPI, radiansPI2, unrollPI 

86 

87__all__ = _ALL_LAZY.heights 

88__version__ = '23.04.11' 

89 

90_error_ = 'error' 

91_insufficient_ = 'insufficient' 

92_llis_ = 'llis' 

93_smoothing_ = 'smoothing' 

94 

95 

96class HeightError(PointsError): 

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

98 ''' 

99 pass 

100 

101 

102def _alist(ais): 

103 # return list of floats, not numpy.float64s 

104 return list(map(float, ais)) 

105 

106 

107def _allis2(llis, m=1, Error=HeightError): # imported by .geoids 

108 # dtermine return type and convert lli C{LatLon}s to list 

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

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

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

112 

113 n = len(llis) 

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

115 llis = llis[0] 

116 try: 

117 n, llis = len2(llis) 

118 _as = _alist # return list of interpolated heights 

119 except TypeError: # single lli 

120 n, llis = 1, [llis] 

121 _as = _ascalar # return single interpolated heights 

122 else: # of 0, 2 or more llis 

123 _as = _atuple # return tuple of interpolated heights 

124 

125 if n < m: 

126 raise _insufficientError(m, Error=Error, llis=n) 

127 return _as, llis 

128 

129 

130def _ascalar(ais): # imported by .geoids 

131 # return single float, not numpy.float64 

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

133 if len(ais) != 1: 

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

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

136 raise _AssertionError(n, txt=t) 

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

138 

139 

140def _atuple(ais): 

141 # return tuple of floats, not numpy.float64s 

142 return tuple(map(float, ais)) 

143 

144 

145def _axyllis4(atype, llis, m=1, off=True): 

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

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

148 _as, llis = _allis2(llis, m=m) 

149 xis, yis, _ = zip(*_xyhs(llis, off=off)) # PYCHOK unzip 

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

151 

152 

153def _insufficientError(need, Error=HeightError, **name_value): # PYCHOK no cover 

154 # create an insufficient Error instance 

155 t = _COMMASPACE_(_insufficient_, str(need)) 

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

157 

158 

159def _ordedup(ts, lo=EPS, hi=PI2-EPS): 

160 # clip, order and remove duplicates 

161 # p, ks = 0, [] 

162 # for k in sorted(max(lo, min(hi, t)) for t in ts): 

163 # if k > p: 

164 # ks.append(k) 

165 # p = k 

166 # return ks 

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

168 

169 

170def _xyhs(lls, off=True, name=_llis_): 

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

172 # x: 0 <= lon <= PI2, y: 0 <= lat <= PI if off is True 

173 # else x: -PI <= lon <= PI, y: -PI_2 <= lat <= PI_2 

174 if off: 

175 xf = yf = _0_0 

176 else: # undo offset 

177 xf, yf = PI, PI_2 

178 try: 

179 for i, ll in enumerate(lls): 

180 yield (max(_0_0, radiansPI2(ll.lon + _180_0)) - xf), \ 

181 (max(_0_0, radiansPI( ll.lat + _90_0)) - yf), ll.height 

182 except AttributeError as x: 

183 i = Fmt.SQUARE(name, i) 

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

185 

186 

187def _xyhs3(atype, m, knots, off=True): 

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

189 xs, ys, hs = zip(*_xyhs(knots, off=off, name=_knots_)) # PYCHOK unzip 

190 n = len(hs) 

191 if n < m: 

192 raise _insufficientError(m, knots=n) 

193 return map1(atype, xs, ys, hs) 

194 

195 

196class _HeightBase(_Named): # imported by .geoids 

197 '''(INTERNAL) Interpolator base class. 

198 ''' 

199 _adjust = None # not applicable 

200 _datum = None # ._height datum 

201 _kmin = 2 # min number of knots 

202 _LLis = LatLon_ # ._height class 

203 _np_sp = None # (numpy, scipy) 

204 _wrap = None # not applicable 

205 

206 def __call__(self, *args): # PYCHOK no cover 

207 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 

208 ''' 

209 notOverloaded(self, callername='__call__', *args) 

210 

211 @Property_RO 

212 def adjust(self): 

213 '''Get the adjust setting (C{bool} or C{None} if not applicable). 

214 ''' 

215 return self._adjust 

216 

217 def _axyllis4(self, llis): 

218 return _axyllis4(self.numpy.array, llis) 

219 

220 @Property_RO 

221 def datum(self): 

222 '''Get the datum (L{Datum} or C{None} if not applicable). 

223 ''' 

224 return self._datum 

225 

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

227 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 

228 ''' 

229 notOverloaded(self, *args) 

230 

231 def _eval(self, llis): # XXX single arg, not *args 

232 _as, xis, yis, _ = self._axyllis4(llis) 

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

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

235 except Exception as x: 

236 raise _SciPyIssue(x) 

237 

238 def _height(self, lats, lons, Error=HeightError): 

239 LLis, d = self._LLis, self.datum 

240 if isscalar(lats) and isscalar(lons): 

241 llis = LLis(lats, lons, datum=d) 

242 else: 

243 n, lats = len2(lats) 

244 m, lons = len2(lons) 

245 if n != m: 

246 # format a LenError, but raise an Error 

247 e = LenError(self.__class__, lats=n, lons=m, txt=None) 

248 raise e if Error is LenError else Error(str(e)) 

249 llis = [LLis(*ll, datum=d) for ll in zip(lats, lons)] 

250 return self(llis) # __call__(lli) or __call__(llis) 

251 

252 @Property_RO 

253 def kmin(self): 

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

255 ''' 

256 return self._kmin 

257 

258 def _np_sp2(self, throwarnings=False): 

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

260 ''' 

261 t = _HeightBase._np_sp 

262 if not t: 

263 # raise SciPyWarnings, but not if 

264 # scipy has been imported already 

265 if throwarnings: # PYCHOK no cover 

266 import sys 

267 if _scipy_ not in sys.modules: 

268 import warnings 

269 warnings.filterwarnings(_error_) 

270 

271 sp = _xscipy(self.__class__, 1, 2) 

272 np = _xnumpy(self.__class__, 1, 9) 

273 

274 _HeightBase._np_sp = t = np, sp 

275 return t 

276 

277 @Property_RO 

278 def numpy(self): 

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

280 ''' 

281 np, _ = self._np_sp2() 

282 return np 

283 

284 @Property_RO 

285 def scipy(self): 

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

287 ''' 

288 _, sp = self._np_sp2() 

289 return sp 

290 

291 @Property_RO 

292 def scipy_interpolate(self): 

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

294 ''' 

295 _ = self.scipy 

296 import scipy.interpolate as spi 

297 return spi 

298 

299 def _xyhs3(self, knots): 

300 return _xyhs3(self.numpy.array, self._kmin, knots) 

301 

302 @Property_RO 

303 def wrap(self): 

304 '''Get the wrap setting (C{bool} or C{None} if not applicable). 

305 ''' 

306 return self._wrap 

307 

308 

309class HeightCubic(_HeightBase): 

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

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

312 C{kind='cubic'}. 

313 ''' 

314 _interp2d = None 

315 _kind = _cubic_ 

316 _kmin = 16 

317 

318 def __init__(self, knots, name=NN): 

319 '''New L{HeightCubic} interpolator. 

320 

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

322 @kwarg name: Optional name for this height interpolator (C{str}). 

323 

324 @raise HeightError: Insufficient number of B{C{knots}} or 

325 invalid B{C{knot}}. 

326 

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

328 or not installed. 

329 

330 @raise SciPyError: A C{scipy.interpolate.interp2d} issue. 

331 

332 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning 

333 as exception. 

334 ''' 

335 spi = self.scipy_interpolate 

336 

337 xs, ys, hs = self._xyhs3(knots) 

338 try: # SciPy.interpolate.interp2d kind 'linear' or 'cubic' 

339 self._interp2d = spi.interp2d(xs, ys, hs, kind=self._kind) 

340 except Exception as x: 

341 raise _SciPyIssue(x) 

342 

343 if name: 

344 self.name = name 

345 

346 def __call__(self, *llis): 

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

348 

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

350 C{LatLon}s). 

351 

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

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

354 

355 @raise HeightError: Insufficient number of B{C{llis}} or 

356 an invalid B{C{lli}}. 

357 

358 @raise SciPyError: A C{scipy.interpolate.interp2d} issue. 

359 

360 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning 

361 as exception. 

362 ''' 

363 return _HeightBase._eval(self, llis) 

364 

365 def _ev(self, yis, xis): # PYCHOK expected 

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

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

368 return map(self._interp2d, xis, yis) 

369 

370 def height(self, lats, lons): 

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

372 

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

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

375 

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

377 interpolated heights (C{float}s). 

378 

379 @raise HeightError: Insufficient or non-matching number of 

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

381 

382 @raise SciPyError: A C{scipy.interpolate.interp2d} issue. 

383 

384 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning 

385 as exception. 

386 ''' 

387 return _HeightBase._height(self, lats, lons) 

388 

389 

390class HeightLinear(HeightCubic): 

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

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

393 C{kind='linear'}. 

394 ''' 

395 _kind = _linear_ 

396 _kmin = 2 

397 

398 def __init__(self, knots, name=NN): 

399 '''New L{HeightLinear} interpolator. 

400 

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

402 @kwarg name: Optional name for this height interpolator (C{str}). 

403 

404 @raise HeightError: Insufficient number of B{C{knots}} or 

405 an invalid B{C{knot}}. 

406 

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

408 or not installed. 

409 

410 @raise SciPyError: A C{scipy.interpolate.interp2d} issue. 

411 

412 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning 

413 as exception. 

414 ''' 

415 HeightCubic.__init__(self, knots, name=name) 

416 

417 if _FOR_DOCS: 

418 __call__ = HeightCubic.__call__ 

419 height = HeightCubic.height 

420 

421 

422class _HeightIDW(_HeightBase): 

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

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

425 height interpolators. 

426 ''' 

427 _beta = 0 # inverse distance power 

428 _hs = () # known heights 

429 _xs = () # knot lons 

430 _ys = () # knot lats 

431 

432 def __init__(self, knots, beta=2, name=NN, **wrap_adjust): 

433 '''New L{_HeightIDW} interpolator. 

434 ''' 

435 self._xs, self._ys, self._hs = _xyhs3(tuple, self._kmin, knots, off=False) 

436 self.beta = beta 

437 if name: 

438 self.name = name 

439 if wrap_adjust: 

440 _boolkwds(self, **wrap_adjust) 

441 

442 def __call__(self, *llis): 

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

444 

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

446 C{LatLon}s). 

447 

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

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

450 

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

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

453 issue. 

454 ''' 

455 _as, xis, yis, _ = _axyllis4(tuple, llis, off=False) 

456 return _as(map(self._hIDW, xis, yis)) 

457 

458 def _datum_setter(self, datum, knots): 

459 '''(INTERNAL) Set the datum. 

460 

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

462 ''' 

463 d = datum or getattr(knots[0], _datum_, datum) 

464 if d not in (None, self._datum): 

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

466 

467 def _distances(self, x, y): # PYCHOK unused (x, y) radians 

468 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 

469 ''' 

470 notOverloaded(self, x, y) 

471 

472 def _distances_angular_(self, func_, x, y): 

473 # return angular distances from func_ 

474 for xk, yk in zip(self._xs, self._ys): 

475 r, _ = unrollPI(xk, x, wrap=self._wrap) 

476 yield func_(yk, y, r) 

477 

478 def _distances_angular_datum_(self, func_, x, y): 

479 # return angular distances from func_ 

480 for xk, yk in zip(self._xs, self._ys): 

481 r, _ = unrollPI(xk, x, wrap=self._wrap) 

482 yield func_(yk, y, r, datum=self._datum) 

483 

484 def _hIDW(self, x, y): 

485 # interpolate height at (x, y) radians or degrees 

486 try: 

487 ds = self._distances(x, y) 

488 return fidw(self._hs, ds, beta=self._beta) 

489 except (TypeError, ValueError) as x: 

490 raise HeightError(str(x), cause=x) 

491 

492 @property 

493 def beta(self): 

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

495 ''' 

496 return self._beta 

497 

498 @beta.setter # PYCHOK setter! 

499 def beta(self, beta): 

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

501 

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

503 ''' 

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

505 

506 def height(self, lats, lons): 

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

508 

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

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

511 

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

513 interpolated heights (C{float}s). 

514 

515 @raise HeightError: Insufficient or non-matching number of 

516 B{C{lats}} and B{C{lons}} or L{pygeodesy.fidw} 

517 issue. 

518 ''' 

519 return _HeightBase._height(self, lats, lons) 

520 

521 

522class HeightIDWcosineAndoyerLambert(_HeightIDW): 

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

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

525 and the I{angular} distance in C{radians} from function 

526 L{pygeodesy.cosineAndoyerLambert_}. 

527 

528 @see: L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWdistanceTo}, 

529 L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas}, 

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

531 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

533 shepard_interp_2d/shepard_interp_2d.html>}. 

534 ''' 

535 _datum = _WGS84 

536 _wrap = False 

537 

538 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN): 

539 '''New L{HeightIDWcosineAndoyerLambert} interpolator. 

540 

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

542 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 

543 and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid}, 

544 L{Ellipsoid2} or L{a_f2Tuple}). 

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

546 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

547 @kwarg name: Optional name for this height interpolator (C{str}). 

548 

549 @raise HeightError: Insufficient number of B{C{knots}} or 

550 an invalid B{C{knot}} or B{C{beta}}. 

551 

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

553 ''' 

554 _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap) 

555 self._datum_setter(datum, knots) 

556 

557 def _distances(self, x, y): # (x, y) radians 

558 return self._distances_angular_datum_(cosineAndoyerLambert_, x, y) 

559 

560 if _FOR_DOCS: 

561 __call__ = _HeightIDW.__call__ 

562 height = _HeightIDW.height 

563 

564 

565class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW): 

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

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

568 and the I{angular} distance in C{radians} from function 

569 L{pygeodesy.cosineForsytheAndoyerLambert_}. 

570 

571 @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWdistanceTo}, 

572 L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas}, 

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

574 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

576 shepard_interp_2d/shepard_interp_2d.html>}. 

577 ''' 

578 _datum = _WGS84 

579 _wrap = False 

580 

581 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN): 

582 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator. 

583 

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

585 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 

586 and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid}, 

587 L{Ellipsoid2} or L{a_f2Tuple}). 

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

589 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

590 @kwarg name: Optional name for this height interpolator (C{str}). 

591 

592 @raise HeightError: Insufficient number of B{C{knots}} or 

593 an invalid B{C{knot}} or B{C{beta}}. 

594 

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

596 ''' 

597 _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap) 

598 self._datum_setter(datum, knots) 

599 

600 def _distances(self, x, y): # (x, y) radians 

601 return self._distances_angular_datum_(cosineForsytheAndoyerLambert_, x, y) 

602 

603 if _FOR_DOCS: 

604 __call__ = _HeightIDW.__call__ 

605 height = _HeightIDW.height 

606 

607 

608class HeightIDWcosineLaw(_HeightIDW): 

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

610 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the 

611 I{angular} distance in C{radians} from function L{pygeodesy.cosineLaw_}. 

612 

613 @see: L{HeightIDWequirectangular}, L{HeightIDWeuclidean}, 

614 L{HeightIDWflatPolar}, L{HeightIDWhaversine}, L{HeightIDWvincentys}, 

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

616 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

618 shepard_interp_2d/shepard_interp_2d.html>}. 

619 

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

621 ''' 

622 _wrap = False 

623 

624 def __init__(self, knots, beta=2, wrap=False, name=NN): 

625 '''New L{HeightIDWcosineLaw} interpolator. 

626 

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

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

629 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

630 @kwarg name: Optional name for this height interpolator (C{str}). 

631 

632 @raise HeightError: Insufficient number of B{C{knots}} or 

633 an invalid B{C{knot}} or B{C{beta}}. 

634 ''' 

635 _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap) 

636 

637 def _distances(self, x, y): # (x, y) radians 

638 return self._distances_angular_(cosineLaw_, x, y) 

639 

640 if _FOR_DOCS: 

641 __call__ = _HeightIDW.__call__ 

642 height = _HeightIDW.height 

643 

644 

645class HeightIDWdistanceTo(_HeightIDW): 

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

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

648 and the distance from the points' C{LatLon.distanceTo} method, 

649 conventionally in C{meter}. 

650 

651 @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert}, 

652 L{HeightIDWflatPolar}, L{HeightIDWkarney}, L{HeightIDWthomas}, 

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

654 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

656 shepard_interp_2d/shepard_interp_2d.html>}. 

657 ''' 

658 _distanceTo_kwds = {} 

659 _ks = () 

660 

661 def __init__(self, knots, beta=2, name=NN, **distanceTo_kwds): 

662 '''New L{HeightIDWdistanceTo} interpolator. 

663 

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

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

666 @kwarg name: Optional name for this height interpolator (C{str}). 

667 @kwarg distanceTo_kwds: Optional keyword arguments for the 

668 B{C{points}}' C{LatLon.distanceTo} 

669 method. 

670 

671 @raise HeightError: Insufficient number of B{C{knots}} or 

672 an invalid B{C{knot}} or B{C{beta}}. 

673 

674 @raise ImportError: Package U{geographiclib 

675 <https://PyPI.org/project/geographiclib>} missing 

676 iff B{C{points}} are L{ellipsoidalKarney.LatLon}s. 

677 

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

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

680 checked however}. 

681 ''' 

682 n, self._ks = len2(knots) 

683 if n < self._kmin: 

684 raise _insufficientError(self._kmin, knots=n) 

685 for i, k in enumerate(self._ks): 

686 if not callable(getattr(k, _distanceTo_, None)): 

687 i = Fmt.SQUARE(_knots_, i) 

688 raise HeightError(i, k, txt=_distanceTo_) 

689 

690 # use knots[0] class and datum to create 

691 # compatible points in _HeightBase._height 

692 # instead of class LatLon_ and datum None 

693 self._datum = self._ks[0].datum 

694 self._LLis = self._ks[0].classof 

695 

696 self.beta = beta 

697 if name: 

698 self.name = name 

699 if distanceTo_kwds: 

700 self._distanceTo_kwds = distanceTo_kwds 

701 

702 def __call__(self, *llis): 

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

704 

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

706 C{LatLon}s). 

707 

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

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

710 

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

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

713 issue. 

714 ''' 

715 _as, llis = _allis2(llis) 

716 return _as(map(self._hIDW, llis)) 

717 

718 if _FOR_DOCS: 

719 height = _HeightIDW.height 

720 

721 def _hIDW(self, lli): # PYCHOK expected 

722 # interpolate height at point lli 

723 try: 

724 kwds = self._distanceTo_kwds 

725 ds = (k.distanceTo(lli, **kwds) for k in self._ks) 

726 return fidw(self._hs, ds, beta=self._beta) 

727 except (TypeError, ValueError) as x: 

728 raise HeightError(str(x)) 

729 

730 @property # NOT _RO, no caching 

731 def _hs(self): # see HeightIDWkarney 

732 for k in self._ks: 

733 yield k.height 

734 

735 

736class HeightIDWequirectangular(_HeightIDW): 

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

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

739 the I{angular} distance in C{radians squared} like function 

740 L{pygeodesy.equirectangular_}. 

741 

742 @see: L{HeightIDWeuclidean}, L{HeightIDWflatPolar}, 

743 L{HeightIDWhaversine}, L{HeightIDWvincentys}, 

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

745 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

747 shepard_interp_2d/shepard_interp_2d.html>}. 

748 ''' 

749 _adjust = True 

750 _wrap = False 

751 

752 def __init__(self, knots, adjust=True, wrap=False, name=NN): 

753 '''New L{HeightIDWequirectangular} interpolator. 

754 

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

756 @kwarg adjust: Adjust the wrapped, unrolled longitudinal 

757 delta by the cosine of the mean latitude (C{bool}). 

758 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

759 @kwarg name: Optional name for this height interpolator (C{str}). 

760 

761 @raise HeightError: Insufficient number of B{C{knots}} or 

762 an invalid B{C{knot}}. 

763 ''' 

764 _HeightIDW.__init__(self, knots, beta=1, name=name, wrap=wrap, 

765 adjust=adjust) 

766 

767 def _distances(self, x, y): # (x, y) radians**2 

768 for xk, yk in zip(self._xs, self._ys): 

769 d, _ = unrollPI(xk, x, wrap=self._wrap) 

770 if self._adjust: 

771 d *= _scale_rad(yk, y) 

772 yield hypot2(d, yk - y) # like equirectangular_ distance2 

773 

774 if _FOR_DOCS: 

775 __call__ = _HeightIDW.__call__ 

776 height = _HeightIDW.height 

777 

778 

779class HeightIDWeuclidean(_HeightIDW): 

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

781 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the 

782 I{angular} distance in C{radians} from function L{pygeodesy.euclidean_}. 

783 

784 @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular}, 

785 L{HeightIDWhaversine}, L{HeightIDWvincentys}, 

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

787 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

789 shepard_interp_2d/shepard_interp_2d.html>}. 

790 ''' 

791 _adjust = True 

792 

793 def __init__(self, knots, adjust=True, beta=2, name=NN): 

794 '''New L{HeightIDWeuclidean} interpolator. 

795 

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

797 @kwarg adjust: Adjust the longitudinal delta by the cosine 

798 of the mean latitude for B{C{adjust}}=C{True}. 

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

800 @kwarg name: Optional name for this height interpolator (C{str}). 

801 

802 @raise HeightError: Insufficient number of B{C{knots}} or 

803 an invalid B{C{knot}} or B{C{beta}}. 

804 ''' 

805 _HeightIDW.__init__(self, knots, beta=beta, name=name, adjust=adjust) 

806 

807 def _distances(self, x, y): # (x, y) radians 

808 for xk, yk in zip(self._xs, self._ys): 

809 yield euclidean_(yk, y, xk - x, adjust=self._adjust) 

810 

811 if _FOR_DOCS: 

812 __call__ = _HeightIDW.__call__ 

813 height = _HeightIDW.height 

814 

815 

816class HeightIDWflatLocal(_HeightIDW): 

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

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

819 and the I{angular} distance in C{radians squared} like function 

820 L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}. 

821 

822 @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert}, 

823 L{HeightIDWdistanceTo}, L{HeightIDWhubeny}, L{HeightIDWthomas}, 

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

825 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

827 shepard_interp_2d/shepard_interp_2d.html>}. 

828 ''' 

829 _datum = _WGS84 

830 _wrap = False 

831 

832 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN): 

833 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator. 

834 

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

836 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 

837 and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid}, 

838 L{Ellipsoid2} or L{a_f2Tuple}). 

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

840 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

841 @kwarg name: Optional name for this height interpolator (C{str}). 

842 

843 @raise HeightError: Insufficient number of B{C{knots}} or 

844 an invalid B{C{knot}} or B{C{beta}}. 

845 

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

847 ''' 

848 _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap) 

849 self._datum_setter(datum, knots) 

850 

851 def _distances(self, x, y): # (x, y) radians 

852 return self._distances_angular_(self._hubeny_2, x, y) # radians**2 

853 

854 @Property_RO 

855 def _hubeny_2(self): 

856 '''(INTERNAL) Get and cache the C{.datum.ellipsoid._hubeny_2} method. 

857 ''' 

858 return self.datum.ellipsoid._hubeny_2 

859 

860 if _FOR_DOCS: 

861 __call__ = _HeightIDW.__call__ 

862 height = _HeightIDW.height 

863 

864 

865class HeightIDWflatPolar(_HeightIDW): 

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

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

868 and the I{angular} distance in C{radians} from function 

869 L{pygeodesy.flatPolar_}. 

870 

871 @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular}, 

872 L{HeightIDWeuclidean}, L{HeightIDWhaversine}, L{HeightIDWvincentys}, 

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

874 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

876 shepard_interp_2d/shepard_interp_2d.html>}. 

877 ''' 

878 _wrap = False 

879 

880 def __init__(self, knots, beta=2, wrap=False, name=NN): 

881 '''New L{HeightIDWflatPolar} interpolator. 

882 

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

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

885 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

886 @kwarg name: Optional name for this height interpolator (C{str}). 

887 

888 @raise HeightError: Insufficient number of B{C{knots}} or 

889 an invalid B{C{knot}} or B{C{beta}}. 

890 ''' 

891 _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap) 

892 

893 def _distances(self, x, y): # (x, y) radians 

894 return self._distances_angular_(flatPolar_, x, y) 

895 

896 if _FOR_DOCS: 

897 __call__ = _HeightIDW.__call__ 

898 height = _HeightIDW.height 

899 

900 

901class HeightIDWhaversine(_HeightIDW): 

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

903 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the 

904 I{angular} distance in C{radians} from function L{pygeodesy.haversine_}. 

905 

906 @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular}, L{HeightIDWeuclidean}, 

907 L{HeightIDWflatPolar}, L{HeightIDWvincentys}, 

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

909 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

911 shepard_interp_2d/shepard_interp_2d.html>}. 

912 

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

914 ''' 

915 _wrap = False 

916 

917 def __init__(self, knots, beta=2, wrap=False, name=NN): 

918 '''New L{HeightIDWhaversine} interpolator. 

919 

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

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

922 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

923 @kwarg name: Optional name for this height interpolator (C{str}). 

924 

925 @raise HeightError: Insufficient number of B{C{knots}} or 

926 an B{C{knot}} or B{C{beta}}. 

927 ''' 

928 _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap) 

929 

930 def _distances(self, x, y): # (x, y) radians 

931 return self._distances_angular_(haversine_, x, y) 

932 

933 if _FOR_DOCS: 

934 __call__ = _HeightIDW.__call__ 

935 height = _HeightIDW.height 

936 

937 

938class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny 

939 if _FOR_DOCS: 

940 __doc__ = HeightIDWflatLocal.__doc__ 

941 __init__ = HeightIDWflatLocal.__init__ 

942 __call__ = HeightIDWflatLocal.__call__ 

943 height = HeightIDWflatLocal.height 

944 

945 

946class HeightIDWkarney(_HeightIDW): 

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

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

949 the I{angular} distance in C{degrees} from I{Karney}'s 

950 U{geographiclib<https://PyPI.org/project/geographiclib>} U{Geodesic 

951 <https://GeographicLib.SourceForge.io/Python/doc/code.html>} 

952 Inverse method. 

953 

954 @see: L{HeightIDWcosineAndoyerLambert}, 

955 L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWdistanceTo}, 

956 L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas}, 

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

958 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

960 shepard_interp_2d/shepard_interp_2d.html>}. 

961 ''' 

962 _datum = _WGS84 

963 _Inverse1 = None 

964 _ks = () 

965 _wrap = False 

966 

967 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN): 

968 '''New L{HeightIDWkarney} interpolator. 

969 

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

971 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 

972 and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid}, 

973 L{Ellipsoid2} or L{a_f2Tuple}). 

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

975 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}). 

976 @kwarg name: Optional name for this height interpolator (C{str}). 

977 

978 @raise HeightError: Insufficient number of B{C{knots}} or 

979 an invalid B{C{knot}}, B{C{datum}} or 

980 B{C{beta}}. 

981 

982 @raise ImportError: Package U{geographiclib 

983 <https://PyPI.org/project/geographiclib>} missing. 

984 

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

986 ''' 

987 n, self._ks = len2(knots) 

988 if n < self._kmin: 

989 raise _insufficientError(self._kmin, knots=n) 

990 self._datum_setter(datum, self._ks) 

991 self._Inverse1 = self.datum.ellipsoid.geodesic.Inverse1 

992 

993 self.beta = beta 

994 if wrap: 

995 self._wrap = True 

996 if name: 

997 self.name = name 

998 

999 def _distances(self, x, y): # (x, y) degrees 

1000 for k in self._ks: 

1001 # non-negative I{angular} distance in C{degrees} 

1002 yield self._Inverse1(y, x, k.lat, k.lon, wrap=self._wrap) 

1003 

1004 @property # NOT _RO, no caching 

1005 def _hs(self): # see HeightIDWdistanceTo 

1006 for k in self._ks: 

1007 yield k.height 

1008 

1009 def __call__(self, *llis): 

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

1011 

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

1013 C{LatLon}s). 

1014 

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

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

1017 

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

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

1020 issue. 

1021 ''' 

1022 def _xy2(lls): 

1023 try: # like _xyhs above, but keeping degrees 

1024 for i, ll in enumerate(lls): 

1025 yield ll.lon, ll.lat 

1026 except AttributeError as x: 

1027 i = Fmt.SQUARE(llis=i) 

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

1029 

1030 _as, llis = _allis2(llis) 

1031 return _as(map(self._hIDW, *zip(*_xy2(llis)))) 

1032 

1033 if _FOR_DOCS: 

1034 height = _HeightIDW.height 

1035 

1036 

1037class HeightIDWthomas(_HeightIDW): 

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

1039 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the 

1040 I{angular} distance in C{radians} from function L{pygeodesy.thomas_}. 

1041 

1042 @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert}, 

1043 L{HeightIDWdistanceTo}, L{HeightIDWflatLocal}, L{HeightIDWhubeny}, 

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

1045 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

1047 shepard_interp_2d/shepard_interp_2d.html>}. 

1048 ''' 

1049 _datum = _WGS84 

1050 _wrap = False 

1051 

1052 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN): 

1053 '''New L{HeightIDWthomas} interpolator. 

1054 

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

1056 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 

1057 and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid}, 

1058 L{Ellipsoid2} or L{a_f2Tuple}). 

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

1060 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

1061 @kwarg name: Optional name for this height interpolator (C{str}). 

1062 

1063 @raise HeightError: Insufficient number of B{C{knots}} or 

1064 an invalid B{C{knot}} or B{C{beta}}. 

1065 

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

1067 ''' 

1068 _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap) 

1069 self._datum_setter(datum, knots) 

1070 

1071 def _distances(self, x, y): # (x, y) radians 

1072 return self._distances_angular_datum_(thomas_, x, y) 

1073 

1074 if _FOR_DOCS: 

1075 __call__ = _HeightIDW.__call__ 

1076 height = _HeightIDW.height 

1077 

1078 

1079class HeightIDWvincentys(_HeightIDW): 

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

1081 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the 

1082 I{angular} distance in C{radians} from function L{pygeodesy.vincentys_}. 

1083 

1084 @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular}, 

1085 L{HeightIDWeuclidean}, L{HeightIDWflatPolar}, L{HeightIDWhaversine}, 

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

1087 geostatistics/Inverse-Distance-Weighting/index.html>} and 

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

1089 shepard_interp_2d/shepard_interp_2d.html>}. 

1090 

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

1092 ''' 

1093 _wrap = False 

1094 

1095 def __init__(self, knots, beta=2, wrap=False, name=NN): 

1096 '''New L{HeightIDWvincentys} interpolator. 

1097 

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

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

1100 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

1101 @kwarg name: Optional name for this height interpolator (C{str}). 

1102 

1103 @raise HeightError: Insufficient number of B{C{knots}} or 

1104 an invalid B{C{knot}} or B{C{beta}}. 

1105 ''' 

1106 _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap) 

1107 

1108 def _distances(self, x, y): # (x, y) radians 

1109 return self._distances_angular_(vincentys_, x, y) 

1110 

1111 if _FOR_DOCS: 

1112 __call__ = _HeightIDW.__call__ 

1113 height = _HeightIDW.height 

1114 

1115 

1116class HeightLSQBiSpline(_HeightBase): 

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

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

1119 interpolate.LSQSphereBivariateSpline.html>}. 

1120 ''' 

1121 _kmin = 16 # k = 3, always 

1122 

1123 def __init__(self, knots, weight=None, name=NN): 

1124 '''New L{HeightLSQBiSpline} interpolator. 

1125 

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

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

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

1129 @kwarg name: Optional name for this height interpolator (C{str}). 

1130 

1131 @raise HeightError: Insufficient number of B{C{knots}} or 

1132 an invalid B{C{knot}} or B{C{weight}}. 

1133 

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

1135 

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

1137 or not installed. 

1138 

1139 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue. 

1140 

1141 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline} 

1142 warning as exception. 

1143 ''' 

1144 np = self.numpy 

1145 spi = self.scipy_interpolate 

1146 

1147 xs, ys, hs = self._xyhs3(knots) 

1148 n = len(hs) 

1149 

1150 w = weight 

1151 if isscalar(w): 

1152 w = float(w) 

1153 if w <= 0: 

1154 raise HeightError(weight=w) 

1155 w = [w] * n 

1156 elif w is not None: 

1157 m, w = len2(w) 

1158 if m != n: 

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

1160 w = map2(float, w) 

1161 m = min(w) 

1162 if m <= 0: # PYCHOK no cover 

1163 w = Fmt.SQUARE(weight=w.find(m)) 

1164 raise HeightError(w, m) 

1165 try: 

1166 T = 1.0e-4 # like SciPy example 

1167 ps = np.array(_ordedup(xs, T, PI2 - T)) 

1168 ts = np.array(_ordedup(ys, T, PI - T)) 

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

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

1171 except Exception as x: 

1172 raise _SciPyIssue(x) 

1173 

1174 if name: 

1175 self.name = name 

1176 

1177 def __call__(self, *llis): 

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

1179 

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

1181 C{LatLon}s). 

1182 

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

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

1185 

1186 @raise HeightError: Insufficient number of B{C{llis}} or 

1187 an invalid B{C{lli}}. 

1188 

1189 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue. 

1190 

1191 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline} 

1192 warning as exception. 

1193 ''' 

1194 return _HeightBase._eval(self, llis) 

1195 

1196 def height(self, lats, lons): 

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

1198 

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

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

1201 

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

1203 interpolated heights (C{float}s). 

1204 

1205 @raise HeightError: Insufficient or non-matching number of 

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

1207 

1208 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue. 

1209 

1210 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline} 

1211 warning as exception. 

1212 ''' 

1213 return _HeightBase._height(self, lats, lons) 

1214 

1215 

1216class HeightSmoothBiSpline(_HeightBase): 

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

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

1219 interpolate.SmoothSphereBivariateSpline.html>}. 

1220 ''' 

1221 _kmin = 16 # k = 3, always 

1222 

1223 def __init__(self, knots, s=4, name=NN): 

1224 '''New L{HeightSmoothBiSpline} interpolator. 

1225 

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

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

1228 @kwarg name: Optional name for this height interpolator (C{str}). 

1229 

1230 @raise HeightError: Insufficient number of B{C{knots}} or 

1231 an invalid B{C{knot}} or B{C{s}}. 

1232 

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

1234 or not installed. 

1235 

1236 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue. 

1237 

1238 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline} 

1239 warning as exception. 

1240 ''' 

1241 spi = self.scipy_interpolate 

1242 

1243 s = Float_(s, name=_smoothing_, Error=HeightError, low=4) 

1244 

1245 xs, ys, hs = self._xyhs3(knots) 

1246 try: 

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

1248 eps=EPS, s=s).ev 

1249 except Exception as x: 

1250 raise _SciPyIssue(x) 

1251 

1252 if name: 

1253 self.name = name 

1254 

1255 def __call__(self, *llis): 

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

1257 

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

1259 C{LatLon}s). 

1260 

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

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

1263 

1264 @raise HeightError: Insufficient number of B{C{llis}} or 

1265 an invalid B{C{lli}}. 

1266 

1267 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue. 

1268 

1269 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline} 

1270 warning as exception. 

1271 ''' 

1272 return _HeightBase._eval(self, llis) 

1273 

1274 def height(self, lats, lons): 

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

1276 

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

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

1279 

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

1281 interpolated heights (C{float}s). 

1282 

1283 @raise HeightError: Insufficient or non-matching number of 

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

1285 

1286 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue. 

1287 

1288 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline} 

1289 warning as exception. 

1290 ''' 

1291 return _HeightBase._height(self, lats, lons) 

1292 

1293 

1294__all__ += _ALL_DOCS(_HeightBase) 

1295 

1296# **) MIT License 

1297# 

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

1299# 

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

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

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

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

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

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

1306# 

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

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

1309# 

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

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

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

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

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

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

1316# OTHER DEALINGS IN THE SOFTWARE.