1

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

3

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

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

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

20 particular way.

21

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

23

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

25

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

27 arguments to select different interpolation options

28

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

30

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

32

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

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

35

36 or

37

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

39

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

41

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

43

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

45

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

47

48 or as 2 lists, 2 tuples, etc.

49

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

51

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

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

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

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

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

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

58

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

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

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

62

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

64 Interpolation.

65'''

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

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

68

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

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

71 from pygeodesy.datums import _ellipsoidal_datum, _WGS84

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

73 _SciPyIssue, _xattr, _xkwds, _xkwds_get

74 from pygeodesy.fmath import fidw, Float_, Int_

75 from pygeodesy.formy import cosineAndoyerLambert, cosineForsytheAndoyerLambert, \

76 cosineLaw, equirectangular_, euclidean, flatLocal, \

77 flatPolar, haversine, thomas, vincentys, radians

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

79 _linear_, _NOTEQUAL_, _PLUS_, _scipy_, _SPACE_, _STAR_

80 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS

81 from pygeodesy.named import _Named, notOverloaded

82 from pygeodesy.points import _distanceTo, LatLon_, Fmt, _Wrap

83 from pygeodesy.props import Property_RO, property_RO

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

85# from pygeodesy.units import Float_, Int_ # from .fmath

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

87

88# from math import radians # from .formy

89

90__all__ = _ALL_LAZY.heights

91__version__ = '23.08.06'

92

93_error_ = 'error'

94_llis_ = 'llis'

95_smoothing_ = 'smoothing'

96

97

98 class HeightError(PointsError):

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

100 '''

101 pass

102

103

104 def _alist(ais):

105 # return list of floats, not numpy.float64s

106 return list(map(float, ais))

107

108

109 def _ascalar(ais): # in .geoids

110 # return single float, not numpy.float64

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

112 if len(ais) != 1:

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

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

115 raise _AssertionError(n, txt=t)

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

117

118

119 def _atuple(ais):

120 # return tuple of floats, not numpy.float64s

121 return tuple(map(float, ais))

122

123

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

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

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

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

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

129

130 n = len(llis)

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

132 llis = llis[0]

133 try:

134 n, llis = len2(llis)

135 _as = _alist # return list of interpolated heights

136 except TypeError: # single lli

137 n, llis = 1, [llis]

138 _as = _ascalar # return single interpolated heights

139 else: # of 0, 2 or more llis

140 _as = _atuple # return tuple of interpolated heights

141

142 if n < m:

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

144 return _as, llis

145

146

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

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

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

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

151 else:

152 n, lats = len2(lats)

153 m, lons = len2(lons)

154 if n != m:

155 # format a LenError, but raise an Error

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

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

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

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

160

161

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

163 # create an insufficient Error instance

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

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

166

167

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

169 # clip, order and remove duplicates

170 # p, ks = 0, []

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

172 # if k > p:

173 # ks.append(k)

174 # p = k

175 # return ks

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

177

178

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

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

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

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

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

184 try:

185 for i, ll in enumerate(lls):

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

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

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

189 except AttributeError as x:

190 i = Fmt.SQUARE(name, i)

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

192

193

194 class _HeightBase(_Named): # in .geoids

195 '''(INTERNAL) Interpolator base class.

196 '''

197 _datum = _WGS84 # default

198 _kmin = 2 # min number of knots

199 _LLis = LatLon_ # ._height class

200 _np_sp = None # (numpy, scipy)

201 _wrap = None # wrap knots and llis

202

203 @property_RO

204 def datum(self):

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

206 '''

207 return self._datum

208

209 @property_RO

210 def kmin(self):

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

212 '''

213 return self._kmin

214

215 @property_RO

216 def wrap(self):

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

218 '''

219 return self._wrap

220

221

222 class _HeightsBase(_HeightBase): # in .geoids

223 '''(INTERNAL) Interpolator base class.

224 '''

225 _np_sp = None # (numpy, scipy)

226

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

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

229

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

231 C{LatLon}s).

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 notOverloaded(self, 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, **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}, see function C{notOverloaded}.

259 '''

