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

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: One or more locations (C{LatLon}s), all positional.

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

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

233 setting.

234

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

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

237

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

239 an invalid B{C{lli}}.

240

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

242

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

244 '''

245 notOverloaded(self, callername='__call__', *llis, **wrap)

246

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

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

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

250 atype = self.numpy.array

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

252 _as, llis = _as_llis2(llis)

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

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

255

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

257 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.

258 '''

259 notOverloaded(self, *args)

260

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

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

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

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

265 except Exception as x:

266 raise _SciPyIssue(x)

267

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

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

270

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

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

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

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

275 the B{C{knots}} setting.

276

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

278 interpolated heights (C{float}s).

279

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

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

282

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

284

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

286 '''

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

288

289 def _np_sp2(self, throwarnings=False):

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

291 '''

292 t = _HeightsBase._np_sp

293 if not t:

294 # raise SciPyWarnings, but not if

295 # scipy has been imported already

296 if throwarnings: # PYCHOK no cover

297 import sys

298 if _scipy_ not in sys.modules:

299 import warnings

300 warnings.filterwarnings(_error_)

301

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

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

304

305 _HeightsBase._np_sp = t = np, sp

306 return t

307

308 @Property_RO

309 def numpy(self):

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

311 '''

312 np, _ = self._np_sp2()

313 return np

314

315 @Property_RO

316 def scipy(self):

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

318 '''

319 _, sp = self._np_sp2()

320 return sp

321

322 @Property_RO

323 def scipy_interpolate(self):

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

325 '''

326 _ = self.scipy

327 import scipy.interpolate as spi

328 return spi

329

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

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

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

333 n = len(hs)

334 if n < self.kmin:

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

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

337

338

339 class HeightCubic(_HeightsBase):

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

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

342 C{kind='cubic'}.

343 '''

344 _interp2d = None

345 _kind = _cubic_

346 _kmin = 16

347

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

349 '''New L{HeightCubic} interpolator.

350

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

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

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

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

355

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

357 invalid B{C{knot}}.

358

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

360 or not installed.

361

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

363

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

365 as exception.

366 '''

367 spi = self.scipy_interpolate

368

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

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

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

372 except Exception as x:

373 raise _SciPyIssue(x)

374

375 if name:

376 self.name = name

377

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

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

380

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

382

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

384 as exception.

385

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

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

388 '''

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

390

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

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

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

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

395

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

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

398

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

400

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

402 as exception.

403

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

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

406 '''

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

408

409

410 class HeightLinear(HeightCubic):

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

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

413 C{kind='linear'}.

414 '''

415 _kind = _linear_

416 _kmin = 2

417

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

419 '''New L{HeightLinear} interpolator.

420

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

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

423 '''

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

425

426 if _FOR_DOCS:

427 __call__ = HeightCubic.__call__

428 height = HeightCubic.height

429

430

431 class HeightLSQBiSpline(_HeightsBase):

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

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

434 interpolate.LSQSphereBivariateSpline.html>}.

435 '''

436 _kmin = 16 # k = 3, always

437

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

439 '''New L{HeightLSQBiSpline} interpolator.

440

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

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

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

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

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

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

447

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

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

450

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

452

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

454 or not installed.

455

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

457

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

459 warning as exception.

460 '''

461 np = self.numpy

462 spi = self.scipy_interpolate

463

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

465 n = len(hs)

466

467 w = weight

468 if isscalar(w):

469 w = float(w)

470 if w <= 0:

471 raise HeightError(weight=w)

472 w = [w] * n

473 elif w is not None:

474 m, w = len2(w)

475 if m != n:

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

477 w = map2(float, w)

478 m = min(w)

479 if m <= 0: # PYCHOK no cover

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

481 raise HeightError(w, m)

482 try:

483 T = 1.0e-4 # like SciPy example

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

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

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

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

488 except Exception as x:

489 raise _SciPyIssue(x)

490

491 if name:

492 self.name = name

493

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

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

496

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

498

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

500 warning as exception.

501

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

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

504 '''

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

506

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

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

509

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

511

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

513 warning as exception.

514

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

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

517 '''

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

519

520

521 class HeightSmoothBiSpline(_HeightsBase):

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

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

524 interpolate.SmoothSphereBivariateSpline.html>}.

525 '''

526 _kmin = 16 # k = 3, always

527

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

529 '''New L{HeightSmoothBiSpline} interpolator.

530

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

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

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

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

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

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

537

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

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

540

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

542 or not installed.

543

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

545

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

547 warning as exception.

548 '''

549 spi = self.scipy_interpolate

550

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

552

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

554 try:

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

556 eps=EPS, s=s).ev

557 except Exception as x:

558 raise _SciPyIssue(x)

559

560 if name:

561 self.name = name

562

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

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

565

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

567

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

569 warning as exception.

570

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

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

573 '''

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

575

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

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

578

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

580

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

582 warning as exception.

583

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

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

