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, _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

80# from pygeodesy.internals import _version2 # _MODS

81 from pygeodesy.interns import NN, _COMMASPACE_, _cubic_, _insufficient_, _linear_, \

82 _NOTEQUAL_, _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

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

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(wrap=False, **name_lls):

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 name, lls = _xkwds_item2(name_lls)

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

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

188 try:

189 for i, ll in enumerate(lls):

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

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

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

193 except Exception as e:

194 i = Fmt.INDEX(name, i)

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

196

197

198 class _HeightBase(_Named): # in .geoids

199 '''(INTERNAL) Interpolator base class.

200 '''

201 _datum = _WGS84 # default

202 _kmin = 2 # min number of knots

203 _LLis = LatLon_ # ._height class

204 _np_sp = None # (numpy, scipy)

205 _wrap = None # wrap knots and llis

206

207 @property_RO

208 def datum(self):

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

210 '''

211 return self._datum

212

213 @property_RO

214 def kmin(self):

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

216 '''

217 return self._kmin

218

219 @property_RO

220 def wrap(self):

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

222 '''

223 return self._wrap

224

225

226 class _HeightsBase(_HeightBase): # in .geoids

227 '''(INTERNAL) Interpolator base class.

228 '''

229 _np_sp = None # (numpy, scipy)

230

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

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

233

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

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

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

237 setting.

238

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

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

241

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

243 an invalid B{C{lli}}.

244

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

246

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

248 '''

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

250

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

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

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

254 atype = self.numpy.array

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

256 _as, llis = _as_llis2(llis)

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

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

259

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

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

262 self._notOverloaded(*args)

263

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

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

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

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

268 except Exception as x:

269 raise _SciPyIssue(x)

270

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

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

273

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

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

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

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

278 the B{C{knots}} setting.

279

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

281 interpolated heights (C{float}s).

282

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

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

285

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

287

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

289 '''

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

291

292 def _np_sp2(self, throwarnings=False):

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

294 '''

295 t = _HeightsBase._np_sp

296 if not t:

297 # raise SciPyWarnings, but not if

298 # scipy has already been imported

299 if throwarnings: # PYCHOK no cover

300 import sys

301 if _scipy_ not in sys.modules:

302 import warnings

303 warnings.filterwarnings(_error_)

304

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

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

307

308 _HeightsBase._np_sp = t = np, sp

309 return t

310

311 @Property_RO

312 def numpy(self):

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

314 '''

315 np, _ = self._np_sp2()

316 return np

317

318 @Property_RO

319 def scipy(self):

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

321 '''

322 _, sp = self._np_sp2()

323 return sp

324

325 @Property_RO

326 def scipy_interpolate(self):

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

328 '''

329 _ = self.scipy

330 import scipy.interpolate as spi # scipy 1.2.2

331 return spi

332

333 def _scipy_version(self, **n):

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

335 '''

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

337

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

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

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

341 n = len(hs)

342 if n < self.kmin:

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

344 if name:

345 self.name = name

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

347

348

349 class HeightCubic(_HeightsBase):

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

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

352 C{kind='cubic'}.

353 '''

354 _interp2d = None

355 _kind = _cubic_

356 _kmin = 16

357

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

359 '''New L{HeightCubic} interpolator.

360

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

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

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

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

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

366

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

368 invalid B{C{knot}}.

369

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

371 or not installed.

372

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

374

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

376 as exception.

377 '''

378 spi = self.scipy_interpolate

379

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

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

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

383 except Exception as x:

384 raise _SciPyIssue(x)

385

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

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

388

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

390

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

392 as exception.

393

394 @see: L{Here<_HeightsBase.__call__>} for further details.

395 '''

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

397

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

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

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

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

402

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

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

405

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

407

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

409 as exception.

410

411 @see: L{Here<_HeightsBase.height>} for further details.

412 '''

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

414

415

416 class HeightLinear(HeightCubic):

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

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

419 C{kind='linear'}.

420 '''

421 _kind = _linear_

422 _kmin = 2

423

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

425 '''New L{HeightLinear} interpolator.

426

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

428 '''

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

430

431 if _FOR_DOCS:

432 __call__ = HeightCubic.__call__

433 height = HeightCubic.height

434

435

436 class HeightLSQBiSpline(_HeightsBase):

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

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

439 interpolate.LSQSphereBivariateSpline.html>}.

440 '''

441 _kmin = 16 # k = 3, always

442

443 def __init__(self, knots, weight=None, **name_wrap):

444 '''New L{HeightLSQBiSpline} interpolator.

445

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

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

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

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

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

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

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

453

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

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

456

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

458

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

460 not installed.

461

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

463

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

465 warning as exception.

