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.into

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

95

96_error_ = 'error'

97_formy = _MODS.into(formy=__name__)

98_llis_ = 'llis'

99_smoothing_ = 'smoothing'

100

101

102 class HeightError(PointsError):

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

104 '''

105 pass

106

107

108 def _alist(ais):

109 # return list of floats, not numpy.float64s

110 return list(map(float, ais))

111

112

113 def _ascalar(ais): # in .geoids

114 # return single float, not numpy.float64

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

116 if len(ais) != 1:

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

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

119 raise _AssertionError(n, txt=t)

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

121

122

123 def _atuple(ais):

124 # return tuple of floats, not numpy.float64s

125 return tuple(map(float, ais))

126

127

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

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

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

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

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

133

134 n = len(llis)

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

136 llis = llis[0]

137 try:

138 n, llis = len2(llis)

139 _as = _alist # return list of interpolated heights

140 except TypeError: # single lli

141 n, llis = 1, [llis]

142 _as = _ascalar # return single interpolated heights

143 else: # of 0, 2 or more llis

144 _as = _atuple # return tuple of interpolated heights

145

146 if n < m:

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

148 return _as, llis

149

150

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

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

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

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

155 else:

156 n, lats = len2(lats)

157 m, lons = len2(lons)

158 if n != m:

159 # format a LenError, but raise an Error

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

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

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

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

164

165

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

167 # create an insufficient Error instance

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

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

170

171

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

173 # clip, order and remove duplicates

174 # p, ks = 0, []

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

176 # if k > p:

177 # ks.append(k)

178 # p = k

179 # return ks

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

181

182

183 def _xyhs(wrap=False, **name_lls):

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

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

186 name, lls = _xkwds_item2(name_lls)

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

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

189 try:

190 for i, ll in enumerate(lls):

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

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

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

194 except Exception as e:

195 i = Fmt.INDEX(name, i)

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

197

198

199 class _HeightBase(_Named): # in .geoids

200 '''(INTERNAL) Interpolator base class.

201 '''

202 _datum = _WGS84 # default

203 _kmin = 2 # min number of knots

204 _LLis = LatLon_ # ._height class

205 _np_sp = None # (numpy, scipy)

206 _wrap = None # wrap knots and llis

207

208 @property_RO

209 def datum(self):

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

211 '''

212 return self._datum

213

214 @property_RO

215 def kmin(self):

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

217 '''

218 return self._kmin

219

220 @property_RO

221 def wrap(self):

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

223 '''

224 return self._wrap

225

226

227 class _HeightsBase(_HeightBase): # in .geoids

228 '''(INTERNAL) Interpolator base class.

229 '''

230 _np_sp = None # (numpy, scipy)

231

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

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

234

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

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

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

238 setting.

239

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

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

242

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

244 an invalid B{C{lli}}.

245

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

247

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

249 '''

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

251

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

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

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

255 atype = self.numpy.array

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

257 _as, llis = _as_llis2(llis)

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

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

260

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

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

263 self._notOverloaded(*args)

264

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

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

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

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

269 except Exception as x:

270 raise _SciPyIssue(x)

271

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

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

274

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

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

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

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

279 the B{C{knots}} setting.

280

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

282 interpolated heights (C{float}s).

283

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

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

286

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

288

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

290 '''

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

292

293 def _np_sp2(self, throwarnings=False):

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

295 '''

296 t = _HeightsBase._np_sp

297 if not t:

298 # raise SciPyWarnings, but not if

299 # scipy has already been imported

300 if throwarnings: # PYCHOK no cover

301 import sys

302 if _scipy_ not in sys.modules:

303 import warnings

304 warnings.filterwarnings(_error_)

305

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

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

308

309 _HeightsBase._np_sp = t = np, sp

310 return t

311

312 @Property_RO

313 def numpy(self):

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

315 '''

316 np, _ = self._np_sp2()

317 return np

318

319 @Property_RO

320 def scipy(self):

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

322 '''

323 _, sp = self._np_sp2()

324 return sp

325

326 @Property_RO

327 def scipy_interpolate(self):

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

329 '''

