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 # _MODS

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

80 _STAR_, _version2

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

82 from pygeodesy.named import _Named, notOverloaded

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

84 from pygeodesy.props import Property_RO, property_RO

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

86 from pygeodesy.units import _isDegrees, Float_, Int_

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

88

89# from math import radians # from .formy

90

91__all__ = _ALL_LAZY.heights

92__version__ = '23.12.03'

93

94_error_ = 'error'

95_llis_ = 'llis'

96_smoothing_ = 'smoothing'

97

98

99 class HeightError(PointsError):

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

101 '''

102 pass

103

104

105 def _alist(ais):

106 # return list of floats, not numpy.float64s

107 return list(map(float, ais))

108

109

110 def _ascalar(ais): # in .geoids

111 # return single float, not numpy.float64

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

113 if len(ais) != 1:

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

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

116 raise _AssertionError(n, txt=t)

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

118

119

120 def _atuple(ais):

121 # return tuple of floats, not numpy.float64s

122 return tuple(map(float, ais))

123

124

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

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

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

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

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

130

131 n = len(llis)

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

133 llis = llis[0]

134 try:

135 n, llis = len2(llis)

136 _as = _alist # return list of interpolated heights

137 except TypeError: # single lli

138 n, llis = 1, [llis]

139 _as = _ascalar # return single interpolated heights

140 else: # of 0, 2 or more llis

141 _as = _atuple # return tuple of interpolated heights

142

143 if n < m:

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

145 return _as, llis

146

147

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

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

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

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

152 else:

153 n, lats = len2(lats)

154 m, lons = len2(lons)

155 if n != m:

156 # format a LenError, but raise an Error

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

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

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

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

161

162

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

164 # create an insufficient Error instance

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

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

167

168

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

170 # clip, order and remove duplicates

171 # p, ks = 0, []

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

173 # if k > p:

174 # ks.append(k)

175 # p = k

176 # return ks

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

178

179

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

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

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

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

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

185 try:

186 for i, ll in enumerate(lls):

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

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

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

190 except AttributeError as x:

191 i = Fmt.SQUARE(name, i)

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

193

194

195 class _HeightBase(_Named): # in .geoids

196 '''(INTERNAL) Interpolator base class.

197 '''

198 _datum = _WGS84 # default

199 _kmin = 2 # min number of knots

200 _LLis = LatLon_ # ._height class

201 _np_sp = None # (numpy, scipy)

202 _wrap = None # wrap knots and llis

203

204 @property_RO

205 def datum(self):

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

207 '''

208 return self._datum

209

210 @property_RO

211 def kmin(self):

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

213 '''

214 return self._kmin

215

216 @property_RO

217 def wrap(self):

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

219 '''

220 return self._wrap

221

222

223 class _HeightsBase(_HeightBase): # in .geoids

224 '''(INTERNAL) Interpolator base class.

225 '''

226 _np_sp = None # (numpy, scipy)

227

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

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

230

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

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

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. I{Must be overloaded}.

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 _scipy_version(self, **n):

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

332 '''

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

334

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

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

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

338 n = len(hs)

339 if n < self.kmin:

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

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

342

343

344 class HeightCubic(_HeightsBase):

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

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

347 C{kind='cubic'}.

348 '''

349 _interp2d = None

350 _kind = _cubic_

351 _kmin = 16

352

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

354 '''New L{HeightCubic} interpolator.

355

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

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

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

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

360

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

362 invalid B{C{knot}}.

363

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

365 or not installed.

366

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

368

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

370 as exception.

371 '''

372 spi = self.scipy_interpolate

373

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

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

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

377 except Exception as x:

378 raise _SciPyIssue(x)

379

380 if name:

381 self.name = name

382

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

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

385

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

387

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

389 as exception.

390

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

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

393 '''

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

395

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

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

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

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

400

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

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

403

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

405

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

407 as exception.

408

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

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

411 '''

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

413

414

415 class HeightLinear(HeightCubic):

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

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

418 C{kind='linear'}.

419 '''

420 _kind = _linear_

421 _kmin = 2

422

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

424 '''New L{HeightLinear} interpolator.

