1

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

3

4 u'''Height interpolations at C{LatLon} points from 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, 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

76# from pygeodesy.fmath import fidw # _MODS

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

78# cosineLaw, equirectangular_, euclidean, flatLocal, \

79# flatPolar, haversine, thomas, vincentys # _MODS

80 from pygeodesy.interns import NN, _COMMASPACE_, _cubic_, _insufficient_, _knots_, \

81 _linear_, _NOTEQUAL_, _PLUS_, _scipy_, _SPACE_, \

82 _STAR_, _version2

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

84 from pygeodesy.named import _Named

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

86 from pygeodesy.props import Property_RO, property_RO

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

95

96_error_ = 'error'

97_llis_ = 'llis'

98_smoothing_ = 'smoothing'

99

100

101 class HeightError(PointsError):

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

103 '''

104 pass

105

106

107 def _alist(ais):

108 # return list of floats, not numpy.float64s

109 return list(map(float, ais))

110

111

112 def _ascalar(ais): # in .geoids

113 # return single float, not numpy.float64

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

115 if len(ais) != 1:

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

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

118 raise _AssertionError(n, txt=t)

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

120

121

122 def _atuple(ais):

123 # return tuple of floats, not numpy.float64s

124 return tuple(map(float, ais))

125

126

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

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

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

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

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

132

133 n = len(llis)

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

135 llis = llis[0]

136 try:

137 n, llis = len2(llis)

138 _as = _alist # return list of interpolated heights

139 except TypeError: # single lli

140 n, llis = 1, [llis]

141 _as = _ascalar # return single interpolated heights

142 else: # of 0, 2 or more llis

143 _as = _atuple # return tuple of interpolated heights

144

145 if n < m:

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

147 return _as, llis

148

149

150 def _height_called(inst, lats, lons, Error=HeightError, **wrap): # in .geoids

151 LLis, d = inst._LLis, inst.datum

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

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

154 else:

155 n, lats = len2(lats)

156 m, lons = len2(lons)

157 if n != m:

158 # format a LenError, but raise an Error

159 e = LenError(inst.__class__, lats=n, lons=m, txt=None)

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

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

162 return inst(llis, **wrap) # __call__(lli) or __call__(llis)

163

164

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

166 # create an insufficient Error instance

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

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

169

170

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

172 # clip, order and remove duplicates

173 # p, ks = 0, []

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

175 # if k > p:

176 # ks.append(k)

177 # p = k

178 # return ks

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

180

181

182 def _xyhs(lls, wrap=False, name=_llis_):

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

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

185 _0, _90, _180 = _0_0, _90_0, _180_0

186 _m, _r, _w = max, radians, _Wrap._latlonop(wrap)

187 try:

188 for i, ll in enumerate(lls):

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

190 yield _m(_0, _r(x + _180)), \

191 _m(_0, _r(y + _90)), ll.height

192 except Exception as e:

193 i = Fmt.INDEX(name, i)

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

195

196

197 class _HeightBase(_Named): # in .geoids

198 '''(INTERNAL) Interpolator base class.

199 '''

200 _datum = _WGS84 # default

201 _kmin = 2 # min number of knots

202 _LLis = LatLon_ # ._height class

203 _np_sp = None # (numpy, scipy)

204 _wrap = None # wrap knots and llis

205

206 @property_RO

207 def datum(self):

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

209 '''

210 return self._datum

211

212 @property_RO

213 def kmin(self):

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

215 '''

216 return self._kmin

217

218 @property_RO

219 def wrap(self):

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

221 '''

222 return self._wrap

223

224

225 class _HeightsBase(_HeightBase): # in .geoids

226 '''(INTERNAL) Interpolator base class.

227 '''

228 _np_sp = None # (numpy, scipy)

229

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

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

232

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

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

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

236 setting.

237

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

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

240

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

242 an invalid B{C{lli}}.

243

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

245

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

247 '''

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

249

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

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

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

253 atype = self.numpy.array

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

