1

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

3

4 u'''Height interpolations at C{LatLon} points from known 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, PI_2, 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, property_ROver

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

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 _orderedup(ts, lo=EPS, hi=PI2-EPS):

173 # clip, order and remove duplicates

174 # p = 0

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

176 # if k > p:

177 # yield k

178 # p = k

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

180

181

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

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

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

185 name, lls = _xkwds_item2(name_lls)

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

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

188 try:

189 for i, ll in enumerate(lls):

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

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

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

193 except Exception as e:

194 i = Fmt.INDEX(name, i)

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

196

197

198 class _HeightBase(_Named): # in .geoids

199 '''(INTERNAL) Interpolator base class.

200 '''

201 _datum = _WGS84 # default

202 _kmin = 2 # min number of knots

203 _LLis = LatLon_ # ._height class

204 _np_sp = None # (numpy, scipy)

205 _wrap = None # wrap knots and llis

206

207 @property_RO

208 def datum(self):

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

210 '''

211 return self._datum

212

213 @property_RO

214 def kmin(self):

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

216 '''

217 return self._kmin

218

219 @property_RO

220 def wrap(self):

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

222 '''

223 return self._wrap

224

225

226 class _HeightsBase(_HeightBase): # in .geoids

227 '''(INTERNAL) Interpolator base class.

228 '''

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

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

231

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

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

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

235 setting.

236

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

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

239

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

241 an invalid B{C{lli}}.

242

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

244

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

246 '''

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

248

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

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

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

252 atype = self.numpy.array

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

254 _as, llis = _as_llis2(llis)

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

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

257

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

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

260 self._notOverloaded(*args)

261

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

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

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

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

266 except Exception as x:

267 raise _SciPyIssue(x)

268

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

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

271

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

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

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

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

276 the B{C{knots}} setting.

277

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

279 interpolated heights (C{float}s).

280

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

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

283

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

285

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

287 '''

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

289

290 def _np_sp2(self, throwarnings=False): # PYCHOK no cover

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

292 '''

293 # raise SciPyWarnings, but not if

294 # scipy has already been imported

295 if throwarnings: # PYCHOK no cover

296 import sys

297 if _scipy_ not in sys.modules:

298 import warnings

299 warnings.filterwarnings(_error_)

300 return self.numpy, self.scipy

301

302 @property_ROver

303 def numpy(self):

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

305 '''

306 return _xnumpy(self.__class__, 1, 9) # overwrite property_ROver

307

308 @property_ROver

309 def scipy(self):

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

311 '''

312 return _xscipy(self.__class__, 1, 2) # overwrite property_ROver

313

314 @property_ROver

315 def scipy_interpolate(self):

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

317 '''

318 _ = self.scipy

319 import scipy.interpolate as spi # scipy 1.2.2

320 return spi # overwrite property_ROver

321

322 def _scipy_version(self, **n):

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

324 '''

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

326

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

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

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

330 n = len(hs)

331 if n < self.kmin:

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

333 if name:

334 self.name = name

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

336

337

338 class HeightCubic(_HeightsBase):

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

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

341 C{kind='cubic'}.

342 '''

343 _interp2d = None

344 _kind = _cubic_

345 _kmin = 16

346

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

348 '''New L{HeightCubic} interpolator.

349

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

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

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

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

354 C{True} (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, **name_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 def __call__(self, *llis, **wrap):

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

377

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

379

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

381 as exception.

382

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

384 '''

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

386

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

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

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

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

391

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

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

394

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

396

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

398 as exception.

399

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

401 '''

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

403

404

405 class HeightLinear(HeightCubic):

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

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

408 C{kind='linear'}.

409 '''

410 _kind = _linear_

411 _kmin = 2

412

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

414 '''New L{HeightLinear} interpolator.

415

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

417 '''

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

419

420 if _FOR_DOCS:

421 __call__ = HeightCubic.__call__

422 height = HeightCubic.height

423

424

425 class HeightLSQBiSpline(_HeightsBase):

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

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

428 interpolate.LSQSphereBivariateSpline.html>}.

429 '''

430 _kmin = 16 # k = 3, always

431

432 def __init__(self, knots, weight=None, low=1e-4, **name_wrap):

433 '''New L{HeightLSQBiSpline} interpolator.

434

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

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

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

438 @kwarg low: Optional lower bound for I{ordered knots} (C{radians}).

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

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

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

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

443

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