425

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

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

428 '''

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

430

431 if _FOR_DOCS:

432 __call__ = HeightCubic.__call__

433 height = HeightCubic.height

434

435

436 class HeightLSQBiSpline(_HeightsBase):

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

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

439 interpolate.LSQSphereBivariateSpline.html>}.

440 '''

441 _kmin = 16 # k = 3, always

442

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

444 '''New L{HeightLSQBiSpline} interpolator.

445

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

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

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

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

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

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

452

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

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

455

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

457

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

459 or not installed.

460

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

462

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

464 warning as exception.

465 '''

466 np = self.numpy

467 spi = self.scipy_interpolate

468

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

470 n = len(hs)

471

472 w = weight

473 if isscalar(w):

474 w = float(w)

475 if w <= 0:

476 raise HeightError(weight=w)

477 w = [w] * n

478 elif w is not None:

479 m, w = len2(w)

480 if m != n:

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

482 w = map2(float, w)

483 m = min(w)

484 if m <= 0: # PYCHOK no cover

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

486 raise HeightError(w, m)

487 try:

488 T = 1.0e-4 # like SciPy example

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

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

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

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

493 except Exception as x:

494 raise _SciPyIssue(x)

495

496 if name:

497 self.name = name

498

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

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

501

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

503

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

505 warning as exception.

506

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

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

509 '''

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

511

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

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

514

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

516

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

518 warning as exception.

519

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

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

522 '''

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

524

525

526 class HeightSmoothBiSpline(_HeightsBase):

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

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

529 interpolate.SmoothSphereBivariateSpline.html>}.

530 '''

531 _kmin = 16 # k = 3, always

532

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

534 '''New L{HeightSmoothBiSpline} interpolator.

535

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

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

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

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

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

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

542

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

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

545

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

547 or not installed.

548

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

550

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

552 warning as exception.

553 '''

554 spi = self.scipy_interpolate

555

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

557

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

559 try:

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

561 eps=EPS, s=s).ev

562 except Exception as x:

563 raise _SciPyIssue(x)

564

565 if name:

566 self.name = name

567

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

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

570

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

572

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

574 warning as exception.

575

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

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

578 '''

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

580

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

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

583

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

585

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

587 warning as exception.

588

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

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

591 '''

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

593

594

595 class _HeightIDW(_HeightBase):

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

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

598 interpolators.

599

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

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

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

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

604 classes.

605 '''

606 _beta = 0 # fidw inverse power

607 _func = None # formy function

608 _ks = () # knots list or tuple

609 _kwds = {} # func_ options

610

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

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

613

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

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

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

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

618 retrievable with property C{kwds}.

619

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

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

622 '''

623 if name:

624 self.name = name

625 n, self._ks = len2(knots)

626 if n < self.kmin:

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

628 self.beta = beta

629 self._kwds = kwds or {}

630

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

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

633

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

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

636 locations (C{bool}).

637

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

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

640

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

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

643 issue.

644 '''

645 def _xy2(lls, wrap=False):

646 _w = _Wrap._latlonop(wrap)

647 try: # like _xyhs above, but degrees

648 for i, ll in enumerate(lls):

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

650 except AttributeError as x:

651 i = Fmt.SQUARE(llis=i)

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

653

654 _as, llis = _as_llis2(llis)

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

656

657 @property_RO

658 def adjust(self):

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

660 '''

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

662

663 @property

664 def beta(self):

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

666 '''

667 return self._beta

668

669 @beta.setter # PYCHOK setter!

670 def beta(self, beta):

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

672

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

674 '''

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

676

677 @property_RO

678 def datum(self):

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

680 '''

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

682

683 def _datum_setter(self, datum):

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

685 '''

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

687 if d and d is not self._datum:

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

689

690 def _distances(self, x, y):

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

692 '''

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

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

695 notOverloaded(self, distance_function=_f)

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

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

698

699 def _hIDW(self, x, y):

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

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