255 _as, llis = _as_llis2(llis)

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

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

258

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

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

261 self._notOverloaded(*args)

262

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

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

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

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

267 except Exception as x:

268 raise _SciPyIssue(x)

269

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

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

272

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

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

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

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

277 the B{C{knots}} setting.

278

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

280 interpolated heights (C{float}s).

281

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

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

284

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

286

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

288 '''

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

290

291 def _np_sp2(self, throwarnings=False):

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

293 '''

294 t = _HeightsBase._np_sp

295 if not t:

296 # raise SciPyWarnings, but not if

297 # scipy has already been imported

298 if throwarnings: # PYCHOK no cover

299 import sys

300 if _scipy_ not in sys.modules:

301 import warnings

302 warnings.filterwarnings(_error_)

303

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

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

306

307 _HeightsBase._np_sp = t = np, sp

308 return t

309

310 @Property_RO

311 def numpy(self):

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

313 '''

314 np, _ = self._np_sp2()

315 return np

316

317 @Property_RO

318 def scipy(self):

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

320 '''

321 _, sp = self._np_sp2()

322 return sp

323

324 @Property_RO

325 def scipy_interpolate(self):

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

327 '''

328 _ = self.scipy

329 import scipy.interpolate as spi # scipy 1.2.2

330 return spi

331

332 def _scipy_version(self, **n):

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

334 '''

335 return _version2(self.scipy.version.version, **n)

336

337 def _xyhs3(self, knots, **wrap):

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

339 xs, ys, hs = zip(*_xyhs(knots, name=_knots_, **wrap)) # PYCHOK yield

340 n = len(hs)

341 if n < self.kmin:

342 raise _insufficientError(self.kmin, knots=n)

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

344

345

346 class HeightCubic(_HeightsBase):

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

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

349 C{kind='cubic'}.

350 '''

351 _interp2d = None

352 _kind = _cubic_

353 _kmin = 16

354

355 def __init__(self, knots, name=NN, **wrap):

356 '''New L{HeightCubic} interpolator.

357

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

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

360 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{knots}}

361 and B{C{llis}} locations (C{bool}).

362

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

364 invalid B{C{knot}}.

365

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

367 or not installed.

368

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

370

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

372 as exception.

373 '''

374 spi = self.scipy_interpolate

375

376 xs, ys, hs = self._xyhs3(knots, **wrap)

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

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

379 except Exception as x:

380 raise _SciPyIssue(x)

381

382 if name:

383 self.name = name

384

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

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

387

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

389

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

391 as exception.

392

393 @see: L{Here<_HeightsBase.__call__>} for other details.

394 '''

395 return _HeightsBase._eval(self, llis, **wrap)

396

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

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

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

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

401

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

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

404

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

406

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

408 as exception.

409

410 @see: L{Here<_HeightsBase.height>} for other details.

411 '''

412 return _height_called(self, lats, lons, **wrap)

413

414

415 class HeightLinear(HeightCubic):

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

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

418 C{kind='linear'}.

419 '''

420 _kind = _linear_

421 _kmin = 2

422

423 def __init__(self, knots, name=NN, **wrap):

424 '''New L{HeightLinear} interpolator.

425

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

427 '''

428 HeightCubic.__init__(self, knots, name=name, **wrap)

429

430 if _FOR_DOCS:

431 __call__ = HeightCubic.__call__

432 height = HeightCubic.height

433

434

435 class HeightLSQBiSpline(_HeightsBase):

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

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

438 interpolate.LSQSphereBivariateSpline.html>}.

439 '''

440 _kmin = 16 # k = 3, always

441

442 def __init__(self, knots, weight=None, name=NN, **wrap):

443 '''New L{HeightLSQBiSpline} interpolator.

444

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

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

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

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

449 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{knots}}

450 and B{C{llis}} locations (C{bool}).

451

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

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

454

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

456

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

458 or not installed.

459

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

461

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

463 warning as exception.