260 notOverloaded(self, *args)

261

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

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

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

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

266 except Exception as x:

267 raise _SciPyIssue(x)

268

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

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

271

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

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

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

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

276 the B{C{knots}} setting.

277

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

279 interpolated heights (C{float}s).

280

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

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

283

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

285

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

287 '''

288 notOverloaded(self, lats, lons, **wrap)

289

290 def _np_sp2(self, throwarnings=False):

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

292 '''

293 t = _HeightsBase._np_sp

294 if not t:

295 # raise SciPyWarnings, but not if

296 # scipy has been imported already

297 if throwarnings: # PYCHOK no cover

298 import sys

299 if _scipy_ not in sys.modules:

300 import warnings

301 warnings.filterwarnings(_error_)

302

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

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

305

306 _HeightsBase._np_sp = t = np, sp

307 return t

308

309 @Property_RO

310 def numpy(self):

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

312 '''

313 np, _ = self._np_sp2()

314 return np

315

316 @Property_RO

317 def scipy(self):

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

319 '''

320 _, sp = self._np_sp2()

321 return sp

322

323 @Property_RO

324 def scipy_interpolate(self):

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

326 '''

327 _ = self.scipy

328 import scipy.interpolate as spi

329 return spi

330

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

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

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

334 n = len(hs)

335 if n < self.kmin:

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

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

338

339

340 class HeightCubic(_HeightsBase):

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

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

343 C{kind='cubic'}.

344 '''

345 _interp2d = None

346 _kind = _cubic_

347 _kmin = 16

348

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

350 '''New L{HeightCubic} interpolator.

351

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

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

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

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

356

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

358 invalid B{C{knot}}.

359

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

361 or not installed.

362

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

364

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

366 as exception.

367 '''

368 spi = self.scipy_interpolate

369

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

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

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

373 except Exception as x:

374 raise _SciPyIssue(x)

375

376 if name:

377 self.name = name

378

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

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

381

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

383

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

385 as exception.

386

387 @see: L{_HeightsBase.__call__} for details about B{C{llis}},

388 B{C{wrap}}, return values and other exceptions.

389 '''

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

391

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

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

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

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

396

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

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

399

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

401

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

403 as exception.

404

405 @see: L{_HeightsBase.height} for details about B{C{lats}},

406 B{C{lons}}, B{C{wrap}}, return values and other exceptions.

407 '''

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

409

410

411 class HeightLinear(HeightCubic):

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

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

414 C{kind='linear'}.

415 '''

416 _kind = _linear_

417 _kmin = 2

418

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

420 '''New L{HeightLinear} interpolator.

421

422 @see: L{HeightCubic.__init__} for details about B{C{knots}},

423 B{C{name}}, B{C{wrap}} and other exceptions.

424 '''

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

426

427 if _FOR_DOCS:

428 __call__ = HeightCubic.__call__

429 height = HeightCubic.height

430

431

432 class HeightLSQBiSpline(_HeightsBase):

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

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

435 interpolate.LSQSphereBivariateSpline.html>}.

436 '''

437 _kmin = 16 # k = 3, always

438

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

440 '''New L{HeightLSQBiSpline} interpolator.

441

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

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

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

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

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

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

448

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

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

451

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

453

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

455 or not installed.

456

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

458

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

460 warning as exception.

461 '''

462 np = self.numpy

463 spi = self.scipy_interpolate

464

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

466 n = len(hs)

467

468 w = weight

469 if isscalar(w):

470 w = float(w)

471 if w <= 0:

472 raise HeightError(weight=w)

473 w = [w] * n

474 elif w is not None:

475 m, w = len2(w)

476 if m != n:

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

478 w = map2(float, w)

479 m = min(w)

480 if m <= 0: # PYCHOK no cover

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

482 raise HeightError(w, m)

483 try:

484 T = 1.0e-4 # like SciPy example

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

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

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

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

489 except Exception as x:

490 raise _SciPyIssue(x)

491

492 if name:

493 self.name = name

494

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

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

497

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

499

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

