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, PI_2, _0_0, _90_0, _180_0

71 from pygeodesy.datums import _ellipsoidal_datum, _WGS84

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

73 from pygeodesy.fmath import fidw, hypot2

74 from pygeodesy.formy import cosineAndoyerLambert_, cosineForsytheAndoyerLambert_, \

75 cosineLaw_, euclidean_, flatPolar_, haversine_, \

76 _scale_rad, thomas_, vincentys_

77 from pygeodesy.interns import NN, _COMMASPACE_, _cubic_, _datum_, _distanceTo_, \

78 _knots_, _linear_, _NOTEQUAL_, _scipy_, _SPACE_, _STAR_

79 from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS

80 from pygeodesy.named import _Named, notOverloaded

81 from pygeodesy.points import LatLon_, Property_RO

82# from pygeodesy.props import Property_RO # from .points

83 from pygeodesy.streprs import _boolkwds, Fmt

84 from pygeodesy.units import Float_, Int_

85 from pygeodesy.utily import radiansPI, radiansPI2, unrollPI

86

87__all__ = _ALL_LAZY.heights

88__version__ = '22.10.06'

89

90_error_ = 'error'

91_insufficient_ = 'insufficient'

92_llis_ = 'llis'

93_smoothing_ = 'smoothing'

94

95

96 class HeightError(PointsError):

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

98 '''

99 pass

100

101

102 def _alist(ais):

103 # return list of floats, not numpy.float64s

104 return list(map(float, ais))

105

106

107 def _allis2(llis, m=1, Error=HeightError): # imported by .geoids

108 # dtermine return type and convert lli C{LatLon}s to list

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

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

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

112

113 n = len(llis)

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

115 llis = llis[0]

116 try:

117 n, llis = len2(llis)

118 _as = _alist # return list of interpolated heights

119 except TypeError: # single lli

120 n, llis = 1, [llis]

121 _as = _ascalar # return single interpolated heights

122 else: # of 0, 2 or more llis

123 _as = _atuple # return tuple of interpolated heights

124

125 if n < m:

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

127 return _as, llis

128

129

130 def _ascalar(ais): # imported by .geoids

131 # return single float, not numpy.float64

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

133 if len(ais) != 1:

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

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

136 raise _AssertionError(n, txt=t)

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

138

139

140 def _atuple(ais):

141 # return tuple of floats, not numpy.float64s

142 return tuple(map(float, ais))

143

144

145 def _axyllis4(atype, llis, m=1, off=True):

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

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

148 _as, llis = _allis2(llis, m=m)

149 xis, yis, _ = zip(*_xyhs(llis, off=off)) # PYCHOK unzip

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

151

152

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

154 # create an insufficient Error instance

155 t = _COMMASPACE_(_insufficient_, str(need))

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

157

158

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

160 # clip, order and remove duplicates

161 # p, ks = 0, []

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

163 # if k > p:

164 # ks.append(k)

165 # p = k

166 # return ks

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

168

169

170 def _xyhs(lls, off=True, name=_llis_):

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

172 # x: 0 <= lon <= PI2, y: 0 <= lat <= PI if off is True

173 # else x: -PI <= lon <= PI, y: -PI_2 <= lat <= PI_2

174 if off:

175 xf = yf = _0_0

176 else: # undo offset

177 xf, yf = PI, PI_2

178 try:

179 for i, ll in enumerate(lls):

180 yield (max(_0_0, radiansPI2(ll.lon + _180_0)) - xf), \

181 (max(_0_0, radiansPI( ll.lat + _90_0)) - yf), ll.height

182 except AttributeError as x:

183 i = Fmt.SQUARE(name, i)

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

185

186

187 def _xyhs3(atype, m, knots, off=True):

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

189 xs, ys, hs = zip(*_xyhs(knots, off=off, name=_knots_)) # PYCHOK unzip

190 n = len(hs)

191 if n < m:

192 raise _insufficientError(m, knots=n)

193 return map1(atype, xs, ys, hs)

194

195

196 class _HeightBase(_Named): # imported by .geoids

197 '''(INTERNAL) Interpolator base class.

