1

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

3

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

5

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

7 L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWcosineLaw},

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

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

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

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

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

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

14

15 Typical usage

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

17

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

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

20

21 C{>>> ...}

22

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

24

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

26

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

28 arguments to select different interpolation options

29

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

31

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

33

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

35

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

37

38 or

39

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

41

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

43

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

45

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

47

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

49

50 or as 2 lists, 2 tuples, etc.

51

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

53

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

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

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

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

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

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

60

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

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

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

64

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

66 Interpolation.

67'''

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

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

70

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

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

73 from pygeodesy.datums import _ellipsoidal_datum, _WGS84

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

75 _SciPyIssue, _xattr, _xkwds, _xkwds_get, _xkwds_item2

76# from pygeodesy.fmath import fidw # _MODS

77# from pygeodesy.formy import cosineAndoyerLambert, cosineForsytheAndoyerLambert, \

78# cosineLaw, equirectangular4, euclidean, flatLocal, \

79# flatPolar, haversine, thomas, vincentys # _MODS.into

80# from pygeodesy.internals import _version2 # _MODS

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

82 _PLUS_, _scipy_, _SPACE_, _STAR_

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

84 from pygeodesy.named import _name2__, _Named

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

86 from pygeodesy.props import Property_RO, property_RO, property_ROver

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

88 from pygeodesy.units import _isDegrees, Float_, Int_

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

90

91# from math import radians # from .points

92

93__all__ = _ALL_LAZY.heights

94__version__ = '24.11.02'

95

96_error_ = 'error'

97_formy = _MODS.into(formy=__name__)

98_linear_ = 'linear'

99_llis_ = 'llis'

100

101

102 class HeightError(PointsError):

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

104 '''

105 pass

106

107

108 def _alist(ais):

109 # return list of floats, not numpy.float64s

110 return list(map(float, ais))

111

112

113 def _ascalar(ais): # in .geoids

114 # return single float, not numpy.float64

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

116 if len(ais) != 1:

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

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

119 raise _AssertionError(n, txt=t)

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

121

122

123 def _atuple(ais):

124 # return tuple of floats, not numpy.float64s

125 return tuple(map(float, ais))

126

127

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

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

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

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

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

133

134 n = len(llis)

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

136 llis = llis[0]

137 try:

138 n, llis = len2(llis)

139 _as = _alist # return list of interpolated heights

140 except TypeError: # single lli

141 n, llis = 1, [llis]

142 _as = _ascalar # return single interpolated heights

143 else: # of 0, 2 or more llis

144 _as = _atuple # return tuple of interpolated heights

145

146 if n < m:

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

148 return _as, llis

149

150

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

152 # create an insufficient Error instance

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

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

155

156

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

158 # clip, order and remove duplicates

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

160

161

162 def _xyhs(wrap=False, _lat=_90_0, _lon=_180_0, **name_lls):

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

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

165 name, lls = _xkwds_item2(name_lls)

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

167 try:

168 for i, ll in enumerate(lls):

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

170 yield max(_0_0, _r(x + _lon)), \

171 max(_0_0, _r(y + _lat)), ll.height

172 except Exception as e:

173 i = Fmt.INDEX(name, i)

174 raise HeightError(i, ll, cause=e)

175

176

177 class _HeightNamed(_Named): # in .geoids

178 '''(INTERNAL) Interpolator base class.

179 '''

180 _datum = _WGS84 # default

181 _Error = HeightError

182 _kmin = 2 # min number of knots

183

184 _LLis = LatLon_ # ._height class

185 _np_sp = None # (numpy, scipy)

186 _wrap = None # wrap knots and llis

187

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

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

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

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

192 else:

193 n, lats = len2(lats)

194 m, lons = len2(lons)

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

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

197 raise self._Error(str(e))

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

199 return llis

200

201 @property_RO

202 def datum(self):

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

204 '''

205 return self._datum

206

207 @property_RO

208 def kmin(self):

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

210 '''

211 return self._kmin

212

213 @property_RO

214 def wrap(self):

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

216 '''

217 return self._wrap

218

219

220 class _HeightBase(_HeightNamed): # in .geoids

221 '''(INTERNAL) Interpolator base class.