464 '''

465 np = self.numpy

466 spi = self.scipy_interpolate

467

468 xs, ys, hs = self._xyhs3(knots, **wrap)

469 n = len(hs)

470

471 w = weight

472 if isscalar(w):

473 w = float(w)

474 if w <= 0:

475 raise HeightError(weight=w)

476 w = [w] * n

477 elif w is not None:

478 m, w = len2(w)

479 if m != n:

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

481 w = map2(float, w)

482 m = min(w)

483 if m <= 0: # PYCHOK no cover

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

485 raise HeightError(w, m)

486 try:

487 T = 1.0e-4 # like SciPy example

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

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

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

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

492 except Exception as x:

493 raise _SciPyIssue(x)

494

495 if name:

496 self.name = name

497

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

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

500

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

502

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

504 warning as exception.

505

506 @see: L{Here<_HeightsBase.__call__>} for other details.

507 '''

508 return _HeightsBase._eval(self, llis, **wrap)

509

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

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

512

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

514

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

516 warning as exception.

517

518 @see: L{Here<_HeightsBase.height>} for other details.

519 '''

520 return _height_called(self, lats, lons, **wrap)

521

522

523 class HeightSmoothBiSpline(_HeightsBase):

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

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

526 interpolate.SmoothSphereBivariateSpline.html>}.

527 '''

528 _kmin = 16 # k = 3, always

529

530 def __init__(self, knots, s=4, name=NN, **wrap):

531 '''New L{HeightSmoothBiSpline} interpolator.

532

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

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

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

536 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{knots}}

537 and any called B{C{llis}} and height B{C{lats}}

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

539

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

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

542

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

544 or not installed.

545

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

547

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

549 warning as exception.

550 '''

551 spi = self.scipy_interpolate

552

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

554

555 xs, ys, hs = self._xyhs3(knots, **wrap)

556 try:

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

558 eps=EPS, s=s).ev

559 except Exception as x:

560 raise _SciPyIssue(x)

561

562 if name:

563 self.name = name

564

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

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

567

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

569

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

571 warning as exception.

572

573 @see: L{Here<_HeightsBase.__call__>} for other details.

574 '''

575 return _HeightsBase._eval(self, llis, **wrap)

576

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

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

579

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

581

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

583 warning as exception.

584

585 @see: L{Here<_HeightsBase.height>} for other details.

586 '''

587 return _height_called(self, lats, lons, **wrap)

588

589

590 class _HeightIDW(_HeightBase):

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

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

593 interpolators.

594

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

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

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

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

599 classes.

600 '''

601 _beta = 0 # fidw inverse power

602 _func = None # formy function

603 _knots = () # knots list or tuple

604 _kwds = {} # func_ options

605

606 def __init__(self, knots, beta=2, name=NN, **kwds):

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

608

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

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

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

612 @kwarg kwds: Optional keyword argument for distance function,

613 retrievable with property C{kwds}.

614

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

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

617 '''

618 if name:

619 self.name = name

620 n, self._knots = len2(knots)

621 if n < self.kmin:

622 raise _insufficientError(self.kmin, knots=n)

623 self.beta = beta

624 self._kwds = kwds or {}

625

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

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

628

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

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

631 locations (C{bool}).

632

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

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

635

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

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

638 issue.

639 '''

640 def _xy2(lls, wrap=False):

641 _w = _Wrap._latlonop(wrap)

642 try: # like _xyhs above, but degrees

643 for i, ll in enumerate(lls):

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

645 except Exception as e:

646 i = Fmt.INDEX(llis=i)

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

648

649 _as, llis = _as_llis2(llis)

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

651

652 @property_RO

653 def adjust(self):

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

655 '''

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

657

658 @property

659 def beta(self):

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

661 '''

662 return self._beta

663

664 @beta.setter # PYCHOK setter!

665 def beta(self, beta):

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

667

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

669 '''

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

671

672 @property_RO

673 def datum(self):

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