198 '''

199 _adjust = None # not applicable

200 _datum = None # ._height datum

201 _kmin = 2 # min number of knots

202 _LLis = LatLon_ # ._height class

203 _np_sp = None # (numpy, scipy)

204 _wrap = None # not applicable

205

206 def __call__(self, *args): # PYCHOK no cover

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

208 '''

209 notOverloaded(self, callername='__call__', *args)

210

211 @Property_RO

212 def adjust(self):

213 '''Get the adjust setting (C{bool} or C{None} if not applicable).

214 '''

215 return self._adjust

216

217 def _axyllis4(self, llis):

218 return _axyllis4(self.numpy.array, llis)

219

220 @Property_RO

221 def datum(self):

222 '''Get the datum (L{Datum} or C{None} if not applicable).

223 '''

224 return self._datum

225

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

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

228 '''

229 notOverloaded(self, *args)

230

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

232 _as, xis, yis, _ = self._axyllis4(llis)

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

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

235 except Exception as x:

236 raise _SciPyIssue(x)

237

238 def _height(self, lats, lons, Error=HeightError):

239 LLis, d = self._LLis, self.datum

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

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

242 else:

243 n, lats = len2(lats)

244 m, lons = len2(lons)

245 if n != m:

246 # format a LenError, but raise an Error

247 e = LenError(self.__class__, lats=n, lons=m, txt=None)

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

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

250 return self(llis) # __call__(lli) or __call__(llis)

251

252 @Property_RO

253 def kmin(self):

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

255 '''

256 return self._kmin

257

258 def _np_sp2(self, throwarnings=False):

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

260 '''

261 t = _HeightBase._np_sp

262 if not t:

263 # raise SciPyWarnings, but not if

264 # scipy has been imported already

265 if throwarnings: # PYCHOK no cover

266 import sys

267 if _scipy_ not in sys.modules:

268 import warnings

269 warnings.filterwarnings(_error_)

270

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

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

273

274 _HeightBase._np_sp = t = np, sp

275 return t

276

277 @Property_RO

278 def numpy(self):

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

280 '''

281 np, _ = self._np_sp2()

282 return np

283

284 @Property_RO

285 def scipy(self):

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

287 '''

288 _, sp = self._np_sp2()

289 return sp

290

291 @Property_RO

292 def scipy_interpolate(self):

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

294 '''

295 _ = self.scipy

296 import scipy.interpolate as spi

297 return spi

298

299 def _xyhs3(self, knots):

300 return _xyhs3(self.numpy.array, self._kmin, knots)

301

302 @Property_RO

303 def wrap(self):

304 '''Get the wrap setting (C{bool} or C{None} if not applicable).

305 '''

306 return self._wrap

307

308

309 class HeightCubic(_HeightBase):

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

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

312 C{kind='cubic'}.

313 '''

314 _interp2d = None

315 _kind = _cubic_

316 _kmin = 16

317

318 def __init__(self, knots, name=NN):

319 '''New L{HeightCubic} interpolator.

320

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

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

323

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

325 invalid B{C{knot}}.

326

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

328 or not installed.

329

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

331

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

333 as exception.

334 '''

335 spi = self.scipy_interpolate

336

337 xs, ys, hs = self._xyhs3(knots)

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

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

340 except Exception as x:

341 raise _SciPyIssue(x)

342

343 if name:

344 self.name = name

345

346 def __call__(self, *llis):

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

348

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

350 C{LatLon}s).

351

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

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

354

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

356 an invalid B{C{lli}}.

357

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

359

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

361 as exception.

362 '''

363 return _HeightBase._eval(self, llis)

364

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

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

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

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

369

370 def height(self, lats, lons):

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

372

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

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

375

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

377 interpolated heights (C{float}s).

378

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

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

381

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

383

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

385 as exception.

386 '''

387 return _HeightBase._height(self, lats, lons)

388

389

390 class HeightLinear(HeightCubic):

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

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

393 C{kind='linear'}.

394 '''

395 _kind = _linear_

396 _kmin = 2

397

398 def __init__(self, knots, name=NN):