222 '''

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

224 -2: _linear_, # for backward compatibility

225 -3: 'cubic',

226 -5: 'quintic'}

227

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

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

230

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

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

233 locations (C{bool}), overriding the B{C{knots}}

234 setting.

235

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

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

238

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

240 an invalid B{C{lli}}.

241

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

243

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

245 '''

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

247

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

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

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

251 atype = self.numpy.array

252 wrap = _xkwds(wrap, wrap=self._wrap)

253 _as, llis = _as_llis2(llis)

254 xis, yis, _ = zip(*_xyhs(llis=llis, **wrap)) # PYCHOK yield

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

256

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

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

259 self._notOverloaded(*args)

260

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

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

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

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

265 except Exception as x:

266 raise _SciPyIssue(x, self._ev_name)

267

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

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

270 self._notOverloaded(x, y)

271

272 @property_RO

273 def _ev_name(self):

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

275 '''

276 _ev = str(self._ev)

277 if _scipy_ not in _ev:

278 _ev = str(self._ev2d)

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

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

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

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

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

284

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

286 '''Interpolate the height for one or several lat-/longitudes. I{Must be overloaded}.

287

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

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

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

291 and B{C{lons}} locationts (C{bool}), overriding

292 the B{C{knots}} setting.

293

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

295 interpolated heights (C{float}s).

296

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

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

299

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

301

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

303 '''

304 lls = self._as_lls(lats, lons)

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

306

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

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

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

310 '''

311 try:

312 spi = self.scipy_interpolate

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

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

315 # DEPRECATED since scipy 1.10, removed altogether in 1.14

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

317

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

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

320 # spi.RectBivariateSpline needs strictly ordered xs and ys

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

322

323 def _bisplev(x, y):

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

325

326 self._ev2d = _bisplev

327

328 except Exception as x:

329 raise _SciPyIssue(x, self._ev_name)

330

331 def _kxky(self, kind):

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

333

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

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

336 '''

337 # raise SciPyWarnings, but not if

338 # scipy has already been imported

339 if throwarnings: # PYCHOK no cover

340 import sys

341 if _scipy_ not in sys.modules:

342 import warnings

343 warnings.filterwarnings(_error_)

344 return self.numpy, self.scipy

345

346 @property_ROver

347 def numpy(self):

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

349 '''

350 return _xnumpy(self.__class__, 1, 9) # overwrite property_ROver

351

352 @property_ROver

353 def scipy(self):

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

355 '''

356 return _xscipy(self.__class__, 1, 2) # overwrite property_ROver

357

358 @property_ROver

359 def scipy_interpolate(self):

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

361 '''

362 _ = self.scipy

363 import scipy.interpolate as spi # scipy 1.2.2

364 return spi # overwrite property_ROver

365

366 def _scipy_version(self, **n):

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

368 '''

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

370

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

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

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

374 n = len(hs)

375 if n < self.kmin:

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

377 if name:

378 self.name = name

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

380

381

382 class HeightCubic(_HeightBase):

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

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

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

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

387 '''

388 _kind = -3

389 _kmin = 16

390

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

392 '''New L{HeightCubic} interpolator.

393

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

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

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

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

398

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

400

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

402

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

404

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

406 '''

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

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

409

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

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

412

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

414 '''

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

416

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

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

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

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

421

422

423 class HeightLinear(HeightCubic):

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

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

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

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

428 '''

429 _kind = -1

430 _kmin = 2

431

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

433 '''New L{HeightLinear} interpolator.

434

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

436 '''

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

438

439 if _FOR_DOCS:

440 __call__ = HeightCubic.__call__

441 height = HeightCubic.height

442

443

444 class HeightLSQBiSpline(_HeightBase):

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

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

447 interpolate.LSQSphereBivariateSpline.html>}.

448 '''

449 _kmin = 16 # k = 3, always

450

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

452 '''New L{HeightLSQBiSpline} interpolator.

453

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

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

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

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

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

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

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

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

462

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

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

465

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

467

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

469 installed.

470

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

472

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

474 '''

475 np = self.numpy

476 spi = self.scipy_interpolate

477

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

479 n = len(hs)

480

481 w = weight

482 if isscalar(w):

483 w = float(w)

484 if w <= 0:

485 raise HeightError(weight=w)

486 w = (w,) * n

487 elif w is not None:

488 m, w = len2(w)

489 if m != n:

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

491 w = map2(float, w)

492 m = min(w)

493 if m <= 0: # PYCHOK no cover

494 w = Fmt.INDEX(weight=w.index(m))

495 raise HeightError(w, m)

496 try:

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

498 raise HeightError(low=low)

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

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

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

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

503 except Exception as x:

504 raise _SciPyIssue(x, self._ev_name)

505

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

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

508

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

510 '''

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