466 '''

467 np = self.numpy

468 spi = self.scipy_interpolate

469

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

471 n = len(hs)

472

473 w = weight

474 if isscalar(w):

475 w = float(w)

476 if w <= 0:

477 raise HeightError(weight=w)

478 w = [w] * n

479 elif w is not None:

480 m, w = len2(w)

481 if m != n:

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

483 w = map2(float, w)

484 m = min(w)

485 if m <= 0: # PYCHOK no cover

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

487 raise HeightError(w, m)

488 try:

489 T = 1.0e-4 # like SciPy example

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

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

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

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

494 except Exception as x:

495 raise _SciPyIssue(x)

496

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

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

499

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

501

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

503 warning as exception.

504

505 @see: L{Here<_HeightsBase.__call__>} for further details.

506 '''

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

508

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

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

511

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

513

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

515 warning as exception.

516

517 @see: L{Here<_HeightsBase.height>} for further details.

518 '''

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

520

521

522 class HeightSmoothBiSpline(_HeightsBase):

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

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

525 interpolate.SmoothSphereBivariateSpline.html>}.

526 '''

527 _kmin = 16 # k = 3, always

528

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

530 '''New L{HeightSmoothBiSpline} interpolator.

531

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

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

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

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

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

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

538

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

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

541

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

543 installed.

544

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

546

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

548 warning as exception.

549 '''

550 spi = self.scipy_interpolate

551

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

553

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

555 try:

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

557 eps=EPS, s=s).ev

558 except Exception as x:

559 raise _SciPyIssue(x)

560

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

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

563

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

565

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

567 warning as exception.

568

569 @see: L{Here<_HeightsBase.__call__>} for further details.

570 '''

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

572

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

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

575

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

577

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

579 warning as exception.

580

581 @see: L{Here<_HeightsBase.height>} for further details.

582 '''

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

584

585

586 class _HeightIDW(_HeightBase):

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

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

589 interpolators.

590

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

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

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

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

595 classes.

596 '''

597 _beta = 0 # fidw inverse power

598 _func = None # formy function

599 _knots = () # knots list or tuple

600 _kwds = {} # func_ options

601

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

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

604

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

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

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

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

609 retrievable with property C{kwds}.

610

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

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

613 '''

614 name, kwds = _name2__(**name__kwds)

615 if name:

616 self.name = name

617

618 n, self._knots = len2(knots)

619 if n < self.kmin:

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

621 self.beta = beta

622 self._kwds = kwds or {}

623

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

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

626

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

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

629 locations (C{bool}).

630

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

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

633

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

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

636 issue.

637 '''

638 def _xy2(wrap=False):

639 _w = _Wrap._latlonop(wrap)

640 try: # like _xyhs above, but degrees

641 for i, ll in enumerate(llis):

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

643 except Exception as x:

644 i = Fmt.INDEX(llis=i)

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

646

647 _as, llis = _as_llis2(llis)

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

649

650 @property_RO

651 def adjust(self):

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

653 '''

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

655

656 @property

657 def beta(self):

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

659 '''

660 return self._beta

661

662 @beta.setter # PYCHOK setter!

663 def beta(self, beta):

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

665

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

667 '''

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

669

670 @property_RO

671 def datum(self):

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

673 '''

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

675

676 def _datum_setter(self, datum):

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

678 '''

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

680 if d and d is not self._datum:

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

682

683 def _distances(self, x, y):

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

685 '''

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

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

688 self._notOverloaded(distance_function=_f)

689 try:

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

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

692 except Exception as e:

693 i = Fmt.INDEX(knots=i)

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

695

696 @property_RO

697 def _fmath(self):

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

699 '''

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

701 return f

702

703 @property_RO

704 def _formy(self):

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

706 '''

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

708 return f

709

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

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

712

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

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

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

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

717

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

719 interpolated heights (C{float}s).

720

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

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

723 '''

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

725

726 @Property_RO

727 def _heights(self):

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

729 '''

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

731

732 def _hIDW(self, x, y):

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

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

735 '''

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

737 try:

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

739 except (TypeError, ValueError) as e:

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

741

742 @property_RO

743 def hypot(self):

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

745 '''

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

747

748 @property_RO

749 def knots(self):

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

751 '''

752 return self._knots

753

754 @property_RO

755 def kwds(self):

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

757 '''

758 return self._kwds

759

760 @property_RO

761 def limit(self):

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

763 '''

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

765

766 @property_RO

767 def radius(self):

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

769 '''

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

771

772 @property_RO

773 def scaled(self):

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

775 '''

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

777

778 @property_RO

779 def wrap(self):

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

781 '''

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

783

784

785 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

789 L{pygeodesy.cosineAndoyerLambert_}.

790 '''

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

792 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

793

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

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

796 function L{pygeodesy.cosineAndoyerLambert}.

797

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

799 '''

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

801 self._func = self._formy.cosineAndoyerLambert

802

803 if _FOR_DOCS:

804 __call__ = _HeightIDW.__call__

805 height = _HeightIDW.height

806

807

808 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

812 L{pygeodesy.cosineForsytheAndoyerLambert_}.

813 '''

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

815 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

816

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

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

819 function L{pygeodesy.cosineForsytheAndoyerLambert}.

820

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

822 '''

823 _HeightIDW.__init__(self, knots, beta=beta, **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__radius_wrap):

839 '''New L{HeightIDWcosineLaw} interpolator.