586 '''

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

588

589

590 class _HeightIDW(_HeightBase):

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

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

593 interpolators.

594

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

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

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

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

599 classes.

600 '''

601 _beta = 0 # fidw inverse power

602 _func = None # formy function

603 _ks = () # knots list or tuple

604 _kwds = {} # func_ options

605

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

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

608

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

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

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

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

613 retrievable with property C{kwds}.

614

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

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

617 '''

618 if name:

619 self.name = name

620 n, self._ks = len2(knots)

621 if n < self.kmin:

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

623 self.beta = beta

624 self._kwds = kwds or {}

625

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

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

628

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

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

631 locations (C{bool}).

632

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

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

635

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

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

638 issue.

639 '''

640 def _xy2(lls, wrap=False):

641 _w = _Wrap._latlonop(wrap)

642 try: # like _xyhs above, but degrees

643 for i, ll in enumerate(lls):

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

645 except AttributeError as x:

646 i = Fmt.SQUARE(llis=i)

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

648

649 _as, llis = _as_llis2(llis)

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

651

652 @property_RO

653 def adjust(self):

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

655 '''

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

657

658 @property

659 def beta(self):

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

661 '''

662 return self._beta

663

664 @beta.setter # PYCHOK setter!

665 def beta(self, beta):

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

667

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

669 '''

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

671

672 @property_RO

673 def datum(self):

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

675 '''

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

677

678 def _datum_setter(self, datum):

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

680 '''

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

682 if d and d is not self._datum:

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

684

685 def _distances(self, x, y):

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

687 '''

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

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

690 notOverloaded(self, distance_function=_f)

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

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

693

694 def _hIDW(self, x, y):

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

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

697 '''

698 try:

699 ds = self._distances(x, y)

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

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

702 except (TypeError, ValueError) as x:

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

704

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

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

707

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

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

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

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

712

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

714 interpolated heights (C{float}s).

715

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

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

718 issue.

719 '''

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

721

722 @property_RO

723 def hypot(self):

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

725 '''

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

727

728 @property_RO

729 def knots(self):

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

731 '''

732 return self._ks

733

734 @property_RO

735 def kwds(self):

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

737 '''

738 return self._kwds

739

740 @property_RO

741 def limit(self):

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

743 '''

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

745

746 @property_RO

747 def radius(self):

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

749 '''

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

751

752 @property_RO

753 def scaled(self):

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

755 '''

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

757

758 @property_RO

759 def wrap(self):

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

761 '''

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

763

764

765 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

769 L{pygeodesy.cosineAndoyerLambert_}.

770 '''

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

772 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

773

774 @kwarg datum_wrap: Optional keyword arguments for function

775 L{pygeodesy.cosineAndoyerLambert}.

776

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

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

779 '''

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

781 self._func = cosineAndoyerLambert

782

783 if _FOR_DOCS:

784 __call__ = _HeightIDW.__call__

785 height = _HeightIDW.height

786

787

788 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

792 L{pygeodesy.cosineForsytheAndoyerLambert_}.

793 '''

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

795 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

796

797 @kwarg datum_wrap: Optional keyword arguments for function

798 L{pygeodesy.cosineAndoyerLambert}.

799

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

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

802 '''

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

804 self._func = cosineForsytheAndoyerLambert

805

806 if _FOR_DOCS:

807 __call__ = _HeightIDW.__call__

808 height = _HeightIDW.height

809

810

811 class HeightIDWcosineLaw(_HeightIDW):

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

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

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

815

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

817 '''

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

819 '''New L{HeightIDWcosineLaw} interpolator.

820

821 @kwarg radius_wrap: Optional keyword arguments for function

822 L{pygeodesy.cosineLaw}.

823

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

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

826 '''

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

828 self._func = cosineLaw

829

830 if _FOR_DOCS:

831 __call__ = _HeightIDW.__call__

832 height = _HeightIDW.height

833

834

835 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

839 conventionally in C{meter}.

840 '''

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

842 '''New L{HeightIDWdistanceTo} interpolator.

843

844 @kwarg distanceTo_kwds: Optional keyword arguments for the

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

846 method.

847

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

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

850

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

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

853 checked however}.

854 '''

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

856 **distanceTo_kwds)

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

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

859 # compatible points in _height_called

860 # instead of class LatLon_ and datum None

861 self._datum = ks[0].datum

862 self._LLis = ks[0].classof

863

864 def _distances(self, x, y):

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

866 '''

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

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

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

870

871 if _FOR_DOCS:

872 __call__ = _HeightIDW.__call__

873 height = _HeightIDW.height

874

875

876 class HeightIDWequirectangular(_HeightIDW):

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

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

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

880 L{pygeodesy.equirectangular_}.

881 '''

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

883 '''New L{HeightIDWequirectangular} interpolator.

884

885 @kwarg adjust_limit_wrap: Optional keyword arguments for

886 function L{pygeodesy.equirectangular_}.

887

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

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