512

513

514 class HeightSmoothBiSpline(_HeightBase):

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

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

517 interpolate.SmoothSphereBivariateSpline.html>}.

518 '''

519 _kmin = 16 # k = 3, always

520

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

522 '''New L{HeightSmoothBiSpline} interpolator.

523

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

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

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

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

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

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

530

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

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

533

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

535 installed.

536

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

538

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

540 '''

541 spi = self.scipy_interpolate

542

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

544

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

546 try:

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

548 eps=EPS, s=s).ev

549 except Exception as x:

550 raise _SciPyIssue(x, self._ev_name)

551

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

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

554

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

556 '''

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

558

559

560 class _HeightIDW(_HeightNamed):

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

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

563 interpolators.

564

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

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

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

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

569 classes.

570 '''

571 _beta = 0 # fidw inverse power

572 _func = None # formy function

573 _knots = () # knots list or tuple

574 _kwds = {} # func_ options

575

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

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

578

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

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

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

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

583 retrievable with property C{kwds}.

584

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

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

587 '''

588 name, kwds = _name2__(**name__kwds)

589 if name:

590 self.name = name

591

592 n, self._knots = len2(knots)

593 if n < self.kmin:

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

595 self.beta = beta

596 self._kwds = kwds or {}

597

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

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

600

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

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

603 locations (C{bool}).

604

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

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

607

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

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

610 issue.

611 '''

612 def _xy2(wrap=False):

613 _w = _Wrap._latlonop(wrap)

614 try: # like _xyhs above, but degrees

615 for i, ll in enumerate(llis):

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

617 except Exception as x:

618 i = Fmt.INDEX(llis=i)

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

620

621 _as, llis = _as_llis2(llis)

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

623

624 @property_RO

625 def adjust(self):

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

627 '''

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

629

630 @property

631 def beta(self):

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

633 '''

634 return self._beta

635

636 @beta.setter # PYCHOK setter!

637 def beta(self, beta):

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

639

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

641 '''

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

643

644 @property_RO

645 def datum(self):

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

647 '''

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

649

650 def _datum_setter(self, datum):

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

652 '''

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

654 if d and d is not self._datum:

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

656

657 def _distances(self, x, y):

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

659 '''

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

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

662 self._notOverloaded(distance_function=_f)

663 try:

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

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

666 except Exception as e:

667 i = Fmt.INDEX(knots=i)

668 raise HeightError(i, k, cause=e)

669

670 def _distancesTo(self, _To):

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

672 '''

673 try:

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

675 yield _To(k)

676 except Exception as e:

677 i = Fmt.INDEX(knots=i)

678 raise HeightError(i, k, cause=e)

679

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

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

682

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

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

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

686 and B{C{lons}} (C{bool}).

687

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

689 interpolated heights (C{float}s).

690

691 @raise HeightError: Insufficient or non-matching number of B{C{lats}}

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

693 '''

694 lls = self._as_lls(lats, lons)

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

696

697 @Property_RO

698 def _heights(self):

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

700 '''

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

702

703 def _hIDW(self, x, y):

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

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

706 '''

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

708 try:

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

710 except (TypeError, ValueError) as e:

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

712

713 @property_RO

714 def hypot(self):

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

716 '''

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

718

719 @property_RO

720 def knots(self):

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

722 '''

723 return self._knots

724

725 @property_RO

726 def kwds(self):

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

728 '''

729 return self._kwds

730

731 @property_RO

732 def limit(self):

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

734 '''

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

736

737 @property_RO

738 def radius(self):

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

740 '''

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

742

743 @property_RO

744 def scaled(self):

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

746 '''

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

748

749 @property_RO

750 def wrap(self):

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

752 '''

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

754

755

756 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

760 L{pygeodesy.cosineAndoyerLambert_}.

761 '''

762 def __init__(self, knots, beta=2, **name__datum_wrap):

763 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

764

765 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height

766 interpolator (C{str}) and any keyword arguments for

767 function L{pygeodesy.cosineAndoyerLambert}.

768

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

770 '''