501 warning as exception.

502

503 @see: L{_HeightsBase.__call__} for details about B{C{llis}},

504 B{C{wrap}}, return values and other exceptions.

505 '''

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

507

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

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

510

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

512

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

514 warning as exception.

515

516 @see: L{_HeightsBase.height} for details about B{C{lats}},

517 B{C{lons}}, B{C{wrap}}, return values and other exceptions.

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=NN, **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: Optional name for this height interpolator (C{str}).

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

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

537 and B{C{lons}} locations (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

543 or not 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, **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 if name:

562 self.name = name

563

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

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

566

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

568

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

570 warning as exception.

571

572 @see: L{_HeightsBase.__call__} for details about B{C{llis}},

573 B{C{wrap}}, return values and other exceptions.

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{_HeightsBase.height} for details about B{C{lats}},

586 B{C{lons}}, B{C{wrap}}, return values and other exceptions.

587 '''

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

589

590

591 class _HeightIDW(_HeightBase):

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

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

594 interpolators.

595

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

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

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

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

600 classes.

601 '''

602 _beta = 0 # fidw inverse power

603 _func = None # formy function

604 _ks = () # knots list or tuple

605 _kwds = {} # func_ options

606

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

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

609

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

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

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

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

614 retrievable with property C{kwds}.

615

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

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

618 '''

619 if name:

620 self.name = name

621 n, self._ks = len2(knots)

622 if n < self.kmin:

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

624 self.beta = beta

625 self._kwds = kwds or {}

626

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

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

629

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

631 C{LatLon}s).

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

633 locations (C{bool}).

634

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

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

637

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

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

640 issue.

641 '''

642 def _xy2(lls, wrap=False):

643 _w = _Wrap._latlonop(wrap)

644 try: # like _xyhs above, but degrees

645 for i, ll in enumerate(lls):

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

647 except AttributeError as x:

648 i = Fmt.SQUARE(llis=i)

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

650

651 _as, llis = _as_llis2(llis)

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

653

654 @property_RO

655 def adjust(self):

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

657 '''

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

659

660 @property

661 def beta(self):

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

663 '''

664 return self._beta

665

666 @beta.setter # PYCHOK setter!

667 def beta(self, beta):

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

669

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

671 '''

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

673

674 @property_RO

675 def datum(self):

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

677 '''

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

679

680 def _datum_setter(self, datum):

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

682 '''

683 d = datum or _xattr(self._ks[0], datum=None)

684 if d and d is not self._datum:

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

686

687 def _distances(self, x, y):

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

689 '''

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

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

692 notOverloaded(self, distance_function=_f)

693 for _, k in enumerate(self._ks):

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

695

696 def _hIDW(self, x, y):

697 '''(INTERNAL) Return the IDW-interpolate height at

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

699 '''

700 try:

701 ds = self._distances(x, y)

702 hs = (k.height for k in self._ks)

703 return fidw(hs, ds, beta=self._beta)

704 except (TypeError, ValueError) as x:

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

706

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

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

709

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

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

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

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

714

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

716 interpolated heights (C{float}s).

717

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

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

720 issue.

721 '''

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

723

724 @property_RO

725 def hypot(self):

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

727 '''

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

729

730 @property_RO

731 def knots(self):

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

733 '''

734 return self._ks

735

736 @property_RO

737 def kwds(self):

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

739 '''

740 return self._kwds

741

742 @property_RO

743 def limit(self):

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

745 '''

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

747

748 @property_RO

749 def radius(self):

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

751 '''

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

753

754 @property_RO

755 def scaled(self):

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

757 '''

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

759

760 @property_RO

761 def wrap(self):

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

763 '''

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

765

766

767 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

771 L{pygeodesy.cosineAndoyerLambert_}.

772 '''

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

774 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

775

776 @kwarg datum_wrap: Optional keyword arguments for function

777 L{pygeodesy.cosineAndoyerLambert}.

778

779 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

780 B{C{beta}}, B{C{name}} and other exceptions.