675 '''

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

677

678 def _datum_setter(self, datum):

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

680 '''

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

682 if d and d is not self._datum:

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

684

685 def _distances(self, x, y):

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

687 '''

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

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

690 self._notOverloaded(distance_function=_f)

691 try:

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

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

694 except Exception as e:

695 i = Fmt.INDEX(knots=i)

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

697

698 @property_RO

699 def _fmath(self):

700 '''(INTERNAL) Get module C{fmath}, I{once}.

701 '''

702 _HeightIDW._fmath = f = _MODS.fmath # overwrite property_RO

703 return f

704

705 @property_RO

706 def _formy(self):

707 '''(INTERNAL) Get module C{formy}, I{once}.

708 '''

709 _HeightIDW._formy = f = _MODS.formy # overwrite property_RO

710 return f

711

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

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

714

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

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

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

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

719

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

721 interpolated heights (C{float}s).

722

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

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

725 '''

726 return _height_called(self, lats, lons, **wrap)

727

728 @Property_RO

729 def _heights(self):

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

731 '''

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

733

734 def _hIDW(self, x, y):

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

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

737 '''

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

739 try:

740 return self._fmath.fidw(hs, ds, beta=self.beta)

741 except (TypeError, ValueError) as e:

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

743

744 @property_RO

745 def hypot(self):

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

747 '''

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

749

750 @property_RO

751 def knots(self):

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

753 '''

754 return self._knots

755

756 @property_RO

757 def kwds(self):

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

759 '''

760 return self._kwds

761

762 @property_RO

763 def limit(self):

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

765 '''

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

767

768 @property_RO

769 def radius(self):

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

771 '''

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

773

774 @property_RO

775 def scaled(self):

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

777 '''

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

779

780 @property_RO

781 def wrap(self):

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

783 '''

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

785

786

787 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

791 L{pygeodesy.cosineAndoyerLambert_}.

792 '''

793 def __init__(self, knots, beta=2, name=NN, **datum_wrap):

794 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

795

796 @kwarg datum_wrap: Optional keyword arguments for function

797 L{pygeodesy.cosineAndoyerLambert}.

798

799 @see: L{Here<_HeightIDW.__init__>} for other details.

800 '''

801 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)

802 self._func = self._formy.cosineAndoyerLambert

803

804 if _FOR_DOCS:

805 __call__ = _HeightIDW.__call__

806 height = _HeightIDW.height

807

808

809 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

813 L{pygeodesy.cosineForsytheAndoyerLambert_}.

814 '''

815 def __init__(self, knots, beta=2, name=NN, **datum_wrap):

816 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

817

818 @kwarg datum_wrap: Optional keyword arguments for function

819 L{pygeodesy.cosineForsytheAndoyerLambert}.

820

821 @see: L{Here<_HeightIDW.__init__>} for other details.

822 '''

823 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)

824 self._func = self._formy.cosineForsytheAndoyerLambert

825

826 if _FOR_DOCS:

827 __call__ = _HeightIDW.__call__

828 height = _HeightIDW.height

829

830

831 class HeightIDWcosineLaw(_HeightIDW):

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

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

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

835

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

837 '''

838 def __init__(self, knots, beta=2, name=NN, **radius_wrap):

839 '''New L{HeightIDWcosineLaw} interpolator.

840

841 @kwarg radius_wrap: Optional keyword arguments for function

842 L{pygeodesy.cosineLaw}.

843

844 @see: L{Here<_HeightIDW.__init__>} for other details.

845 '''

846 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)

847 self._func = self._formy.cosineLaw

848

849 if _FOR_DOCS:

850 __call__ = _HeightIDW.__call__

851 height = _HeightIDW.height

852

853

854 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

858 conventionally in C{meter}.

859 '''

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

861 '''New L{HeightIDWdistanceTo} interpolator.

862

863 @kwarg distanceTo_kwds: Optional keyword arguments for the

864 B{C{knots}}' C{LatLon.distanceTo}

865 method.