330 _ = self.scipy

331 import scipy.interpolate as spi # scipy 1.2.2

332 return spi

333

334 def _scipy_version(self, **n):

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

336 '''

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

338

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

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

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

342 n = len(hs)

343 if n < self.kmin:

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

345 if name:

346 self.name = name

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

348

349

350 class HeightCubic(_HeightsBase):

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

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

353 C{kind='cubic'}.

354 '''

355 _interp2d = None

356 _kind = _cubic_

357 _kmin = 16

358

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

360 '''New L{HeightCubic} interpolator.

361

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

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

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

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

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

367

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

369 invalid B{C{knot}}.

370

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

372 or not installed.

373

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

375

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

377 as exception.

378 '''

379 spi = self.scipy_interpolate

380

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

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

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

384 except Exception as x:

385 raise _SciPyIssue(x)

386

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

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

389

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

391

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

393 as exception.

394

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

396 '''

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

398

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

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

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

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

403

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

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

406

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

408

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

410 as exception.

411

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

413 '''

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

415

416

417 class HeightLinear(HeightCubic):

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

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

420 C{kind='linear'}.

421 '''

422 _kind = _linear_

423 _kmin = 2

424

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

426 '''New L{HeightLinear} interpolator.

427

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

429 '''

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

431

432 if _FOR_DOCS:

433 __call__ = HeightCubic.__call__

434 height = HeightCubic.height

435

436

437 class HeightLSQBiSpline(_HeightsBase):

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

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

440 interpolate.LSQSphereBivariateSpline.html>}.

441 '''

442 _kmin = 16 # k = 3, always

443

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

445 '''New L{HeightLSQBiSpline} interpolator.

446

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

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

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

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

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

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

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

454

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

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

457

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

459

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

461 not installed.

462

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

464

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

466 warning as exception.

467 '''

468 np = self.numpy

469 spi = self.scipy_interpolate

470

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

472 n = len(hs)

473

474 w = weight

475 if isscalar(w):

476 w = float(w)

477 if w <= 0:

478 raise HeightError(weight=w)

479 w = [w] * n

480 elif w is not None:

481 m, w = len2(w)

482 if m != n:

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

484 w = map2(float, w)

485 m = min(w)

486 if m <= 0: # PYCHOK no cover

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

488 raise HeightError(w, m)

489 try:

490 T = 1.0e-4 # like SciPy example

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

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

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

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

495 except Exception as x:

496 raise _SciPyIssue(x)

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 further 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 further 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_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_wrap: Optional C{B{name}=NN} for this height interpolator

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

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