771 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)

772 self._func = _formy.cosineAndoyerLambert

773

774 if _FOR_DOCS:

775 __call__ = _HeightIDW.__call__

776 height = _HeightIDW.height

777

778

779 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

783 L{pygeodesy.cosineForsytheAndoyerLambert_}.

784 '''

785 def __init__(self, knots, beta=2, **name__datum_wrap):

786 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

787

788 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height

789 interpolator (C{str}) and any keyword arguments for

790 function L{pygeodesy.cosineForsytheAndoyerLambert}.

791

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

793 '''

794 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)

795 self._func = _formy.cosineForsytheAndoyerLambert

796

797 if _FOR_DOCS:

798 __call__ = _HeightIDW.__call__

799 height = _HeightIDW.height

800

801

802 class HeightIDWcosineLaw(_HeightIDW):

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

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

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

806

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

808 '''

809 def __init__(self, knots, beta=2, **name__radius_wrap):

810 '''New L{HeightIDWcosineLaw} interpolator.

811

812 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this height

813 interpolator (C{str}) and any keyword arguments for

814 function L{pygeodesy.cosineLaw}.

815

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

817 '''

818 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)

819 self._func = _formy.cosineLaw

820

821 if _FOR_DOCS:

822 __call__ = _HeightIDW.__call__

823 height = _HeightIDW.height

824

825

826 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

830 conventionally in C{meter}.

831 '''

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

833 '''New L{HeightIDWdistanceTo} interpolator.

834

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

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

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

838

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

840

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

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

843 checked}.

844 '''

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

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

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

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

849 self._datum = ks0.datum

850 self._LLis = ks0.classof # type(ks0)

851

852 def _distances(self, x, y):

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

854 '''

855 kwds, ll = self._kwds, self._LLis(y, x)

856

857 def _To(k):

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

859

860 return self._distancesTo(_To)

861

862 if _FOR_DOCS:

863 __call__ = _HeightIDW.__call__

864 height = _HeightIDW.height

865

866

867 class HeightIDWequirectangular(_HeightIDW):

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

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

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

871 L{pygeodesy.equirectangular4}.

872 '''

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

874 '''New L{HeightIDWequirectangular} interpolator.

875

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

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

878 for function L{pygeodesy.equirectangular4}.

879

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

881 '''

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

883

884 def _distances(self, x, y):

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

886 '''

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

888

889 def _To(k):

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

891

892 return self._distancesTo(_To)

893

894 if _FOR_DOCS:

895 __call__ = _HeightIDW.__call__

896 height = _HeightIDW.height

897

898

899 class HeightIDWeuclidean(_HeightIDW):

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

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

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

903 '''

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

905 '''New L{HeightIDWeuclidean} interpolator.

906

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

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

909 for function function L{pygeodesy.euclidean}.

910

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

912 '''

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

914 self._func = _formy.euclidean

915

916 if _FOR_DOCS:

917 __call__ = _HeightIDW.__call__

918 height = _HeightIDW.height

919

920

921 class HeightIDWexact(_HeightIDW):

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

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

924 and the I{angular} distance in C{degrees} from method

925 L{GeodesicExact.Inverse}.

926 '''

927 def __init__(self, knots, beta=2, datum=None, **name__wrap):

928 '''New L{HeightIDWexact} interpolator.

929

930 @kwarg datum: Datum to override the default C{Datums.WGS84} and

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

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

933 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator

934 (C{str}) and a keyword argument for method C{Inverse1} of

935 class L{geodesicx.GeodesicExact}.

936

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

938

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

940 '''

941 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)

942 self._datum_setter(datum)

943 self._func = self.datum.ellipsoid.geodesicx.Inverse1

944

945 if _FOR_DOCS:

946 __call__ = _HeightIDW.__call__

947 height = _HeightIDW.height

948

949

950 class HeightIDWflatLocal(_HeightIDW):

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

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

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

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

955 '''

956 def __init__(self, knots, beta=2, **name__datum_hypot_scaled_wrap):

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

958

959 @kwarg name__datum_hypot_scaled_wrap: Optional C{B{name}=NN}

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

961 keyword arguments for L{pygeodesy.flatLocal}.

962

963 @see: L{HeightIDW<_HeightIDW.__init__>} for further details.