866

867 @see: L{Here<_HeightIDW.__init__>} for other details.

868

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

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

871 checked}.

872 '''

873 _HeightIDW.__init__(self, knots, beta=beta, name=name,

874 **distanceTo_kwds)

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

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

877 # in _height_called instead of class LatLon_ and datum None

878 self._datum = ks0.datum

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

880

881 def _distances(self, x, y):

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

883 '''

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

885 try:

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

887 yield k.distanceTo(ll, **kwds)

888 except Exception as e:

889 i = Fmt.INDEX(knots=i)

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

891

892 if _FOR_DOCS:

893 __call__ = _HeightIDW.__call__

894 height = _HeightIDW.height

895

896

897 class HeightIDWequirectangular(_HeightIDW):

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

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

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

901 L{pygeodesy.equirectangular_}.

902 '''

903 def __init__(self, knots, beta=2, name=NN, **adjust_limit_wrap): # XXX beta=1

904 '''New L{HeightIDWequirectangular} interpolator.

905

906 @kwarg adjust_limit_wrap: Optional keyword arguments for

907 function L{pygeodesy.equirectangular_}.

908

909 @see: L{Here<_HeightIDW.__init__>} for other details.

910 '''

911 _HeightIDW.__init__(self, knots, beta=beta, name=name,

912 **adjust_limit_wrap)

913

914 def _distances(self, x, y):

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

916 '''

917 _f, kwds = self._formy.equirectangular_, self._kwds

918 try:

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

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

921 except Exception as e:

922 i = Fmt.INDEX(knots=i)

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

924

925 if _FOR_DOCS:

926 __call__ = _HeightIDW.__call__

927 height = _HeightIDW.height

928

929

930 class HeightIDWeuclidean(_HeightIDW):

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

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

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

934 '''

935 def __init__(self, knots, beta=2, name=NN, **adjust_radius_wrap):

936 '''New L{HeightIDWeuclidean} interpolator.

937

938 @kwarg adjust_radius_wrap: Optional keyword arguments for

939 function L{pygeodesy.euclidean}.

940

941 @see: L{Here<_HeightIDW.__init__>} for other details.

942 '''

943 _HeightIDW.__init__(self, knots, beta=beta, name=name,

944 **adjust_radius_wrap)

945 self._func = self._formy.euclidean

946

947 if _FOR_DOCS:

948 __call__ = _HeightIDW.__call__

949 height = _HeightIDW.height

950

951

952 class HeightIDWexact(_HeightIDW):

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

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

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

956 L{GeodesicExact.Inverse}.

957 '''

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

959 '''New L{HeightIDWexact} interpolator.

960

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

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

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

964 @kwarg wrap: Optional keyword argument for method C{Inverse1}

965 of class L{geodesicx.GeodesicExact}.

966

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

968

969 @see: L{Here<_HeightIDW.__init__>} for other details.

970 '''

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

972 self._datum_setter(datum)

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

974

975 if _FOR_DOCS:

976 __call__ = _HeightIDW.__call__

977 height = _HeightIDW.height

978

979

980 class HeightIDWflatLocal(_HeightIDW):

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

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

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

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

985 '''

986 def __init__(self, knots, beta=2, name=NN, **datum_hypot_scaled_wrap):

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

988

989 @kwarg datum_hypot_scaled_wrap: Optional keyword arguments for

990 function L{pygeodesy.flatLocal}.

991

992 @see: L{HeightIDW<_HeightIDW.__init__>} for other details.

993 '''

994 _HeightIDW.__init__(self, knots, beta=beta, name=name,

995 **datum_hypot_scaled_wrap)

996 self._func = self._formy.flatLocal

997

998 if _FOR_DOCS:

999 __call__ = _HeightIDW.__call__

1000 height = _HeightIDW.height

1001

1002

1003 class HeightIDWflatPolar(_HeightIDW):

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

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

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

1007 L{pygeodesy.flatPolar_}.