399 '''New L{HeightLinear} interpolator.

400

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

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

403

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

405 an invalid B{C{knot}}.

406

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

408 or not installed.

409

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

411

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

413 as exception.

414 '''

415 HeightCubic.__init__(self, knots, name=name)

416

417 if _FOR_DOCS:

418 __call__ = HeightCubic.__call__

419 height = HeightCubic.height

420

421

422 class _HeightIDW(_HeightBase):

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

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

425 height interpolators.

426 '''

427 _beta = 0 # inverse distance power

428 _hs = () # known heights

429 _xs = () # knot lons

430 _ys = () # knot lats

431

432 def __init__(self, knots, beta=2, name=NN, **wrap_adjust):

433 '''New L{_HeightIDW} interpolator.

434 '''

435 self._xs, self._ys, self._hs = _xyhs3(tuple, self._kmin, knots, off=False)

436 self.beta = beta

437 if name:

438 self.name = name

439 if wrap_adjust:

440 _boolkwds(self, **wrap_adjust)

441

442 def __call__(self, *llis):

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

444

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

446 C{LatLon}s).

447

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

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

450

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

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

453 issue.

454 '''

455 _as, xis, yis, _ = _axyllis4(tuple, llis, off=False)

456 return _as(map(self._hIDW, xis, yis))

457

458 def _datum_setter(self, datum, knots):

459 '''(INTERNAL) Set the datum.

460

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

462 '''

463 d = datum or getattr(knots[0], _datum_, datum)

464 if d not in (None, self._datum):

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

466

467 def _distances(self, x, y): # PYCHOK unused (x, y) radians

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

469 '''

470 notOverloaded(self, x, y)

471

472 def _distances_angular_(self, func_, x, y):

473 # return angular distances from func_

474 for xk, yk in zip(self._xs, self._ys):

475 r, _ = unrollPI(xk, x, wrap=self._wrap)

476 yield func_(yk, y, r)

477

478 def _distances_angular_datum_(self, func_, x, y):

479 # return angular distances from func_

480 for xk, yk in zip(self._xs, self._ys):

481 r, _ = unrollPI(xk, x, wrap=self._wrap)

482 yield func_(yk, y, r, datum=self._datum)

483

484 def _hIDW(self, x, y):

485 # interpolate height at (x, y) radians or degrees

486 try:

487 ds = self._distances(x, y)

488 return fidw(self._hs, ds, beta=self._beta)

489 except (TypeError, ValueError) as x:

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

491

492 @property

493 def beta(self):

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

495 '''

496 return self._beta

497

498 @beta.setter # PYCHOK setter!

499 def beta(self, beta):

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

501

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

503 '''

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

505

506 def height(self, lats, lons):

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

508

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

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

511

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

513 interpolated heights (C{float}s).

514

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

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

517 issue.

518 '''

519 return _HeightBase._height(self, lats, lons)

520

521

522 class HeightIDWcosineAndoyerLambert(_HeightIDW):

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

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

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

526 L{pygeodesy.cosineAndoyerLambert_}.

527

528 @see: L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWdistanceTo},

529 L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas},

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

531 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

533 shepard_interp_2d/shepard_interp_2d.html>}.

534 '''

535 _datum = _WGS84

536 _wrap = False

537

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

539 '''New L{HeightIDWcosineAndoyerLambert} interpolator.

540

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

542 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}

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

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

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

546 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

548

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

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

551

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

553 '''

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

555 self._datum_setter(datum, knots)

556

557 def _distances(self, x, y): # (x, y) radians

558 return self._distances_angular_datum_(cosineAndoyerLambert_, x, y)

559

560 if _FOR_DOCS:

561 __call__ = _HeightIDW.__call__

562 height = _HeightIDW.height

563

564

565 class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):

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

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

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

569 L{pygeodesy.cosineForsytheAndoyerLambert_}.

570

571 @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWdistanceTo},

572 L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas},

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

574 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

576 shepard_interp_2d/shepard_interp_2d.html>}.

577 '''

578 _datum = _WGS84

579 _wrap = False

580

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

582 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.

583

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

585 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}

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

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

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

589 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

591

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

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