781 '''

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

783 self._func = cosineAndoyerLambert

784

785 if _FOR_DOCS:

786 __call__ = _HeightIDW.__call__

787 height = _HeightIDW.height

788

789

790 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

794 L{pygeodesy.cosineForsytheAndoyerLambert_}.

795 '''

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

797 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

798

799 @kwarg datum_wrap: Optional keyword arguments for function

800 L{pygeodesy.cosineAndoyerLambert}.

801

802 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

803 B{C{beta}}, B{C{name}} and other exceptions.

804 '''

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

806 self._func = cosineForsytheAndoyerLambert

807

808 if _FOR_DOCS:

809 __call__ = _HeightIDW.__call__

810 height = _HeightIDW.height

811

812

813 class HeightIDWcosineLaw(_HeightIDW):

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

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

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

817

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

819 '''

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

821 '''New L{HeightIDWcosineLaw} interpolator.

822

823 @kwarg radius_wrap: Optional keyword arguments for function

824 L{pygeodesy.cosineLaw}.

825

826 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

827 B{C{beta}}, B{C{name}} and other exceptions.

828 '''

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

830 self._func = cosineLaw

831

832 if _FOR_DOCS:

833 __call__ = _HeightIDW.__call__

834 height = _HeightIDW.height

835

836

837 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

841 conventionally in C{meter}.

842 '''

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

844 '''New L{HeightIDWdistanceTo} interpolator.

845

846 @kwarg distanceTo_kwds: Optional keyword arguments for the

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

848 method.

849

850 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

851 B{C{beta}}, B{C{name}} and other exceptions.

852

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

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

855 checked however}.

856 '''

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

858 **distanceTo_kwds)

859 ks = _distanceTo(HeightError, knots=self._ks)

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

861 # compatible points in _height_called

862 # instead of class LatLon_ and datum None

863 self._datum = ks[0].datum

864 self._LLis = ks[0].classof

865

866 def _distances(self, x, y):

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

868 '''

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

870 for _, k in enumerate(self._ks):

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

872

873 if _FOR_DOCS:

874 __call__ = _HeightIDW.__call__

875 height = _HeightIDW.height

876

877

878 class HeightIDWequirectangular(_HeightIDW):

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

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

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

882 L{pygeodesy.equirectangular_}.

883 '''

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

885 '''New L{HeightIDWequirectangular} interpolator.

886

887 @kwarg adjust_limit_wrap: Optional keyword arguments for

888 function L{pygeodesy.equirectangular_}.

889

890 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

891 B{C{beta}}, B{C{name}} and other exceptions.

892 '''

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

894 **adjust_limit_wrap)

895

896 def _distances(self, x, y):

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

898 '''

899 _f, kwds = equirectangular_, self._kwds

900 for _, k in enumerate(self._ks):

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

902

903 if _FOR_DOCS:

904 __call__ = _HeightIDW.__call__

905 height = _HeightIDW.height

906

907

908 class HeightIDWeuclidean(_HeightIDW):

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

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

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

912 '''

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

914 '''New L{HeightIDWeuclidean} interpolator.

915

916 @kwarg adjust_radius_wrap: Optional keyword arguments for

917 function L{pygeodesy.euclidean_}.

918

919 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

920 B{C{beta}}, B{C{name}} and other exceptions.

921 '''

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

923 **adjust_radius_wrap)

924 self._func = euclidean

925

926 if _FOR_DOCS:

927 __call__ = _HeightIDW.__call__

928 height = _HeightIDW.height

929

930

931 class HeightIDWexact(_HeightIDW):

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

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

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

935 L{GeodesicExact}C{.Inverse}.

936 '''

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

938 '''New L{HeightIDWexact} interpolator.

939

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

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

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

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

944 of class L{geodesicx.GeodesicExact}.

945

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

947

948 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

949 B{C{beta}}, B{C{name}} and other exceptions.

950 '''

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

952 self._datum_setter(datum)

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

954

955 if _FOR_DOCS:

956 __call__ = _HeightIDW.__call__

957 height = _HeightIDW.height

958

959

960 class HeightIDWflatLocal(_HeightIDW):

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

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

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

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