445 B{C{knot}}, B{C{weight}} or B{C{eps}}.

446

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

448

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

450 not installed.

451

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

453

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

455 warning as exception.

456 '''

457 np = self.numpy

458 spi = self.scipy_interpolate

459

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

461 n = len(hs)

462

463 w = weight

464 if isscalar(w):

465 w = float(w)

466 if w <= 0:

467 raise HeightError(weight=w)

468 w = [w] * n

469 elif w is not None:

470 m, w = len2(w)

471 if m != n:

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

473 w = map2(float, w)

474 m = min(w)

475 if m <= 0: # PYCHOK no cover

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

477 raise HeightError(w, m)

478 try:

479 if not EPS < low < (PI_2 - EPS): # 1e-4 like SciPy example

480 raise HeightError(low=low)

481 ps = np.array(_orderedup(xs, low, PI2 - low))

482 ts = np.array(_orderedup(ys, low, PI - low))

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

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

485 except Exception as x:

486 raise _SciPyIssue(x)

487

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

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

490

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

492

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

494 warning as exception.

495

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

497 '''

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

499

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

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

502

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

504

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

506 warning as exception.

507

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

509 '''

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

511

512

513 class HeightSmoothBiSpline(_HeightsBase):

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

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

516 interpolate.SmoothSphereBivariateSpline.html>}.

517 '''

518 _kmin = 16 # k = 3, always

519

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

521 '''New L{HeightSmoothBiSpline} interpolator.

522

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

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

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

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

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

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

529

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

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

532

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

534 installed.

535

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

537

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

539 warning as exception.

540 '''

541 spi = self.scipy_interpolate

542

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

544

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

546 try:

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

548 eps=EPS, s=s).ev

549 except Exception as x:

550 raise _SciPyIssue(x)

551

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

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

554

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

556

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

558 warning as exception.

559

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

561 '''

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

563

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

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

566

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

568

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

570 warning as exception.

571

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

573 '''

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

575

576

577 class _HeightIDW(_HeightBase):

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

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

580 interpolators.

581

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

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

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

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

586 classes.

587 '''

588 _beta = 0 # fidw inverse power

589 _func = None # formy function

590 _knots = () # knots list or tuple

591 _kwds = {} # func_ options

592

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

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

595

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

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

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

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

600 retrievable with property C{kwds}.

601

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

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

604 '''

605 name, kwds = _name2__(**name__kwds)

606 if name:

607 self.name = name

608

609 n, self._knots = len2(knots)

610 if n < self.kmin:

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

612 self.beta = beta

613 self._kwds = kwds or {}

614

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

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

617

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

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

620 locations (C{bool}).

621

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

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

624

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

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

627 issue.

628 '''

629 def _xy2(wrap=False):

630 _w = _Wrap._latlonop(wrap)

631 try: # like _xyhs above, but degrees

632 for i, ll in enumerate(llis):

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

634 except Exception as x:

635 i = Fmt.INDEX(llis=i)

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

637

638 _as, llis = _as_llis2(llis)

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

640

641 @property_RO

642 def adjust(self):

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

644 '''

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

646

647 @property

648 def beta(self):

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

650 '''

651 return self._beta

652

653 @beta.setter # PYCHOK setter!

654 def beta(self, beta):

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

656

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

658 '''

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

660

661 @property_RO

662 def datum(self):

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

664 '''

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

666

667 def _datum_setter(self, datum):

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

669 '''

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

671 if d and d is not self._datum:

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

673

674 def _distances(self, x, y):

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

676 '''

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

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

679 self._notOverloaded(distance_function=_f)

680 try:

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

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

683 except Exception as e:

684 i = Fmt.INDEX(knots=i)

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

686

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

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

689

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

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

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

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

694

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

696 interpolated heights (C{float}s).

697

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

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

700 '''

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

702

703 @Property_RO

704 def _heights(self):

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

706 '''

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

708

709 def _hIDW(self, x, y):

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

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

712 '''

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

714 try:

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

716 except (TypeError, ValueError) as e:

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

718

719 @property_RO

720 def hypot(self):

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

722 '''

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

724

725 @property_RO

726 def knots(self):

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

728 '''

729 return self._knots

730

731 @property_RO

732 def kwds(self):

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

734 '''

735 return self._kwds

736

737 @property_RO

738 def limit(self):

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

740 '''

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

742

743 @property_RO

744 def radius(self):

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