594

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

596 '''

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

598 self._datum_setter(datum, knots)

599

600 def _distances(self, x, y): # (x, y) radians

601 return self._distances_angular_datum_(cosineForsytheAndoyerLambert_, x, y)

602

603 if _FOR_DOCS:

604 __call__ = _HeightIDW.__call__

605 height = _HeightIDW.height

606

607

608 class HeightIDWcosineLaw(_HeightIDW):

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

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

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

612

613 @see: L{HeightIDWequirectangular}, L{HeightIDWeuclidean},

614 L{HeightIDWflatPolar}, L{HeightIDWhaversine}, L{HeightIDWvincentys},

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

616 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

618 shepard_interp_2d/shepard_interp_2d.html>}.

619

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

621 '''

622 _wrap = False

623

624 def __init__(self, knots, beta=2, wrap=False, name=NN):

625 '''New L{HeightIDWcosineLaw} interpolator.

626

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

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

629 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

631

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

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

634 '''

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

636

637 def _distances(self, x, y): # (x, y) radians

638 return self._distances_angular_(cosineLaw_, x, y)

639

640 if _FOR_DOCS:

641 __call__ = _HeightIDW.__call__

642 height = _HeightIDW.height

643

644

645 class HeightIDWdistanceTo(_HeightIDW):

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

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

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

649 conventionally in C{meter}.

650

651 @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert},

652 L{HeightIDWflatPolar}, L{HeightIDWkarney}, L{HeightIDWthomas},

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

654 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

656 shepard_interp_2d/shepard_interp_2d.html>}.

657 '''

658 _distanceTo_kwds = {}

659 _ks = ()

660

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

662 '''New L{HeightIDWdistanceTo} interpolator.

663

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

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

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

667 @kwarg distanceTo_kwds: Optional keyword arguments for the

668 B{C{points}}' C{LatLon.distanceTo}

669 method.

670

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

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

673

674 @raise ImportError: Package U{geographiclib

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

676 iff B{C{points}} are L{ellipsoidalKarney.LatLon}s.

677

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

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

680 checked however}.

681 '''

682 n, self._ks = len2(knots)

683 if n < self._kmin:

684 raise _insufficientError(self._kmin, knots=n)

685 for i, k in enumerate(self._ks):

686 if not callable(getattr(k, _distanceTo_, None)):

687 i = Fmt.SQUARE(_knots_, i)

688 raise HeightError(i, k, txt=_distanceTo_)

689

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

691 # compatible points in _HeightBase._height

692 # instead of class LatLon_ and datum None

693 self._datum = self._ks[0].datum

694 self._LLis = self._ks[0].classof

695

696 self.beta = beta

697 if name:

698 self.name = name

699 if distanceTo_kwds:

700 self._distanceTo_kwds = distanceTo_kwds

701

702 def __call__(self, *llis):

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

704

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

706 C{LatLon}s).

707

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

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

710

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

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

713 issue.

714 '''

715 _as, llis = _allis2(llis)

716 return _as(map(self._hIDW, llis))

717

718 if _FOR_DOCS:

719 height = _HeightIDW.height

720

721 def _hIDW(self, lli): # PYCHOK expected

722 # interpolate height at point lli

723 try:

724 kwds = self._distanceTo_kwds

725 ds = (k.distanceTo(lli, **kwds) for k in self._ks)

726 return fidw(self._hs, ds, beta=self._beta)

727 except (TypeError, ValueError) as x:

728 raise HeightError(str(x))

729

730 @property # NOT _RO, no caching

731 def _hs(self): # see HeightIDWkarney

732 for k in self._ks:

733 yield k.height

734

735

736 class HeightIDWequirectangular(_HeightIDW):

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

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

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

740 L{pygeodesy.equirectangular_}.

741

742 @see: L{HeightIDWeuclidean}, L{HeightIDWflatPolar},

743 L{HeightIDWhaversine}, L{HeightIDWvincentys},

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

745 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

747 shepard_interp_2d/shepard_interp_2d.html>}.