965 '''

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

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

968

969 @kwarg datum_hypot_scaled_wrap: Optional keyword arguments for

970 function L{pygeodesy.flatLocal_}.

971

972 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

973 B{C{beta}}, B{C{name}} and other exceptions.

974 '''

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

976 **datum_hypot_scaled_wrap)

977 self._func = flatLocal

978

979 if _FOR_DOCS:

980 __call__ = _HeightIDW.__call__

981 height = _HeightIDW.height

982

983

984 class HeightIDWflatPolar(_HeightIDW):

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

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

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

988 L{pygeodesy.flatPolar_}.

989 '''

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

991 '''New L{HeightIDWflatPolar} interpolator.

992

993 @kwarg radius_wrap: Optional keyword arguments for function

994 L{pygeodesy.flatPolar}.

995

996 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

997 B{C{beta}}, B{C{name}} and other exceptions.

998 '''

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

1000 self._func = flatPolar

1001

1002 if _FOR_DOCS:

1003 __call__ = _HeightIDW.__call__

1004 height = _HeightIDW.height

1005

1006

1007 class HeightIDWhaversine(_HeightIDW):

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

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

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

1011

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

1013 '''

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

1015 '''New L{HeightIDWhaversine} interpolator.

1016

1017 @kwarg radius_wrap: Optional keyword arguments for function

1018 L{pygeodesy.haversine}.

1019

1020 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

1021 B{C{beta}}, B{C{name}} and other exceptions.

1022 '''

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

1024 self._func = haversine

1025

1026 if _FOR_DOCS:

1027 __call__ = _HeightIDW.__call__

1028 height = _HeightIDW.height

1029

1030

1031 class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny

1032 if _FOR_DOCS:

1033 __doc__ = HeightIDWflatLocal.__doc__

1034 __init__ = HeightIDWflatLocal.__init__

1035 __call__ = HeightIDWflatLocal.__call__

1036 height = HeightIDWflatLocal.height

1037

1038

1039 class HeightIDWkarney(_HeightIDW):

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

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

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

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

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

1045 Inverse method.

1046 '''

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

1048 '''New L{HeightIDWkarney} interpolator.

1049

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

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

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

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

1054 of class L{geodesicw.Geodesic}.

1055

1056 @raise ImportError: Package U{geographiclib

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

1058

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

1060

1061 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

1062 B{C{beta}}, B{C{name}} and other exceptions.

1063 '''

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

1065 self._datum_setter(datum)

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

1067

1068 if _FOR_DOCS:

1069 __call__ = _HeightIDW.__call__

1070 height = _HeightIDW.height

1071

1072

1073 class HeightIDWthomas(_HeightIDW):

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

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

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

1077 '''

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

1079 '''New L{HeightIDWthomas} interpolator.

1080

1081 @kwarg datum_wrap: Optional keyword argument for function

1082 L{pygeodesy.thomas}.

1083

1084 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

1085 B{C{beta}}, B{C{name}} and other exceptions.

1086 '''

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

1088 self._func = thomas

1089

1090 if _FOR_DOCS:

1091 __call__ = _HeightIDW.__call__

1092 height = _HeightIDW.height

1093

1094

1095 class HeightIDWvincentys(_HeightIDW):

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

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

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

1099

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

1101 '''

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

1103 '''New L{HeightIDWvincentys} interpolator.

1104

1105 @kwarg radius_wrap: Optional keyword arguments for function

1106 L{pygeodesy.vincentys}.

1107

1108 @see: L{_HeightIDW.__init__} for details about B{C{knots}},

1109 B{C{beta}}, B{C{name}} and other exceptions.

1110 '''

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

1112 self._func = vincentys

1113

1114 if _FOR_DOCS:

1115 __call__ = _HeightIDW.__call__

1116 height = _HeightIDW.height

1117

1118

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

1120

1121# **) MIT License

1122#

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

1124#

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

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

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

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

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

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

1131#

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

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

1134#

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

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

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

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

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

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

1141# OTHER DEALINGS IN THE SOFTWARE.