538 C{True} (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 or not

544 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, **name_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 def __call__(self, *llis, **wrap):

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

564

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

566

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

568 warning as exception.

569

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

571 '''

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

573

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

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

576

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

578

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

580 warning as exception.

581

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

583 '''

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

585

586

587 class _HeightIDW(_HeightBase):

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

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

590 interpolators.

591

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

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

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

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

596 classes.

597 '''

598 _beta = 0 # fidw inverse power

599 _func = None # formy function

600 _knots = () # knots list or tuple

601 _kwds = {} # func_ options

602

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

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

605

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

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

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

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

610 retrievable with property C{kwds}.

611

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

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

614 '''

615 name, kwds = _name2__(**name__kwds)

616 if name:

617 self.name = name

618

619 n, self._knots = len2(knots)

620 if n < self.kmin:

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

622 self.beta = beta

623 self._kwds = kwds or {}

624

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

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

627

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

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

630 locations (C{bool}).

631

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

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

634

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

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

637 issue.

638 '''

639 def _xy2(wrap=False):

640 _w = _Wrap._latlonop(wrap)

641 try: # like _xyhs above, but degrees

642 for i, ll in enumerate(llis):

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

644 except Exception as x:

645 i = Fmt.INDEX(llis=i)

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

647

648 _as, llis = _as_llis2(llis)

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

650

651 @property_RO

652 def adjust(self):

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

654 '''

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

656

657 @property

658 def beta(self):

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

660 '''

661 return self._beta

662

663 @beta.setter # PYCHOK setter!

664 def beta(self, beta):

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

666

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

668 '''

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

670

671 @property_RO

672 def datum(self):

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

674 '''

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

676

677 def _datum_setter(self, datum):

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

679 '''

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

681 if d and d is not self._datum:

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

683

684 def _distances(self, x, y):

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

686 '''

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

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

689 self._notOverloaded(distance_function=_f)

690 try:

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

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

693 except Exception as e:

694 i = Fmt.INDEX(knots=i)

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

696

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

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

699

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

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

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

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

704

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

706 interpolated heights (C{float}s).

707

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

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

710 '''

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

712

713 @Property_RO

714 def _heights(self):

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

716 '''

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

718

719 def _hIDW(self, x, y):

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

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

722 '''

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

724 try:

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

726 except (TypeError, ValueError) as e:

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

728

729 @property_RO

730 def hypot(self):

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

732 '''

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

734

735 @property_RO

736 def knots(self):

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

738 '''

739 return self._knots

740

741 @property_RO

742 def kwds(self):

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

744 '''

745 return self._kwds

746

747 @property_RO

748 def limit(self):

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

750 '''

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

752

753 @property_RO

754 def radius(self):

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

756 '''

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

758

759 @property_RO

760 def scaled(self):

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

762 '''

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

764

765 @property_RO

766 def wrap(self):

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

768 '''

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

770

771

772 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

776 L{pygeodesy.cosineAndoyerLambert_}.

777 '''

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

779 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

780

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

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

783 function L{pygeodesy.cosineAndoyerLambert}.

784

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

786 '''

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

788 self._func = _formy.cosineAndoyerLambert

789

790 if _FOR_DOCS:

791 __call__ = _HeightIDW.__call__

792 height = _HeightIDW.height

793

794

795 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

799 L{pygeodesy.cosineForsytheAndoyerLambert_}.

800 '''

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

802 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

803

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

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

806 function L{pygeodesy.cosineForsytheAndoyerLambert}.

807

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

809 '''

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

811 self._func = _formy.cosineForsytheAndoyerLambert

812

813 if _FOR_DOCS:

814 __call__ = _HeightIDW.__call__

815 height = _HeightIDW.height

816

817

818 class HeightIDWcosineLaw(_HeightIDW):

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

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

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

822

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

824 '''

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

826 '''New L{HeightIDWcosineLaw} interpolator.

827

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

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

830 function L{pygeodesy.cosineLaw}.

831

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

833 '''

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

835 self._func = _formy.cosineLaw

836

837 if _FOR_DOCS:

838 __call__ = _HeightIDW.__call__

839 height = _HeightIDW.height

840

841

842 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

846 conventionally in C{meter}.

847 '''

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

849 '''New L{HeightIDWdistanceTo} interpolator.

850

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

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

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

854

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

856

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

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

859 checked}.

860 '''

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

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

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

864 # in _height_called instead of class LatLon_ and datum None

865 self._datum = ks0.datum

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

867

868 def _distances(self, x, y):

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

870 '''

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

872 try:

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

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

875 except Exception as e:

876 i = Fmt.INDEX(knots=i)

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

878

879 if _FOR_DOCS:

880 __call__ = _HeightIDW.__call__

881 height = _HeightIDW.height

882

883

884 class HeightIDWequirectangular(_HeightIDW):

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

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

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

888 L{pygeodesy.equirectangular4}.

889 '''

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

891 '''New L{HeightIDWequirectangular} interpolator.

892

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

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

895 for function L{pygeodesy.equirectangular4}.

896

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

898 '''

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

900

901 def _distances(self, x, y):

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

903 '''

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

905 try:

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

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

908 except Exception as e:

909 i = Fmt.INDEX(knots=i)

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

911

912 if _FOR_DOCS:

913 __call__ = _HeightIDW.__call__

914 height = _HeightIDW.height

915

916

917 class HeightIDWeuclidean(_HeightIDW):

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

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

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

921 '''

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

923 '''New L{HeightIDWeuclidean} interpolator.

924

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

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

927 for function function L{pygeodesy.euclidean}.

928

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