748 '''

749 _adjust = True

750 _wrap = False

751

752 def __init__(self, knots, adjust=True, wrap=False, name=NN):

753 '''New L{HeightIDWequirectangular} interpolator.

754

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

756 @kwarg adjust: Adjust the wrapped, unrolled longitudinal

757 delta by the cosine of the mean latitude (C{bool}).

758 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

760

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

762 an invalid B{C{knot}}.

763 '''

764 _HeightIDW.__init__(self, knots, beta=1, name=name, wrap=wrap,

765 adjust=adjust)

766

767 def _distances(self, x, y): # (x, y) radians**2

768 for xk, yk in zip(self._xs, self._ys):

769 d, _ = unrollPI(xk, x, wrap=self._wrap)

770 if self._adjust:

771 d *= _scale_rad(yk, y)

772 yield hypot2(d, yk - y) # like equirectangular_ distance2

773

774 if _FOR_DOCS:

775 __call__ = _HeightIDW.__call__

776 height = _HeightIDW.height

777

778

779 class HeightIDWeuclidean(_HeightIDW):

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

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

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

783

784 @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular},

785 L{HeightIDWhaversine}, L{HeightIDWvincentys},

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

787 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

789 shepard_interp_2d/shepard_interp_2d.html>}.

790 '''

791 _adjust = True

792

793 def __init__(self, knots, adjust=True, beta=2, name=NN):

794 '''New L{HeightIDWeuclidean} interpolator.

795

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

797 @kwarg adjust: Adjust the longitudinal delta by the cosine

798 of the mean latitude for B{C{adjust}}=C{True}.

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

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

801

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

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

804 '''

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

806

807 def _distances(self, x, y): # (x, y) radians

808 for xk, yk in zip(self._xs, self._ys):

809 yield euclidean_(yk, y, xk - x, adjust=self._adjust)

810

811 if _FOR_DOCS:

812 __call__ = _HeightIDW.__call__

813 height = _HeightIDW.height

814

815

816 class HeightIDWflatLocal(_HeightIDW):

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

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

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

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

821

822 @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert},

823 L{HeightIDWdistanceTo}, L{HeightIDWhubeny}, L{HeightIDWthomas},

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

825 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

827 shepard_interp_2d/shepard_interp_2d.html>}.

828 '''

829 _datum = _WGS84

830 _wrap = False

831

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

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

834

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

836 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}

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

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

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

840 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

842

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

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

845

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

847 '''

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

849 self._datum_setter(datum, knots)

850

851 def _distances(self, x, y): # (x, y) radians

852 return self._distances_angular_(self._hubeny_2, x, y) # radians**2

853

854 @Property_RO

855 def _hubeny_2(self):

856 '''(INTERNAL) Get and cache the C{.datum.ellipsoid._hubeny_2} method.

857 '''

858 return self.datum.ellipsoid._hubeny_2

859

860 if _FOR_DOCS:

861 __call__ = _HeightIDW.__call__

862 height = _HeightIDW.height

863

864

865 class HeightIDWflatPolar(_HeightIDW):

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

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

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

869 L{pygeodesy.flatPolar_}.

870

871 @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular},

872 L{HeightIDWeuclidean}, L{HeightIDWhaversine}, L{HeightIDWvincentys},

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

874 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

876 shepard_interp_2d/shepard_interp_2d.html>}.

877 '''

878 _wrap = False

879

880 def __init__(self, knots, beta=2, wrap=False, name=NN):

881 '''New L{HeightIDWflatPolar} interpolator.

882

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

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

885 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

887

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

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

890 '''

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

892

893 def _distances(self, x, y): # (x, y) radians

894 return self._distances_angular_(flatPolar_, x, y)

895

896 if _FOR_DOCS:

897 __call__ = _HeightIDW.__call__

898 height = _HeightIDW.height

899

900

901 class HeightIDWhaversine(_HeightIDW):

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

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

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

905

906 @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular}, L{HeightIDWeuclidean},

907 L{HeightIDWflatPolar}, L{HeightIDWvincentys},

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

909 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

911 shepard_interp_2d/shepard_interp_2d.html>}.

912

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

914 '''

915 _wrap = False

916