840

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

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

843 function L{pygeodesy.cosineLaw}.

844

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

846 '''

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

848 self._func = self._formy.cosineLaw

849

850 if _FOR_DOCS:

851 __call__ = _HeightIDW.__call__

852 height = _HeightIDW.height

853

854

855 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

859 conventionally in C{meter}.

860 '''

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

862 '''New L{HeightIDWdistanceTo} interpolator.

863

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

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

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

867

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

869

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

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

872 checked}.

873 '''

874 _HeightIDW.__init__(self, knots, beta=beta, **name__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.equirectangular4}.

902 '''

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

904 '''New L{HeightIDWequirectangular} interpolator.

905

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

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

908 for function L{pygeodesy.equirectangular4}.

909

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

911 '''

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

913

914 def _distances(self, x, y):

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

916 '''

917 _f, kwds = self._formy.equirectangular4, 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__adjust_radius_wrap):

936 '''New L{HeightIDWeuclidean} interpolator.

937

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

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

940 for function function L{pygeodesy.euclidean}.

941

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

943 '''

944 _HeightIDW.__init__(self, knots, beta=beta, **name__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, datum=None, **name__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 name__wrap: Optional C{B{name}=NN} for this height interpolator

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

966 class L{geodesicx.GeodesicExact}.

967

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

969

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

971 '''

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

973 self._datum_setter(datum)

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

975

976 if _FOR_DOCS:

977 __call__ = _HeightIDW.__call__

978 height = _HeightIDW.height

979

980

981 class HeightIDWflatLocal(_HeightIDW):

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

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

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

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

986 '''

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

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

989

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

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

992 keyword arguments for L{pygeodesy.flatLocal}.

993

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

995 '''

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

997 **name__datum_hypot_scaled_wrap)

998 self._func = self._formy.flatLocal

999

1000 if _FOR_DOCS:

1001 __call__ = _HeightIDW.__call__

1002 height = _HeightIDW.height

1003

1004

1005 class HeightIDWflatPolar(_HeightIDW):

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

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

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

1009 L{pygeodesy.flatPolar_}.

1010 '''

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

1012 '''New L{HeightIDWflatPolar} interpolator.

1013

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

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

1016 arguments for function L{pygeodesy.flatPolar}.

1017

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

1019 '''

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

1021 self._func = self._formy.flatPolar

1022

1023 if _FOR_DOCS:

1024 __call__ = _HeightIDW.__call__

1025 height = _HeightIDW.height

1026

1027

1028 class HeightIDWhaversine(_HeightIDW):

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

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

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

1032

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

1034 '''

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

1036 '''New L{HeightIDWhaversine} interpolator.

1037

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

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

1040 arguments for function L{pygeodesy.haversine}.

1041

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

1043 '''

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

1045 self._func = self._formy.haversine

1046

1047 if _FOR_DOCS:

1048 __call__ = _HeightIDW.__call__

1049 height = _HeightIDW.height

1050

1051

1052 class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny

1053 if _FOR_DOCS:

1054 __doc__ = HeightIDWflatLocal.__doc__

1055 __init__ = HeightIDWflatLocal.__init__

1056 __call__ = HeightIDWflatLocal.__call__

1057 height = HeightIDWflatLocal.height

1058

1059

1060 class HeightIDWkarney(_HeightIDW):

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

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

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

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

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

1066 geographiclib.geodesic.Geodesic.Inverse>}.

1067 '''

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

1069 '''New L{HeightIDWkarney} interpolator.

1070

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

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

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

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

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

1076 class L{geodesicw.Geodesic}.

1077

1078 @raise ImportError: Package U{geographiclib

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

1080

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

1082

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

1084 '''

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

1086 self._datum_setter(datum)

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

1088

1089 if _FOR_DOCS:

1090 __call__ = _HeightIDW.__call__

1091 height = _HeightIDW.height

1092

1093

1094 class HeightIDWthomas(_HeightIDW):

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

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

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

1098 '''

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

1100 '''New L{HeightIDWthomas} interpolator.

1101

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

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

1104 arguments for function L{pygeodesy.thomas}.

1105

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

1107 '''

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

1109 self._func = self._formy.thomas

1110

1111 if _FOR_DOCS:

1112 __call__ = _HeightIDW.__call__

1113 height = _HeightIDW.height

1114

1115

1116 class HeightIDWvincentys(_HeightIDW):

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

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

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

1120

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

1122 '''

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

1124 '''New L{HeightIDWvincentys} interpolator.

1125

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

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

1128 arguments for function L{pygeodesy.vincentys}.

1129

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

1131 '''

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

1133 self._func = self._formy.vincentys

1134

1135 if _FOR_DOCS:

1136 __call__ = _HeightIDW.__call__

1137 height = _HeightIDW.height

1138

1139

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

1141

1142# **) MIT License

1143#

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

1145#

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

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

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

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

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

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

1152#

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

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

1155#

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

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

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

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

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

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

1162# OTHER DEALINGS IN THE SOFTWARE.