Coverage for pygeodesy/heights.py: 96%
361 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
2# -*- coding: utf-8 -*-
4u'''Height interpolations of C{LatLon} points.
6Classes L{HeightCubic}, L{HeightIDWcosineAndoyerLambert},
7L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWcosineLaw},
8L{HeightIDWdistanceTo}, L{HeightIDWequirectangular}, L{HeightIDWeuclidean},
9L{HeightIDWflatLocal}, L{HeightIDWflatPolar}, L{HeightIDWhaversine},
10L{HeightIDWhubeny}, L{HeightIDWkarney}, L{HeightIDWthomas}, L{HeightIDWvincentys},
11L{HeightLinear}, L{HeightLSQBiSpline} and L{HeightSmoothBiSpline}
12to interpolate the height of C{LatLon} locations or separate
13lat-/longitudes from a set of C{LatLon} points with I{known heights}.
15Typical usage
16=============
181. Get or create a set of C{LatLon} points with I{known heights},
19called C{knots}. The C{knots} do not need to be ordered in any
20particular way.
222. Select one of the C{Height} classes for height interpolation
24C{>>> from pygeodesy import HeightCubic # or other Height... as HeightXyz}
263. Instantiate a height interpolator with the C{knots} and use keyword
27arguments to select different interpolation options
29C{>>> hinterpolator = HeightXyz(knots, **options)}
314. Get the interpolated height of other C{LatLon} location(s) with
33C{>>> ll = LatLon(1, 2, ...)}
34C{>>> h = hinterpolator(ll)}
36or
38C{>>> h0, h1, h2, ... = hinterpolator(ll0, ll1, ll2, ...)}
40or a list, tuple, generator, etc. of C{LatLon}s
42C{>>> hs = hinterpolator(lls)}
445. For separate lat- and longitudes invoke the C{.height} method
46C{>>> h = hinterpolator.height(lat, lon)}
48or as 2 lists, 2 tuples, etc.
50C{>>> hs = hinterpolator.height(lats, lons)}
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.
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}.
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
67from __future__ import division as _; del _ # PYCHOK semicolon
69from pygeodesy.basics import isscalar, len2, map1, map2, _xnumpy, _xscipy
70from pygeodesy.constants import EPS, PI, PI2, PI_2, _0_0, _90_0, _180_0
71from pygeodesy.datums import _ellipsoidal_datum, _WGS84
72from pygeodesy.errors import _AssertionError, LenError, PointsError, _SciPyIssue
73from pygeodesy.fmath import fidw, hypot2
74from pygeodesy.formy import cosineAndoyerLambert_, cosineForsytheAndoyerLambert_, \
75 cosineLaw_, euclidean_, flatPolar_, haversine_, \
76 _scale_rad, thomas_, vincentys_
77from pygeodesy.interns import NN, _COMMASPACE_, _cubic_, _datum_, _distanceTo_, \
78 _knots_, _linear_, _NOTEQUAL_, _scipy_, _SPACE_, _STAR_
79from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS
80from pygeodesy.named import _Named, notOverloaded
81from pygeodesy.points import LatLon_, Property_RO
82# from pygeodesy.props import Property_RO # from .points
83from pygeodesy.streprs import _boolkwds, Fmt
84from pygeodesy.units import Float_, Int_
85from pygeodesy.utily import radiansPI, radiansPI2, unrollPI
87__all__ = _ALL_LAZY.heights
88__version__ = '23.04.11'
90_error_ = 'error'
91_insufficient_ = 'insufficient'
92_llis_ = 'llis'
93_smoothing_ = 'smoothing'
96class HeightError(PointsError):
97 '''Height interpolator C{Height...} or interpolation issue.
98 '''
99 pass
102def _alist(ais):
103 # return list of floats, not numpy.float64s
104 return list(map(float, ais))
107def _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))
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
125 if n < m:
126 raise _insufficientError(m, Error=Error, llis=n)
127 return _as, llis
130def _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>
140def _atuple(ais):
141 # return tuple of floats, not numpy.float64s
142 return tuple(map(float, ais))
145def _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
153def _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)
159def _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
170def _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)
187def _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)
196class _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
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)
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
217 def _axyllis4(self, llis):
218 return _axyllis4(self.numpy.array, llis)
220 @Property_RO
221 def datum(self):
222 '''Get the datum (L{Datum} or C{None} if not applicable).
223 '''
224 return self._datum
226 def _ev(self, *args): # PYCHOK no cover
227 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
228 '''
229 notOverloaded(self, *args)
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)
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)
252 @Property_RO
253 def kmin(self):
254 '''Get the minimum number of knots (C{int}).
255 '''
256 return self._kmin
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_)
271 sp = _xscipy(self.__class__, 1, 2)
272 np = _xnumpy(self.__class__, 1, 9)
274 _HeightBase._np_sp = t = np, sp
275 return t
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
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
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
299 def _xyhs3(self, knots):
300 return _xyhs3(self.numpy.array, self._kmin, knots)
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
309class 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
318 def __init__(self, knots, name=NN):
319 '''New L{HeightCubic} interpolator.
321 @arg knots: The points with known height (C{LatLon}s).
322 @kwarg name: Optional name for this height interpolator (C{str}).
324 @raise HeightError: Insufficient number of B{C{knots}} or
325 invalid B{C{knot}}.
327 @raise ImportError: Package C{numpy} or C{scipy} not found
328 or not installed.
330 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
332 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
333 as exception.
334 '''
335 spi = self.scipy_interpolate
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)
343 if name:
344 self.name = name
346 def __call__(self, *llis):
347 '''Interpolate the height for one or several locations.
349 @arg llis: The location or locations (C{LatLon}, ... or
350 C{LatLon}s).
352 @return: A single interpolated height (C{float}) or a list
353 or tuple of interpolated heights (C{float}s).
355 @raise HeightError: Insufficient number of B{C{llis}} or
356 an invalid B{C{lli}}.
358 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
360 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
361 as exception.
362 '''
363 return _HeightBase._eval(self, llis)
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)
370 def height(self, lats, lons):
371 '''Interpolate the height for one or several lat-/longitudes.
373 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
374 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
376 @return: A single interpolated height (C{float}) or a list of
377 interpolated heights (C{float}s).
379 @raise HeightError: Insufficient or non-matching number of
380 B{C{lats}} and B{C{lons}}.
382 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
384 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
385 as exception.
386 '''
387 return _HeightBase._height(self, lats, lons)
390class 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
398 def __init__(self, knots, name=NN):
399 '''New L{HeightLinear} interpolator.
401 @arg knots: The points with known height (C{LatLon}s).
402 @kwarg name: Optional name for this height interpolator (C{str}).
404 @raise HeightError: Insufficient number of B{C{knots}} or
405 an invalid B{C{knot}}.
407 @raise ImportError: Package C{numpy} or C{scipy} not found
408 or not installed.
410 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
412 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
413 as exception.
414 '''
415 HeightCubic.__init__(self, knots, name=name)
417 if _FOR_DOCS:
418 __call__ = HeightCubic.__call__
419 height = HeightCubic.height
422class _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
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)
442 def __call__(self, *llis):
443 '''Interpolate the height for one or several locations.
445 @arg llis: The location or locations (C{LatLon}, ... or
446 C{LatLon}s).
448 @return: A single interpolated height (C{float}) or a list
449 or tuple of interpolated heights (C{float}s).
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))
458 def _datum_setter(self, datum, knots):
459 '''(INTERNAL) Set the datum.
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)
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)
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)
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)
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)
492 @property
493 def beta(self):
494 '''Get the inverse distance power (C{int}).
495 '''
496 return self._beta
498 @beta.setter # PYCHOK setter!
499 def beta(self, beta):
500 '''Set the inverse distance power (C{int} 1, 2, or 3).
502 @raise HeightError: Invalid B{C{beta}}.
503 '''
504 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
506 def height(self, lats, lons):
507 '''Interpolate the height for one or several lat-/longitudes.
509 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
510 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
512 @return: A single interpolated height (C{float}) or a list of
513 interpolated heights (C{float}s).
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)
522class 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_}.
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
538 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
539 '''New L{HeightIDWcosineAndoyerLambert} interpolator.
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}).
549 @raise HeightError: Insufficient number of B{C{knots}} or
550 an invalid B{C{knot}} or B{C{beta}}.
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)
557 def _distances(self, x, y): # (x, y) radians
558 return self._distances_angular_datum_(cosineAndoyerLambert_, x, y)
560 if _FOR_DOCS:
561 __call__ = _HeightIDW.__call__
562 height = _HeightIDW.height
565class 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_}.
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
581 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
582 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
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}).
592 @raise HeightError: Insufficient number of B{C{knots}} or
593 an invalid B{C{knot}} or B{C{beta}}.
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)
600 def _distances(self, x, y): # (x, y) radians
601 return self._distances_angular_datum_(cosineForsytheAndoyerLambert_, x, y)
603 if _FOR_DOCS:
604 __call__ = _HeightIDW.__call__
605 height = _HeightIDW.height
608class 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_}.
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>}.
620 @note: See note at function L{pygeodesy.vincentys_}.
621 '''
622 _wrap = False
624 def __init__(self, knots, beta=2, wrap=False, name=NN):
625 '''New L{HeightIDWcosineLaw} interpolator.
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}).
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)
637 def _distances(self, x, y): # (x, y) radians
638 return self._distances_angular_(cosineLaw_, x, y)
640 if _FOR_DOCS:
641 __call__ = _HeightIDW.__call__
642 height = _HeightIDW.height
645class 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}.
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 = ()
661 def __init__(self, knots, beta=2, name=NN, **distanceTo_kwds):
662 '''New L{HeightIDWdistanceTo} interpolator.
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.
671 @raise HeightError: Insufficient number of B{C{knots}} or
672 an invalid B{C{knot}} or B{C{beta}}.
674 @raise ImportError: Package U{geographiclib
675 <https://PyPI.org/project/geographiclib>} missing
676 iff B{C{points}} are L{ellipsoidalKarney.LatLon}s.
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_)
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
696 self.beta = beta
697 if name:
698 self.name = name
699 if distanceTo_kwds:
700 self._distanceTo_kwds = distanceTo_kwds
702 def __call__(self, *llis):
703 '''Interpolate the height for one or several locations.
705 @arg llis: The location or locations (C{LatLon}, ... or
706 C{LatLon}s).
708 @return: A single interpolated height (C{float}) or a list
709 or tuple of interpolated heights (C{float}s).
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))
718 if _FOR_DOCS:
719 height = _HeightIDW.height
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))
730 @property # NOT _RO, no caching
731 def _hs(self): # see HeightIDWkarney
732 for k in self._ks:
733 yield k.height
736class 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_}.
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
752 def __init__(self, knots, adjust=True, wrap=False, name=NN):
753 '''New L{HeightIDWequirectangular} interpolator.
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}).
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)
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
774 if _FOR_DOCS:
775 __call__ = _HeightIDW.__call__
776 height = _HeightIDW.height
779class 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_}.
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
793 def __init__(self, knots, adjust=True, beta=2, name=NN):
794 '''New L{HeightIDWeuclidean} interpolator.
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}).
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)
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)
811 if _FOR_DOCS:
812 __call__ = _HeightIDW.__call__
813 height = _HeightIDW.height
816class 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_}.
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
832 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
833 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
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}).
843 @raise HeightError: Insufficient number of B{C{knots}} or
844 an invalid B{C{knot}} or B{C{beta}}.
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)
851 def _distances(self, x, y): # (x, y) radians
852 return self._distances_angular_(self._hubeny_2, x, y) # radians**2
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
860 if _FOR_DOCS:
861 __call__ = _HeightIDW.__call__
862 height = _HeightIDW.height
865class 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_}.
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
880 def __init__(self, knots, beta=2, wrap=False, name=NN):
881 '''New L{HeightIDWflatPolar} interpolator.
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}).
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)
893 def _distances(self, x, y): # (x, y) radians
894 return self._distances_angular_(flatPolar_, x, y)
896 if _FOR_DOCS:
897 __call__ = _HeightIDW.__call__
898 height = _HeightIDW.height
901class 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_}.
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>}.
913 @note: See note at function L{pygeodesy.vincentys_}.
914 '''
915 _wrap = False
917 def __init__(self, knots, beta=2, wrap=False, name=NN):
918 '''New L{HeightIDWhaversine} interpolator.
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}).
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)
930 def _distances(self, x, y): # (x, y) radians
931 return self._distances_angular_(haversine_, x, y)
933 if _FOR_DOCS:
934 __call__ = _HeightIDW.__call__
935 height = _HeightIDW.height
938class 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
946class 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/Python/doc/code.html>}
952 Inverse method.
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
967 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
968 '''New L{HeightIDWkarney} interpolator.
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}).
978 @raise HeightError: Insufficient number of B{C{knots}} or
979 an invalid B{C{knot}}, B{C{datum}} or
980 B{C{beta}}.
982 @raise ImportError: Package U{geographiclib
983 <https://PyPI.org/project/geographiclib>} missing.
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
993 self.beta = beta
994 if wrap:
995 self._wrap = True
996 if name:
997 self.name = name
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)
1004 @property # NOT _RO, no caching
1005 def _hs(self): # see HeightIDWdistanceTo
1006 for k in self._ks:
1007 yield k.height
1009 def __call__(self, *llis):
1010 '''Interpolate the height for one or several locations.
1012 @arg llis: The location or locations (C{LatLon}, ... or
1013 C{LatLon}s).
1015 @return: A single interpolated height (C{float}) or a list
1016 or tuple of interpolated heights (C{float}s).
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)
1030 _as, llis = _allis2(llis)
1031 return _as(map(self._hIDW, *zip(*_xy2(llis))))
1033 if _FOR_DOCS:
1034 height = _HeightIDW.height
1037class 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_}.
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
1052 def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
1053 '''New L{HeightIDWthomas} interpolator.
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}).
1063 @raise HeightError: Insufficient number of B{C{knots}} or
1064 an invalid B{C{knot}} or B{C{beta}}.
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)
1071 def _distances(self, x, y): # (x, y) radians
1072 return self._distances_angular_datum_(thomas_, x, y)
1074 if _FOR_DOCS:
1075 __call__ = _HeightIDW.__call__
1076 height = _HeightIDW.height
1079class 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_}.
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>}.
1091 @note: See note at function L{pygeodesy.vincentys_}.
1092 '''
1093 _wrap = False
1095 def __init__(self, knots, beta=2, wrap=False, name=NN):
1096 '''New L{HeightIDWvincentys} interpolator.
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}).
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)
1108 def _distances(self, x, y): # (x, y) radians
1109 return self._distances_angular_(vincentys_, x, y)
1111 if _FOR_DOCS:
1112 __call__ = _HeightIDW.__call__
1113 height = _HeightIDW.height
1116class 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
1123 def __init__(self, knots, weight=None, name=NN):
1124 '''New L{HeightLSQBiSpline} interpolator.
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}).
1131 @raise HeightError: Insufficient number of B{C{knots}} or
1132 an invalid B{C{knot}} or B{C{weight}}.
1134 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
1136 @raise ImportError: Package C{numpy} or C{scipy} not found
1137 or not installed.
1139 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
1141 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
1142 warning as exception.
1143 '''
1144 np = self.numpy
1145 spi = self.scipy_interpolate
1147 xs, ys, hs = self._xyhs3(knots)
1148 n = len(hs)
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)
1174 if name:
1175 self.name = name
1177 def __call__(self, *llis):
1178 '''Interpolate the height for one or several locations.
1180 @arg llis: The location or locations (C{LatLon}, ... or
1181 C{LatLon}s).
1183 @return: A single interpolated height (C{float}) or a list
1184 or tuple of interpolated heights (C{float}s).
1186 @raise HeightError: Insufficient number of B{C{llis}} or
1187 an invalid B{C{lli}}.
1189 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
1191 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
1192 warning as exception.
1193 '''
1194 return _HeightBase._eval(self, llis)
1196 def height(self, lats, lons):
1197 '''Interpolate the height for one or several lat-/longitudes.
1199 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
1200 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
1202 @return: A single interpolated height (C{float}) or a list of
1203 interpolated heights (C{float}s).
1205 @raise HeightError: Insufficient or non-matching number of
1206 B{C{lats}} and B{C{lons}}.
1208 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
1210 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
1211 warning as exception.
1212 '''
1213 return _HeightBase._height(self, lats, lons)
1216class 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
1223 def __init__(self, knots, s=4, name=NN):
1224 '''New L{HeightSmoothBiSpline} interpolator.
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}).
1230 @raise HeightError: Insufficient number of B{C{knots}} or
1231 an invalid B{C{knot}} or B{C{s}}.
1233 @raise ImportError: Package C{numpy} or C{scipy} not found
1234 or not installed.
1236 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
1238 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
1239 warning as exception.
1240 '''
1241 spi = self.scipy_interpolate
1243 s = Float_(s, name=_smoothing_, Error=HeightError, low=4)
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)
1252 if name:
1253 self.name = name
1255 def __call__(self, *llis):
1256 '''Interpolate the height for one or several locations.
1258 @arg llis: The location or locations (C{LatLon}, ... or
1259 C{LatLon}s).
1261 @return: A single interpolated height (C{float}) or a list
1262 or tuple of interpolated heights (C{float}s).
1264 @raise HeightError: Insufficient number of B{C{llis}} or
1265 an invalid B{C{lli}}.
1267 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
1269 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
1270 warning as exception.
1271 '''
1272 return _HeightBase._eval(self, llis)
1274 def height(self, lats, lons):
1275 '''Interpolate the height for one or several lat-/longitudes.
1277 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
1278 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
1280 @return: A single interpolated height (C{float}) or a list of
1281 interpolated heights (C{float}s).
1283 @raise HeightError: Insufficient or non-matching number of
1284 B{C{lats}} and B{C{lons}}.
1286 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
1288 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
1289 warning as exception.
1290 '''
1291 return _HeightBase._height(self, lats, lons)
1294__all__ += _ALL_DOCS(_HeightBase)
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.