964 '''

965 _HeightIDW.__init__(self, knots, beta=beta,

966 **name__datum_hypot_scaled_wrap)

967 self._func = _formy.flatLocal

968

969 if _FOR_DOCS:

970 __call__ = _HeightIDW.__call__

971 height = _HeightIDW.height

972

973

974 class HeightIDWflatPolar(_HeightIDW):

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

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

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

978 L{pygeodesy.flatPolar_}.

979 '''

980 def __init__(self, knots, beta=2, **name__radius_wrap):

981 '''New L{HeightIDWflatPolar} interpolator.

982

983 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this

984 height interpolator (C{str}) and any keyword

985 arguments for function L{pygeodesy.flatPolar}.

986

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

988 '''

989 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)

990 self._func = _formy.flatPolar

991

992 if _FOR_DOCS:

993 __call__ = _HeightIDW.__call__

994 height = _HeightIDW.height

995

996

997 class HeightIDWhaversine(_HeightIDW):

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

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

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

1001

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

1003 '''

1004 def __init__(self, knots, beta=2, **name__radius_wrap):

1005 '''New L{HeightIDWhaversine} interpolator.

1006

1007 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this

1008 height interpolator (C{str}) and any keyword

1009 arguments for function L{pygeodesy.haversine}.

1010

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

1012 '''

1013 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)

1014 self._func = _formy.haversine

1015

1016 if _FOR_DOCS:

1017 __call__ = _HeightIDW.__call__

1018 height = _HeightIDW.height

1019

1020

1021 class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny

1022 if _FOR_DOCS:

1023 __doc__ = HeightIDWflatLocal.__doc__

1024 __init__ = HeightIDWflatLocal.__init__

1025 __call__ = HeightIDWflatLocal.__call__

1026 height = HeightIDWflatLocal.height

1027

1028

1029 class HeightIDWkarney(_HeightIDW):

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

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

1032 I{angular} distance in C{degrees} from I{Karney}'s U{geographiclib

1033 <https://PyPI.org/project/geographiclib>} method U{Geodesic.Inverse

1034 <https://GeographicLib.SourceForge.io/Python/doc/code.html#

1035 geographiclib.geodesic.Geodesic.Inverse>}.

1036 '''

1037 def __init__(self, knots, beta=2, datum=None, **name__wrap):

1038 '''New L{HeightIDWkarney} interpolator.

1039

1040 @kwarg datum: Datum to override the default C{Datums.WGS84} and

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

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

1043 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator

1044 (C{str}) and a keyword argument for method C{Inverse1} of

1045 class L{geodesicw.Geodesic}.

1046

1047 @raise ImportError: Package U{geographiclib

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

1049

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

1051

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

1053 '''

1054 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)

1055 self._datum_setter(datum)

1056 self._func = self.datum.ellipsoid.geodesic.Inverse1

1057

1058 if _FOR_DOCS:

1059 __call__ = _HeightIDW.__call__

1060 height = _HeightIDW.height

1061

1062

1063 class HeightIDWthomas(_HeightIDW):

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

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

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

1067 '''

1068 def __init__(self, knots, beta=2, **name__datum_wrap):

1069 '''New L{HeightIDWthomas} interpolator.

1070

1071 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this

1072 height interpolator (C{str}) and any keyword

1073 arguments for function L{pygeodesy.thomas}.

1074

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

1076 '''

1077 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)

1078 self._func = _formy.thomas

1079

1080 if _FOR_DOCS:

1081 __call__ = _HeightIDW.__call__

1082 height = _HeightIDW.height

1083

1084

1085 class HeightIDWvincentys(_HeightIDW):

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

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

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

1089

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

1091 '''

1092 def __init__(self, knots, beta=2, **name__radius_wrap):

1093 '''New L{HeightIDWvincentys} interpolator.

1094

1095 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this

1096 height interpolator (C{str}) and any keyword

1097 arguments for function L{pygeodesy.vincentys}.

1098

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

1100 '''

1101 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)

1102 self._func = _formy.vincentys

1103

1104 if _FOR_DOCS:

1105 __call__ = _HeightIDW.__call__

1106 height = _HeightIDW.height

1107

1108

1109__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightNamed)

1110

1111# **) MIT License

1112#

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

1114#

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

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

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

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

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

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

1121#

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

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

1124#

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

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

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

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

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

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

1131# OTHER DEALINGS IN THE SOFTWARE.