1008 '''

1009 def __init__(self, knots, beta=2, name=NN, **radius_wrap):

1010 '''New L{HeightIDWflatPolar} interpolator.

1011

1012 @kwarg radius_wrap: Optional keyword arguments for function

1013 L{pygeodesy.flatPolar}.

1014

1015 @see: L{Here<_HeightIDW.__init__>} for other details.

1016 '''

1017 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)

1018 self._func = self._formy.flatPolar

1019

1020 if _FOR_DOCS:

1021 __call__ = _HeightIDW.__call__

1022 height = _HeightIDW.height

1023

1024

1025 class HeightIDWhaversine(_HeightIDW):

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

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

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

1029

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

1031 '''

1032 def __init__(self, knots, beta=2, name=NN, **radius_wrap):

1033 '''New L{HeightIDWhaversine} interpolator.

1034

1035 @kwarg radius_wrap: Optional keyword arguments for function

1036 L{pygeodesy.haversine}.

1037

1038 @see: L{Here<_HeightIDW.__init__>} for other details.

1039 '''

1040 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)

1041 self._func = self._formy.haversine

1042

1043 if _FOR_DOCS:

1044 __call__ = _HeightIDW.__call__

1045 height = _HeightIDW.height

1046

1047

1048 class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny

1049 if _FOR_DOCS:

1050 __doc__ = HeightIDWflatLocal.__doc__

1051 __init__ = HeightIDWflatLocal.__init__

1052 __call__ = HeightIDWflatLocal.__call__

1053 height = HeightIDWflatLocal.height

1054

1055

1056 class HeightIDWkarney(_HeightIDW):

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

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

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

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

1061 <https://geographiclib.sourceforge.io/Python/doc/code.html#

1062 geographiclib.geodesic.Geodesic.Inverse>}.

1063 '''

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

1065 '''New L{HeightIDWkarney} interpolator.

1066

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

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

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

1070 @kwarg wrap: Optional keyword argument for method C{Inverse1}

1071 of class L{geodesicw.Geodesic}.

1072

1073 @raise ImportError: Package U{geographiclib

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

1075

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

1077

1078 @see: L{Here<_HeightIDW.__init__>} for other details.

1079 '''

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

1081 self._datum_setter(datum)

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

1083

1084 if _FOR_DOCS:

1085 __call__ = _HeightIDW.__call__

1086 height = _HeightIDW.height

1087

1088

1089 class HeightIDWthomas(_HeightIDW):

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

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

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

1093 '''

1094 def __init__(self, knots, beta=2, name=NN, **datum_wrap):

1095 '''New L{HeightIDWthomas} interpolator.

1096

1097 @kwarg datum_wrap: Optional keyword argument for function

1098 L{pygeodesy.thomas}.

1099

1100 @see: L{Here<_HeightIDW.__init__>} for other details.

1101 '''

1102 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)

1103 self._func = self._formy.thomas

1104

1105 if _FOR_DOCS:

1106 __call__ = _HeightIDW.__call__

1107 height = _HeightIDW.height

1108

1109

1110 class HeightIDWvincentys(_HeightIDW):

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

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

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

1114

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

1116 '''

1117 def __init__(self, knots, beta=2, name=NN, **radius_wrap):

1118 '''New L{HeightIDWvincentys} interpolator.

1119

1120 @kwarg radius_wrap: Optional keyword arguments for function

1121 L{pygeodesy.vincentys}.

1122

1123 @see: L{Here<_HeightIDW.__init__>} for other details.

1124 '''

1125 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)

1126 self._func = self._formy.vincentys

1127

1128 if _FOR_DOCS:

1129 __call__ = _HeightIDW.__call__

1130 height = _HeightIDW.height

1131

1132

1133__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightsBase)

1134

1135# **) MIT License

1136#

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

1138#

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

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

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

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

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

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

1145#

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

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

1148#

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

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

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

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

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

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

1155# OTHER DEALINGS IN THE SOFTWARE.