917 def __init__(self, knots, beta=2, wrap=False, name=NN):

918 '''New L{HeightIDWhaversine} interpolator.

919

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

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

922 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

924

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

926 an B{C{knot}} or B{C{beta}}.

927 '''

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

929

930 def _distances(self, x, y): # (x, y) radians

931 return self._distances_angular_(haversine_, x, y)

932

933 if _FOR_DOCS:

934 __call__ = _HeightIDW.__call__

935 height = _HeightIDW.height

936

937

938 class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny

939 if _FOR_DOCS:

940 __doc__ = HeightIDWflatLocal.__doc__

941 __init__ = HeightIDWflatLocal.__init__

942 __call__ = HeightIDWflatLocal.__call__

943 height = HeightIDWflatLocal.height

944

945

946 class HeightIDWkarney(_HeightIDW):

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

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

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

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

951 <https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}

952 Inverse method.

953

954 @see: L{HeightIDWcosineAndoyerLambert},

955 L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWdistanceTo},

956 L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas},

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

958 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

960 shepard_interp_2d/shepard_interp_2d.html>}.

961 '''

962 _datum = _WGS84

963 _Inverse1 = None

964 _ks = ()

965 _wrap = False

966

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

968 '''New L{HeightIDWkarney} interpolator.

969

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

971 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}

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

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

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

975 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).

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

977

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

979 an invalid B{C{knot}}, B{C{datum}} or

980 B{C{beta}}.

981

982 @raise ImportError: Package U{geographiclib

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

984

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

986 '''

987 n, self._ks = len2(knots)

988 if n < self._kmin:

989 raise _insufficientError(self._kmin, knots=n)

990 self._datum_setter(datum, self._ks)

991 self._Inverse1 = self.datum.ellipsoid.geodesic.Inverse1

992

993 self.beta = beta

994 if wrap:

995 self._wrap = True

996 if name:

997 self.name = name

998

999 def _distances(self, x, y): # (x, y) degrees

1000 for k in self._ks:

1001 # non-negative I{angular} distance in C{degrees}

1002 yield self._Inverse1(y, x, k.lat, k.lon, wrap=self._wrap)

1003

1004 @property # NOT _RO, no caching

1005 def _hs(self): # see HeightIDWdistanceTo

1006 for k in self._ks:

1007 yield k.height

1008

1009 def __call__(self, *llis):

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

1011

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

1013 C{LatLon}s).

1014

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

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

1017

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

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

1020 issue.

1021 '''

1022 def _xy2(lls):

1023 try: # like _xyhs above, but keeping degrees

1024 for i, ll in enumerate(lls):

1025 yield ll.lon, ll.lat

1026 except AttributeError as x:

1027 i = Fmt.SQUARE(llis=i)

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

1029

1030 _as, llis = _allis2(llis)

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

1032

1033 if _FOR_DOCS:

1034 height = _HeightIDW.height

1035

1036

1037 class HeightIDWthomas(_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{radians} from function L{pygeodesy.thomas_}.

1041

1042 @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert},

1043 L{HeightIDWdistanceTo}, L{HeightIDWflatLocal}, L{HeightIDWhubeny},

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

1045 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

1047 shepard_interp_2d/shepard_interp_2d.html>}.

1048 '''

1049 _datum = _WGS84

1050 _wrap = False

1051

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

1053 '''New L{HeightIDWthomas} interpolator.

1054

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

1056 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}

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

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

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

1060 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

1062

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

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

1065

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

1067 '''

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

1069 self._datum_setter(datum, knots)

1070

1071 def _distances(self, x, y): # (x, y) radians

1072 return self._distances_angular_datum_(thomas_, x, y)

1073

1074 if _FOR_DOCS:

1075 __call__ = _HeightIDW.__call__

1076 height = _HeightIDW.height

1077

1078

1079 class HeightIDWvincentys(_HeightIDW):

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

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

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

1083

1084 @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular},

1085 L{HeightIDWeuclidean}, L{HeightIDWflatPolar}, L{HeightIDWhaversine},

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

1087 geostatistics/Inverse-Distance-Weighting/index.html>} and

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

1089 shepard_interp_2d/shepard_interp_2d.html>}.

1090

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

1092 '''