890 '''

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

892 **adjust_limit_wrap)

893

894 def _distances(self, x, y):

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

896 '''

897 _f, kwds = equirectangular_, self._kwds

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

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

900

901 if _FOR_DOCS:

902 __call__ = _HeightIDW.__call__

903 height = _HeightIDW.height

904

905

906 class HeightIDWeuclidean(_HeightIDW):

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

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

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

910 '''

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

912 '''New L{HeightIDWeuclidean} interpolator.

913

914 @kwarg adjust_radius_wrap: Optional keyword arguments for

915 function L{pygeodesy.euclidean_}.

916

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

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

919 '''

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

921 **adjust_radius_wrap)

922 self._func = euclidean

923

924 if _FOR_DOCS:

925 __call__ = _HeightIDW.__call__

926 height = _HeightIDW.height

927

928

929 class HeightIDWexact(_HeightIDW):

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

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

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

933 L{GeodesicExact}C{.Inverse}.

934 '''

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

936 '''New L{HeightIDWexact} interpolator.

937

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

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

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

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

942 of class L{geodesicx.GeodesicExact}.

943

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

945

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

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

948 '''

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

950 self._datum_setter(datum)

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

952

953 if _FOR_DOCS:

954 __call__ = _HeightIDW.__call__

955 height = _HeightIDW.height

956

957

958 class HeightIDWflatLocal(_HeightIDW):

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

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

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

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

963 '''

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

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

966

967 @kwarg datum_hypot_scaled_wrap: Optional keyword arguments for

968 function L{pygeodesy.flatLocal_}.

969

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

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

972 '''

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

974 **datum_hypot_scaled_wrap)

975 self._func = flatLocal

976

977 if _FOR_DOCS:

978 __call__ = _HeightIDW.__call__

979 height = _HeightIDW.height

980

981

982 class HeightIDWflatPolar(_HeightIDW):

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

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

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

986 L{pygeodesy.flatPolar_}.

987 '''

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

989 '''New L{HeightIDWflatPolar} interpolator.

990

991 @kwarg radius_wrap: Optional keyword arguments for function

992 L{pygeodesy.flatPolar}.

993

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

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

996 '''

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

998 self._func = flatPolar

999

1000 if _FOR_DOCS:

1001 __call__ = _HeightIDW.__call__

1002 height = _HeightIDW.height

1003

1004

1005 class HeightIDWhaversine(_HeightIDW):

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

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

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

1009

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

1011 '''

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

1013 '''New L{HeightIDWhaversine} interpolator.

1014

1015 @kwarg radius_wrap: Optional keyword arguments for function

1016 L{pygeodesy.haversine}.

1017

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

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

1020 '''

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

1022 self._func = haversine

1023

1024 if _FOR_DOCS:

1025 __call__ = _HeightIDW.__call__

1026 height = _HeightIDW.height

1027

1028

1029 class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny

1030 if _FOR_DOCS:

1031 __doc__ = HeightIDWflatLocal.__doc__

1032 __init__ = HeightIDWflatLocal.__init__

1033 __call__ = HeightIDWflatLocal.__call__

1034 height = HeightIDWflatLocal.height

1035

1036

1037 class HeightIDWkarney(_HeightIDW):

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

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

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

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

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

1043 Inverse method.

1044 '''

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

1046 '''New L{HeightIDWkarney} interpolator.

1047

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

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

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

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

1052 of class L{geodesicw.Geodesic}.

1053

1054 @raise ImportError: Package U{geographiclib

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

1056

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

1058

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

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

1061 '''

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

1063 self._datum_setter(datum)

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

1065

1066 if _FOR_DOCS:

1067 __call__ = _HeightIDW.__call__

1068 height = _HeightIDW.height

1069

1070

1071 class HeightIDWthomas(_HeightIDW):

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

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

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

1075 '''

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

1077 '''New L{HeightIDWthomas} interpolator.

1078

1079 @kwarg datum_wrap: Optional keyword argument for function

1080 L{pygeodesy.thomas}.

1081

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

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

1084 '''

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

1086 self._func = thomas

1087

1088 if _FOR_DOCS:

1089 __call__ = _HeightIDW.__call__

1090 height = _HeightIDW.height

1091

1092

1093 class HeightIDWvincentys(_HeightIDW):

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

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

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

1097

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

1099 '''

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

1101 '''New L{HeightIDWvincentys} interpolator.

1102

1103 @kwarg radius_wrap: Optional keyword arguments for function

1104 L{pygeodesy.vincentys}.

1105

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

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

1108 '''

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

1110 self._func = vincentys

1111

1112 if _FOR_DOCS:

1113 __call__ = _HeightIDW.__call__

1114 height = _HeightIDW.height

1115

1116

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

1118

1119# **) MIT License

1120#

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

1122#

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

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

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

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

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

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

1129#

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

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

1132#

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

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

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

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

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

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

1139# OTHER DEALINGS IN THE SOFTWARE.