702 '''

703 try:

704 ds = self._distances(x, y)

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

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

707 except (TypeError, ValueError) as x:

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

709

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

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

712

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

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

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

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

717

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

719 interpolated heights (C{float}s).

720

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

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

723 issue.

724 '''

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

726

727 @property_RO

728 def hypot(self):

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

730 '''

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

732

733 @property_RO

734 def knots(self):

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

736 '''

737 return self._ks

738

739 @property_RO

740 def kwds(self):

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

742 '''

743 return self._kwds

744

745 @property_RO

746 def limit(self):

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

748 '''

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

750

751 @property_RO

752 def radius(self):

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

754 '''

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

756

757 @property_RO

758 def scaled(self):

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

760 '''

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

762

763 @property_RO

764 def wrap(self):

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

766 '''

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

768

769

770 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

774 L{pygeodesy.cosineAndoyerLambert_}.

775 '''

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

777 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

778

779 @kwarg datum_wrap: Optional keyword arguments for function

780 L{pygeodesy.cosineAndoyerLambert}.

781

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

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

784 '''

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

786 self._func = cosineAndoyerLambert

787

788 if _FOR_DOCS:

789 __call__ = _HeightIDW.__call__

790 height = _HeightIDW.height

791

792

793 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

797 L{pygeodesy.cosineForsytheAndoyerLambert_}.

798 '''

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

800 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

801

802 @kwarg datum_wrap: Optional keyword arguments for function

803 L{pygeodesy.cosineAndoyerLambert}.

804

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

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

807 '''

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

809 self._func = cosineForsytheAndoyerLambert

810

811 if _FOR_DOCS:

812 __call__ = _HeightIDW.__call__

813 height = _HeightIDW.height

814

815

816 class HeightIDWcosineLaw(_HeightIDW):

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

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

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

820

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

822 '''

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

824 '''New L{HeightIDWcosineLaw} interpolator.

825

826 @kwarg radius_wrap: Optional keyword arguments for function

827 L{pygeodesy.cosineLaw}.

828

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

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

831 '''

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

833 self._func = cosineLaw

834

835 if _FOR_DOCS:

836 __call__ = _HeightIDW.__call__

837 height = _HeightIDW.height

838

839

840 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

844 conventionally in C{meter}.

845 '''

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

847 '''New L{HeightIDWdistanceTo} interpolator.

848

849 @kwarg distanceTo_kwds: Optional keyword arguments for the

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

851 method.

852

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

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

855

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

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

858 checked however}.

859 '''

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

861 **distanceTo_kwds)

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

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

864 # compatible points in _height_called

865 # instead of class LatLon_ and datum None

866 self._datum = ks[0].datum

867 self._LLis = ks[0].classof

868

869 def _distances(self, x, y):

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

871 '''

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

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

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

875

876 if _FOR_DOCS:

877 __call__ = _HeightIDW.__call__

878 height = _HeightIDW.height

879

880

881 class HeightIDWequirectangular(_HeightIDW):

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

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

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

885 L{pygeodesy.equirectangular_}.

886 '''

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

888 '''New L{HeightIDWequirectangular} interpolator.

889

890 @kwarg adjust_limit_wrap: Optional keyword arguments for

891 function L{pygeodesy.equirectangular_}.

892

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

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

895 '''

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

897 **adjust_limit_wrap)

898

899 def _distances(self, x, y):

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

901 '''

902 _f, kwds = equirectangular_, self._kwds

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

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

905

906 if _FOR_DOCS:

907 __call__ = _HeightIDW.__call__

908 height = _HeightIDW.height

909

910

911 class HeightIDWeuclidean(_HeightIDW):

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

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

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

915 '''

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

917 '''New L{HeightIDWeuclidean} interpolator.

918

919 @kwarg adjust_radius_wrap: Optional keyword arguments for

920 function L{pygeodesy.euclidean_}.

921

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

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

924 '''

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

926 **adjust_radius_wrap)

927 self._func = euclidean

928

929 if _FOR_DOCS:

930 __call__ = _HeightIDW.__call__

931 height = _HeightIDW.height

932

933

934 class HeightIDWexact(_HeightIDW):

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

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

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

938 L{GeodesicExact}C{.Inverse}.

939 '''

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