1093 _wrap = False

1094

1095 def __init__(self, knots, beta=2, wrap=False, name=NN):

1096 '''New L{HeightIDWvincentys} interpolator.

1097

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

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

1100 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).

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

1102

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

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

1105 '''

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

1107

1108 def _distances(self, x, y): # (x, y) radians

1109 return self._distances_angular_(vincentys_, x, y)

1110

1111 if _FOR_DOCS:

1112 __call__ = _HeightIDW.__call__

1113 height = _HeightIDW.height

1114

1115

1116 class HeightLSQBiSpline(_HeightBase):

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

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

1119 interpolate.LSQSphereBivariateSpline.html>}.

1120 '''

1121 _kmin = 16 # k = 3, always

1122

1123 def __init__(self, knots, weight=None, name=NN):

1124 '''New L{HeightLSQBiSpline} interpolator.

1125

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

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

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

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

1130

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

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

1133

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

1135

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

1137 or not installed.

1138

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

1140

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

1142 warning as exception.

1143 '''

1144 np = self.numpy

1145 spi = self.scipy_interpolate

1146

1147 xs, ys, hs = self._xyhs3(knots)

1148 n = len(hs)

1149

1150 w = weight

1151 if isscalar(w):

1152 w = float(w)

1153 if w <= 0:

1154 raise HeightError(weight=w)

1155 w = [w] * n

1156 elif w is not None:

1157 m, w = len2(w)

1158 if m != n:

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

1160 w = map2(float, w)

1161 m = min(w)

1162 if m <= 0: # PYCHOK no cover

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

1164 raise HeightError(w, m)

1165 try:

1166 T = 1.0e-4 # like SciPy example

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

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

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

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

1171 except Exception as x:

1172 raise _SciPyIssue(x)

1173

1174 if name:

1175 self.name = name

1176

1177 def __call__(self, *llis):

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

1179

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

1181 C{LatLon}s).

1182

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

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

1185

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

1187 an invalid B{C{lli}}.

1188

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

1190

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

1192 warning as exception.

1193 '''

1194 return _HeightBase._eval(self, llis)

1195

1196 def height(self, lats, lons):

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

1198

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

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

1201

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

1203 interpolated heights (C{float}s).

1204

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

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

1207

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

1209

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

1211 warning as exception.

1212 '''

1213 return _HeightBase._height(self, lats, lons)

1214

1215

1216 class HeightSmoothBiSpline(_HeightBase):

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

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

1219 interpolate.SmoothSphereBivariateSpline.html>}.

1220 '''

1221 _kmin = 16 # k = 3, always

1222

1223 def __init__(self, knots, s=4, name=NN):

1224 '''New L{HeightSmoothBiSpline} interpolator.

1225

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

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

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

1229

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

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

1232

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

1234 or not installed.

1235

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

1237

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

1239 warning as exception.

1240 '''

1241 spi = self.scipy_interpolate

1242

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

1244

1245 xs, ys, hs = self._xyhs3(knots)

1246 try:

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

1248 eps=EPS, s=s).ev

1249 except Exception as x:

1250 raise _SciPyIssue(x)

1251

1252 if name:

1253 self.name = name

1254

1255 def __call__(self, *llis):

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

1257

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

1259 C{LatLon}s).

1260

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

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

1263

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

1265 an invalid B{C{lli}}.

1266

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

1268

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

1270 warning as exception.

1271 '''

1272 return _HeightBase._eval(self, llis)

1273

1274 def height(self, lats, lons):

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

1276

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

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

1279

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

1281 interpolated heights (C{float}s).

1282

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

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

1285

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

1287

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

1289 warning as exception.

1290 '''

1291 return _HeightBase._height(self, lats, lons)

1292

1293

1294__all__ += _ALL_DOCS(_HeightBase)

1295

1296# **) MIT License

1297#

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

1299#

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

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

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

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

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

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

1306#

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

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

1309#

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

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

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

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

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

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

1316# OTHER DEALINGS IN THE SOFTWARE.