746 '''

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

748

749 @property_RO

750 def scaled(self):

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

752 '''

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

754

755 @property_RO

756 def wrap(self):

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

758 '''

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

760

761

762 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

766 L{pygeodesy.cosineAndoyerLambert_}.

767 '''

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

769 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

770

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

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

773 function L{pygeodesy.cosineAndoyerLambert}.

774

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

776 '''

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

778 self._func = _formy.cosineAndoyerLambert

779

780 if _FOR_DOCS:

781 __call__ = _HeightIDW.__call__

782 height = _HeightIDW.height

783

784

785 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

789 L{pygeodesy.cosineForsytheAndoyerLambert_}.

790 '''

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

792 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

793

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

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

796 function L{pygeodesy.cosineForsytheAndoyerLambert}.

797

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

799 '''

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

801 self._func = _formy.cosineForsytheAndoyerLambert

802

803 if _FOR_DOCS:

804 __call__ = _HeightIDW.__call__

805 height = _HeightIDW.height

806

807

808 class HeightIDWcosineLaw(_HeightIDW):

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

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

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

812

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

814 '''

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

816 '''New L{HeightIDWcosineLaw} interpolator.

817

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

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

820 function L{pygeodesy.cosineLaw}.

821

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

823 '''

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

825 self._func = _formy.cosineLaw

826

827 if _FOR_DOCS:

828 __call__ = _HeightIDW.__call__

829 height = _HeightIDW.height

830

831

832 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

836 conventionally in C{meter}.

837 '''

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

839 '''New L{HeightIDWdistanceTo} interpolator.

840

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

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

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

844

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

846

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

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

849 checked}.

850 '''

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

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

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

854 # in _height_called instead of class LatLon_ and datum None

855 self._datum = ks0.datum

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

857

858 def _distances(self, x, y):

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

860 '''

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

862 try:

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

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

865 except Exception as e:

866 i = Fmt.INDEX(knots=i)

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

868

869 if _FOR_DOCS:

870 __call__ = _HeightIDW.__call__

871 height = _HeightIDW.height

872

873

874 class HeightIDWequirectangular(_HeightIDW):

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

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

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

878 L{pygeodesy.equirectangular4}.

879 '''

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

881 '''New L{HeightIDWequirectangular} interpolator.

882

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

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

885 for function L{pygeodesy.equirectangular4}.

886

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

888 '''

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

890

891 def _distances(self, x, y):

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

893 '''

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

895 try:

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

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

898 except Exception as e:

899 i = Fmt.INDEX(knots=i)

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

901

902 if _FOR_DOCS:

903 __call__ = _HeightIDW.__call__

904 height = _HeightIDW.height

905

906

907 class HeightIDWeuclidean(_HeightIDW):

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

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

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

911 '''

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

913 '''New L{HeightIDWeuclidean} interpolator.

914

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

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

917 for function function L{pygeodesy.euclidean}.

918

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

920 '''

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

922 self._func = _formy.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.Inverse}.

934 '''

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

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

943 class L{geodesicx.GeodesicExact}.

944

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

946

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

948 '''

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

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

966

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

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

969 keyword arguments for L{pygeodesy.flatLocal}.

970

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

972 '''

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

974 **name__datum_hypot_scaled_wrap)

975 self._func = _formy.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__radius_wrap):

989 '''New L{HeightIDWflatPolar} interpolator.

990

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

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

993 arguments for function L{pygeodesy.flatPolar}.

994

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

996 '''

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

998 self._func = _formy.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__radius_wrap):

1013 '''New L{HeightIDWhaversine} interpolator.

1014

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

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

1017 arguments for function L{pygeodesy.haversine}.

1018

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

1020 '''

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

1022 self._func = _formy.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 the

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

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

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

1043 geographiclib.geodesic.Geodesic.Inverse>}.

1044 '''

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

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

1053 class L{geodesicw.Geodesic}.

1054

1055 @raise ImportError: Package U{geographiclib

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

1057

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

1059

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

1061 '''

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

1077 '''New L{HeightIDWthomas} interpolator.

1078

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

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

1081 arguments for function L{pygeodesy.thomas}.

1082

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

1084 '''

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

1086 self._func = _formy.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__radius_wrap):

1101 '''New L{HeightIDWvincentys} interpolator.

1102

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

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

1105 arguments for function L{pygeodesy.vincentys}.

1106

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

1108 '''

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

1110 self._func = _formy.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-2024 -- 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.