941 '''New L{HeightIDWexact} interpolator.

942

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

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

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

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

947 of class L{geodesicx.GeodesicExact}.

948

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

950

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

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

953 '''

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

955 self._datum_setter(datum)

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

957

958 if _FOR_DOCS:

959 __call__ = _HeightIDW.__call__

960 height = _HeightIDW.height

961

962

963 class HeightIDWflatLocal(_HeightIDW):

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

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

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

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

968 '''

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

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

971

972 @kwarg datum_hypot_scaled_wrap: Optional keyword arguments for

973 function L{pygeodesy.flatLocal_}.

974

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

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

977 '''

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

979 **datum_hypot_scaled_wrap)

980 self._func = flatLocal

981

982 if _FOR_DOCS:

983 __call__ = _HeightIDW.__call__

984 height = _HeightIDW.height

985

986

987 class HeightIDWflatPolar(_HeightIDW):

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

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

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

991 L{pygeodesy.flatPolar_}.

992 '''

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

994 '''New L{HeightIDWflatPolar} interpolator.

995

996 @kwarg radius_wrap: Optional keyword arguments for function

997 L{pygeodesy.flatPolar}.

998

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

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

1001 '''

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

1003 self._func = flatPolar

1004

1005 if _FOR_DOCS:

1006 __call__ = _HeightIDW.__call__

1007 height = _HeightIDW.height

1008

1009

1010 class HeightIDWhaversine(_HeightIDW):

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

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

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

1014

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

1016 '''

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

1018 '''New L{HeightIDWhaversine} interpolator.

1019

1020 @kwarg radius_wrap: Optional keyword arguments for function

1021 L{pygeodesy.haversine}.

1022

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

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

1025 '''

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

1027 self._func = haversine

1028

1029 if _FOR_DOCS:

1030 __call__ = _HeightIDW.__call__

1031 height = _HeightIDW.height

1032

1033

1034 class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny

1035 if _FOR_DOCS:

1036 __doc__ = HeightIDWflatLocal.__doc__

1037 __init__ = HeightIDWflatLocal.__init__

1038 __call__ = HeightIDWflatLocal.__call__

1039 height = HeightIDWflatLocal.height

1040

1041

1042 class HeightIDWkarney(_HeightIDW):

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

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

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

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

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

1048 Inverse method.

1049 '''

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

1051 '''New L{HeightIDWkarney} interpolator.

1052

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

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

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

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

1057 of class L{geodesicw.Geodesic}.

1058

1059 @raise ImportError: Package U{geographiclib

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

1061

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

1063

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

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

1066 '''

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

1068 self._datum_setter(datum)

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

1070

1071 if _FOR_DOCS:

1072 __call__ = _HeightIDW.__call__

1073 height = _HeightIDW.height

1074

1075

1076 class HeightIDWthomas(_HeightIDW):

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

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

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

1080 '''

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

1082 '''New L{HeightIDWthomas} interpolator.

1083

1084 @kwarg datum_wrap: Optional keyword argument for function

1085 L{pygeodesy.thomas}.

1086

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

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

1089 '''

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

1091 self._func = thomas

1092

1093 if _FOR_DOCS:

1094 __call__ = _HeightIDW.__call__

1095 height = _HeightIDW.height

1096

1097

1098 class HeightIDWvincentys(_HeightIDW):

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

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

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

1102

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

1104 '''

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

1106 '''New L{HeightIDWvincentys} interpolator.

1107

1108 @kwarg radius_wrap: Optional keyword arguments for function

1109 L{pygeodesy.vincentys}.

1110

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

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

1113 '''

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

1115 self._func = vincentys

1116

1117 if _FOR_DOCS:

1118 __call__ = _HeightIDW.__call__

1119 height = _HeightIDW.height

1120

1121

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

1123

1124# **) MIT License

1125#

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

1127#

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

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

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

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

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

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

1134#

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

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

1137#

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

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

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

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

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

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

1144# OTHER DEALINGS IN THE SOFTWARE.