930 '''

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

932 self._func = _formy.euclidean

933

934 if _FOR_DOCS:

935 __call__ = _HeightIDW.__call__

936 height = _HeightIDW.height

937

938

939 class HeightIDWexact(_HeightIDW):

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

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

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

943 L{GeodesicExact.Inverse}.

944 '''

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

946 '''New L{HeightIDWexact} interpolator.

947

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

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

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

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

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

953 class L{geodesicx.GeodesicExact}.

954

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

956

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

958 '''

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

960 self._datum_setter(datum)

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

962

963 if _FOR_DOCS:

964 __call__ = _HeightIDW.__call__

965 height = _HeightIDW.height

966

967

968 class HeightIDWflatLocal(_HeightIDW):

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

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

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

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

973 '''

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

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

976

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

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

979 keyword arguments for L{pygeodesy.flatLocal}.

980

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

982 '''

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

984 **name__datum_hypot_scaled_wrap)

985 self._func = _formy.flatLocal

986

987 if _FOR_DOCS:

988 __call__ = _HeightIDW.__call__

989 height = _HeightIDW.height

990

991

992 class HeightIDWflatPolar(_HeightIDW):

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

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

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

996 L{pygeodesy.flatPolar_}.

997 '''

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

999 '''New L{HeightIDWflatPolar} interpolator.

1000

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

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

1003 arguments for function L{pygeodesy.flatPolar}.

1004

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

1006 '''

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

1008 self._func = _formy.flatPolar

1009

1010 if _FOR_DOCS:

1011 __call__ = _HeightIDW.__call__

1012 height = _HeightIDW.height

1013

1014

1015 class HeightIDWhaversine(_HeightIDW):

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

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

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

1019

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

1021 '''

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

1023 '''New L{HeightIDWhaversine} interpolator.

1024

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

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

1027 arguments for function L{pygeodesy.haversine}.

1028

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

1030 '''

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

1032 self._func = _formy.haversine

1033

1034 if _FOR_DOCS:

1035 __call__ = _HeightIDW.__call__

1036 height = _HeightIDW.height

1037

1038

1039 class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny

1040 if _FOR_DOCS:

1041 __doc__ = HeightIDWflatLocal.__doc__

1042 __init__ = HeightIDWflatLocal.__init__

1043 __call__ = HeightIDWflatLocal.__call__

1044 height = HeightIDWflatLocal.height

1045

1046

1047 class HeightIDWkarney(_HeightIDW):

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

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

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

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

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

1053 geographiclib.geodesic.Geodesic.Inverse>}.

1054 '''

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

1056 '''New L{HeightIDWkarney} interpolator.

1057

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

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

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

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

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

1063 class L{geodesicw.Geodesic}.

1064

1065 @raise ImportError: Package U{geographiclib

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

1067

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

1069

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

1071 '''

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

1073 self._datum_setter(datum)

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

1075

1076 if _FOR_DOCS:

1077 __call__ = _HeightIDW.__call__

1078 height = _HeightIDW.height

1079

1080

1081 class HeightIDWthomas(_HeightIDW):

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

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

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

1085 '''

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

1087 '''New L{HeightIDWthomas} interpolator.

1088

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

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

1091 arguments for function L{pygeodesy.thomas}.

1092

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

1094 '''

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

1096 self._func = _formy.thomas

1097

1098 if _FOR_DOCS:

1099 __call__ = _HeightIDW.__call__

1100 height = _HeightIDW.height

1101

1102

1103 class HeightIDWvincentys(_HeightIDW):

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

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

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

1107

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

1109 '''

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

1111 '''New L{HeightIDWvincentys} interpolator.

1112

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

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

1115 arguments for function L{pygeodesy.vincentys}.

1116

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

1118 '''

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

1120 self._func = _formy.vincentys

1121

1122 if _FOR_DOCS:

1123 __call__ = _HeightIDW.__call__

1124 height = _HeightIDW.height

1125

1126

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

1128

1129# **) MIT License

1130#

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

1132#

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

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

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

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

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

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

1139#

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

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

1142#

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

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

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

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

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

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

1149# OTHER DEALINGS IN THE SOFTWARE.