Coverage for pygeodesy/heights.py: 95%
318 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-11-12 16:17 -0500
« prev ^ index » next coverage.py v7.6.1, created at 2024-11-12 16:17 -0500
2# -*- coding: utf-8 -*-
4u'''Height interpolations at C{LatLon} points from known C{knots}.
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}, called
19C{knots}. The C{knots} do not need to be ordered in any particular way.
21C{>>> ...}
232. Select one of the C{Height} classes for height interpolation
25C{>>> from pygeodesy import HeightCubic # or other Height... as HeightXyz}
273. Instantiate a height interpolator with the C{knots} and use keyword
28arguments to select different interpolation options
30C{>>> hinterpolator = HeightXyz(knots, **options)}
324. Get the interpolated height of other C{LatLon} location(s) with
34C{>>> ll = LatLon(1, 2, ...)}
36C{>>> h = hinterpolator(ll)}
38or
40C{>>> h0, h1, h2, ... = hinterpolator(ll0, ll1, ll2, ...)}
42or a list, tuple, generator, etc. of C{LatLon}s
44C{>>> hs = hinterpolator(lls)}
465. For separate lat- and longitudes invoke the C{.height} method
48C{>>> h = hinterpolator.height(lat, lon)}
50or as 2 lists, 2 tuples, etc.
52C{>>> hs = hinterpolator.height(lats, lons)}
54@note: Classes L{HeightCubic} and L{HeightLinear} require package U{numpy
55 <https://PyPI.org/project/numpy>}, classes L{HeightLSQBiSpline} and
56 L{HeightSmoothBiSpline} require package U{scipy<https://SciPy.org>}.
57 Classes L{HeightIDWkarney} and L{HeightIDWdistanceTo} -if used with
58 L{ellipsoidalKarney.LatLon} points- require I{Karney}'s U{geographiclib
59 <https://PyPI.org/project/geographiclib>} to be installed.
61@note: Errors from C{scipy} are raised as L{SciPyError}s. Warnings issued
62 by C{scipy} can be thrown as L{SciPyWarning} exceptions, provided
63 Python C{warnings} are filtered accordingly, see L{SciPyWarning}.
65@see: U{SciPy<https://docs.SciPy.org/doc/scipy/reference/interpolate.html>}
66 Interpolation.
67'''
68# make sure int/int division yields float quotient, see .basics
69from __future__ import division as _; del _ # PYCHOK semicolon
71from pygeodesy.basics import isscalar, len2, map1, map2, _xnumpy, _xscipy
72from pygeodesy.constants import EPS, PI, PI_2, PI2, _0_0, _90_0, _180_0
73from pygeodesy.datums import _ellipsoidal_datum, _WGS84
74from pygeodesy.errors import _AssertionError, LenError, PointsError, \
75 _SciPyIssue, _xattr, _xkwds, _xkwds_get, _xkwds_item2
76# from pygeodesy.fmath import fidw # _MODS
77# from pygeodesy.formy import cosineAndoyerLambert, cosineForsytheAndoyerLambert, \
78# cosineLaw, equirectangular4, euclidean, flatLocal, \
79# flatPolar, haversine, thomas, vincentys # _MODS.into
80# from pygeodesy.internals import _version2 # _MODS
81from pygeodesy.interns import NN, _COMMASPACE_, _insufficient_, _NOTEQUAL_, \
82 _PLUS_, _scipy_, _SPACE_, _STAR_
83from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, _FOR_DOCS
84from pygeodesy.named import _name2__, _Named
85from pygeodesy.points import _distanceTo, LatLon_, Fmt, radians, _Wrap
86from pygeodesy.props import Property_RO, property_RO, property_ROver
87# from pygeodesy.streprs import Fmt # from .points
88from pygeodesy.units import _isDegrees, Float_, Int_
89# from pygeodesy.utily import _Wrap # from .points
91# from math import radians # from .points
93__all__ = _ALL_LAZY.heights
94__version__ = '24.11.02'
96_error_ = 'error'
97_formy = _MODS.into(formy=__name__)
98_linear_ = 'linear'
99_llis_ = 'llis'
102class HeightError(PointsError):
103 '''Height interpolator C{Height...} or interpolation issue.
104 '''
105 pass
108def _alist(ais):
109 # return list of floats, not numpy.float64s
110 return list(map(float, ais))
113def _ascalar(ais): # in .geoids
114 # return single float, not numpy.float64
115 ais = list(ais) # np.array, etc. to list
116 if len(ais) != 1:
117 n = Fmt.PAREN(len=repr(ais))
118 t = _SPACE_(len(ais), _NOTEQUAL_, 1)
119 raise _AssertionError(n, txt=t)
120 return float(ais[0]) # remove np.<type>
123def _atuple(ais):
124 # return tuple of floats, not numpy.float64s
125 return tuple(map(float, ais))
128def _as_llis2(llis, m=1, Error=HeightError): # in .geoids
129 # determine return type and convert lli C{LatLon}s to list
130 if not isinstance(llis, tuple): # llis are *args
131 n = Fmt.PAREN(type_=_STAR_(NN, _llis_))
132 raise _AssertionError(n, txt=repr(llis))
134 n = len(llis)
135 if n == 1: # convert single lli to 1-item list
136 llis = llis[0]
137 try:
138 n, llis = len2(llis)
139 _as = _alist # return list of interpolated heights
140 except TypeError: # single lli
141 n, llis = 1, [llis]
142 _as = _ascalar # return single interpolated heights
143 else: # of 0, 2 or more llis
144 _as = _atuple # return tuple of interpolated heights
146 if n < m:
147 raise _InsufficientError(m, Error=Error, llis=n)
148 return _as, llis
151def _InsufficientError(need, Error=HeightError, **name_value): # PYCHOK no cover
152 # create an insufficient Error instance
153 t = _COMMASPACE_(_insufficient_, str(need) + _PLUS_)
154 return Error(txt=t, **name_value)
157def _orderedup(ts, lo=EPS, hi=PI2-EPS):
158 # clip, order and remove duplicates
159 return sorted(set(max(lo, min(hi, t)) for t in ts)) # list
162def _xyhs(wrap=False, _lat=_90_0, _lon=_180_0, **name_lls):
163 # map (lat, lon, h) to (x, y, h) in radians, offset
164 # x as 0 <= lon <= PI2 and y as 0 <= lat <= PI
165 name, lls = _xkwds_item2(name_lls)
166 _r, _w = radians, _Wrap._latlonop(wrap)
167 try:
168 for i, ll in enumerate(lls):
169 y, x = _w(ll.lat, ll.lon)
170 yield max(_0_0, _r(x + _lon)), \
171 max(_0_0, _r(y + _lat)), ll.height
172 except Exception as e:
173 i = Fmt.INDEX(name, i)
174 raise HeightError(i, ll, cause=e)
177class _HeightNamed(_Named): # in .geoids
178 '''(INTERNAL) Interpolator base class.
179 '''
180 _datum = _WGS84 # default
181 _Error = HeightError
182 _kmin = 2 # min number of knots
184 _LLis = LatLon_ # ._height class
185 _np_sp = None # (numpy, scipy)
186 _wrap = None # wrap knots and llis
188 def _as_lls(self, lats, lons): # in .geoids
189 LLis, d = self._LLis, self.datum
190 if _isDegrees(lats) and _isDegrees(lons):
191 llis = LLis(lats, lons, datum=d)
192 else:
193 n, lats = len2(lats)
194 m, lons = len2(lons)
195 if n != m: # format a LenError, but raise self._Error
196 e = LenError(self.__class__, lats=n, lons=m, txt=None)
197 raise self._Error(str(e))
198 llis = [LLis(*t, datum=d) for t in zip(lats, lons)]
199 return llis
201 @property_RO
202 def datum(self):
203 '''Get the C{datum} setting or the default (L{Datum}).
204 '''
205 return self._datum
207 @property_RO
208 def kmin(self):
209 '''Get the minimum number of knots (C{int}).
210 '''
211 return self._kmin
213 @property_RO
214 def wrap(self):
215 '''Get the C{wrap} setting (C{bool}) or C{None}.
216 '''
217 return self._wrap
220class _HeightBase(_HeightNamed): # in .geoids
221 '''(INTERNAL) Interpolator base class.
222 '''
223 _k2interp2d = {-1: _linear_, # in .geoids._GeoidBase.__init__
224 -2: _linear_, # for backward compatibility
225 -3: 'cubic',
226 -5: 'quintic'}
228 def __call__(self, *llis, **wrap): # PYCHOK no cover
229 '''Interpolate the height for one or several locations. I{Must be overloaded}.
231 @arg llis: One or more locations (C{LatLon}s), all positional.
232 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
233 locations (C{bool}), overriding the B{C{knots}}
234 setting.
236 @return: A single interpolated height (C{float}) or a list
237 or tuple of interpolated heights (C{float}s).
239 @raise HeightError: Insufficient number of B{C{llis}} or
240 an invalid B{C{lli}}.
242 @raise SciPyError: A C{scipy} issue.
244 @raise SciPyWarning: A C{scipy} warning as exception.
245 '''
246 self._notOverloaded(callername='__call__', *llis, **wrap)
248 def _as_xyllis4(self, llis, **wrap):
249 # convert lli C{LatLon}s to tuples or C{NumPy} arrays of
250 # C{SciPy} sphericals and determine the return type
251 atype = self.numpy.array
252 wrap = _xkwds(wrap, wrap=self._wrap)
253 _as, llis = _as_llis2(llis)
254 xis, yis, _ = zip(*_xyhs(llis=llis, **wrap)) # PYCHOK yield
255 return _as, atype(xis), atype(yis), llis
257 def _ev(self, *args): # PYCHOK no cover
258 '''(INTERNAL) I{Must be overloaded}.'''
259 self._notOverloaded(*args)
261 def _evalls(self, llis, **wrap): # XXX single arg, not *args
262 _as, xis, yis, _ = self._as_xyllis4(llis, **wrap)
263 try: # SciPy .ev signature: y first, then x!
264 return _as(self._ev(yis, xis))
265 except Exception as x:
266 raise _SciPyIssue(x, self._ev_name)
268 def _ev2d(self, x, y): # PYCHOK no cover
269 '''(INTERNAL) I{Must be overloaded}.'''
270 self._notOverloaded(x, y)
272 @property_RO
273 def _ev_name(self):
274 '''(INTERNAL) Get the name of the C{.ev} method.
275 '''
276 _ev = str(self._ev)
277 if _scipy_ not in _ev:
278 _ev = str(self._ev2d)
279 # '<scipy.interpolate._interpolate.interp2d object at ...>
280 # '<function _HeightBase._interp2d.<locals>._bisplev at ...>
281 # '<bound method BivariateSpline.ev of ... object at ...>
282 _ev = _ev[1:].split(None, 4)
283 return Fmt.PAREN(_ev['sfb'.index(_ev[0][0])])
285 def height(self, lats, lons, **wrap): # PYCHOK no cover
286 '''Interpolate the height for one or several lat-/longitudes. I{Must be overloaded}.
288 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
289 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
290 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
291 and B{C{lons}} locationts (C{bool}), overriding
292 the B{C{knots}} setting.
294 @return: A single interpolated height (C{float}) or a list of
295 interpolated heights (C{float}s).
297 @raise HeightError: Insufficient or non-matching number of
298 B{C{lats}} and B{C{lons}}.
300 @raise SciPyError: A C{scipy} issue.
302 @raise SciPyWarning: A C{scipy} warning as exception.
303 '''
304 lls = self._as_lls(lats, lons)
305 return self(lls, **wrap) # __call__(ll) or __call__(lls)
307 def _interp2d(self, xs, ys, hs, kind=-3):
308 '''Create a C{scipy.interpolate.interp2d} or C{-.bisplrep/-ev}
309 interpolator before, respectively since C{SciPy} version 1.14.
310 '''
311 try:
312 spi = self.scipy_interpolate
313 if self._scipy_version() < (1, 14) and kind in self._k2interp2d:
314 # SciPy.interpolate.interp2d kind 'linear', 'cubic' or 'quintic'
315 # DEPRECATED since scipy 1.10, removed altogether in 1.14
316 self._ev2d = spi.interp2d(xs, ys, hs, kind=self._k2interp2d[kind])
318 else: # <https://scipy.GitHub.io/devdocs/tutorial/interpolate/interp_transition_guide.html>
319 k = self._kxky(abs(kind))
320 # spi.RectBivariateSpline needs strictly ordered xs and ys
321 r = spi.bisplrep(xs, ys, hs.T, kx=k, ky=k)
323 def _bisplev(x, y):
324 return spi.bisplev(x, y, r) # .T
326 self._ev2d = _bisplev
328 except Exception as x:
329 raise _SciPyIssue(x, self._ev_name)
331 def _kxky(self, kind):
332 return Int_(kind=kind, low=1, high=5, Error=self._Error)
334 def _np_sp2(self, throwarnings=False): # PYCHOK no cover
335 '''(INTERNAL) Import C{numpy} and C{scipy}, once.
336 '''
337 # raise SciPyWarnings, but not if
338 # scipy has already been imported
339 if throwarnings: # PYCHOK no cover
340 import sys
341 if _scipy_ not in sys.modules:
342 import warnings
343 warnings.filterwarnings(_error_)
344 return self.numpy, self.scipy
346 @property_ROver
347 def numpy(self):
348 '''Get the C{numpy} module or C{None}.
349 '''
350 return _xnumpy(self.__class__, 1, 9) # overwrite property_ROver
352 @property_ROver
353 def scipy(self):
354 '''Get the C{scipy} module or C{None}.
355 '''
356 return _xscipy(self.__class__, 1, 2) # overwrite property_ROver
358 @property_ROver
359 def scipy_interpolate(self):
360 '''Get the C{scipy.interpolate} module or C{None}.
361 '''
362 _ = self.scipy
363 import scipy.interpolate as spi # scipy 1.2.2
364 return spi # overwrite property_ROver
366 def _scipy_version(self, **n):
367 '''Get the C{scipy} version as 2- or 3-tuple C{(major, minor, micro)}.
368 '''
369 return _MODS.internals._version2(self.scipy.version.version, **n)
371 def _xyhs3(self, knots, wrap=False, **name):
372 # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals
373 xs, ys, hs = zip(*_xyhs(knots=knots, wrap=wrap)) # PYCHOK yield
374 n = len(hs)
375 if n < self.kmin:
376 raise _InsufficientError(self.kmin, knots=n)
377 if name:
378 self.name = name
379 return map1(self.numpy.array, xs, ys, hs)
382class HeightCubic(_HeightBase):
383 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
384 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
385 C{kind='cubic'} or U{bisplrep/-ev<https://docs.SciPy.org/doc/scipy/
386 reference/generated/scipy.interpolate.interp2d.html>} C{kx=ky=3}.
387 '''
388 _kind = -3
389 _kmin = 16
391 def __init__(self, knots, **name_wrap):
392 '''New L{HeightCubic} interpolator.
394 @arg knots: The points with known height (C{LatLon}s).
395 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator (C{str})
396 and keyword argument C{b{wrap}=False} to wrap or I{normalize} all
397 B{C{knots}} and B{C{llis}} locations iff C{True} (C{bool}).
399 @raise HeightError: Insufficient number of B{C{knots}} or invalid B{C{knot}}.
401 @raise ImportError: Package C{numpy} or C{scipy} not found or not installed.
403 @raise SciPyError: A C{scipy} issue.
405 @raise SciPyWarning: A C{scipy} warning as exception.
406 '''
407 xs_yx_hs = self._xyhs3(knots, **name_wrap)
408 self._interp2d(*xs_yx_hs, kind=self._kind)
410 def __call__(self, *llis, **wrap):
411 '''Interpolate the height for one or several locations.
413 @see: L{Here<_HeightBase.__call__>} for further details.
414 '''
415 return self._evalls(llis, **wrap)
417 def _ev(self, yis, xis): # PYCHOK overwritten with .RectBivariateSpline.ev
418 # to make SciPy .interp2d single (x, y) signature
419 # match SciPy .ev signature(ys, xs), flipped multiples
420 return map(self._ev2d, xis, yis)
423class HeightLinear(HeightCubic):
424 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
425 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
426 C{kind='linear'} or U{bisplrep/-ev<https://docs.SciPy.org/doc/scipy/
427 reference/generated/scipy.interpolate.interp2d.html>} C{kx=ky=1}.
428 '''
429 _kind = -1
430 _kmin = 2
432 def __init__(self, knots, **name_wrap):
433 '''New L{HeightLinear} interpolator.
435 @see: L{Here<HeightCubic.__init__>} for all details.
436 '''
437 HeightCubic.__init__(self, knots, **name_wrap)
439 if _FOR_DOCS:
440 __call__ = HeightCubic.__call__
441 height = HeightCubic.height
444class HeightLSQBiSpline(_HeightBase):
445 '''Height interpolator using C{SciPy} U{LSQSphereBivariateSpline
446 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
447 interpolate.LSQSphereBivariateSpline.html>}.
448 '''
449 _kmin = 16 # k = 3, always
451 def __init__(self, knots, weight=None, low=1e-4, **name_wrap):
452 '''New L{HeightLSQBiSpline} interpolator.
454 @arg knots: The points with known height (C{LatLon}s).
455 @kwarg weight: Optional weight or weights for each B{C{knot}}
456 (C{scalar} or C{scalar}s).
457 @kwarg low: Optional lower bound for I{ordered knots} (C{radians}).
458 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
459 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
460 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
461 C{True} (C{bool}).
463 @raise HeightError: Insufficient number of B{C{knots}} or an invalid
464 B{C{knot}}, B{C{weight}} or B{C{eps}}.
466 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
468 @raise ImportError: Package C{numpy} or C{scipy} not found or not
469 installed.
471 @raise SciPyError: A C{scipy} issue.
473 @raise SciPyWarning: A C{scipy} warning as exception.
474 '''
475 np = self.numpy
476 spi = self.scipy_interpolate
478 xs, ys, hs = self._xyhs3(knots, **name_wrap)
479 n = len(hs)
481 w = weight
482 if isscalar(w):
483 w = float(w)
484 if w <= 0:
485 raise HeightError(weight=w)
486 w = (w,) * n
487 elif w is not None:
488 m, w = len2(w)
489 if m != n:
490 raise LenError(HeightLSQBiSpline, weight=m, knots=n)
491 w = map2(float, w)
492 m = min(w)
493 if m <= 0: # PYCHOK no cover
494 w = Fmt.INDEX(weight=w.index(m))
495 raise HeightError(w, m)
496 try:
497 if not EPS < low < (PI_2 - EPS): # 1e-4 like SciPy example
498 raise HeightError(low=low)
499 ps = np.array(_orderedup(xs, low, PI2 - low))
500 ts = np.array(_orderedup(ys, low, PI - low))
501 self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs,
502 ts, ps, eps=EPS, w=w).ev
503 except Exception as x:
504 raise _SciPyIssue(x, self._ev_name)
506 def __call__(self, *llis, **wrap):
507 '''Interpolate the height for one or several locations.
509 @see: L{Here<_HeightBase.__call__>} for further details.
510 '''
511 return self._evalls(llis, **wrap)
514class HeightSmoothBiSpline(_HeightBase):
515 '''Height interpolator using C{SciPy} U{SmoothSphereBivariateSpline
516 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
517 interpolate.SmoothSphereBivariateSpline.html>}.
518 '''
519 _kmin = 16 # k = 3, always
521 def __init__(self, knots, s=4, **name_wrap):
522 '''New L{HeightSmoothBiSpline} interpolator.
524 @arg knots: The points with known height (C{LatLon}s).
525 @kwarg s: The spline smoothing factor (C{scalar}), default C{4}.
526 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
527 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
528 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
529 C{True} (C{bool}).
531 @raise HeightError: Insufficient number of B{C{knots}} or an invalid
532 B{C{knot}} or B{C{s}}.
534 @raise ImportError: Package C{numpy} or C{scipy} not found or not
535 installed.
537 @raise SciPyError: A C{scipy} issue.
539 @raise SciPyWarning: A C{scipy} warning as exception.
540 '''
541 spi = self.scipy_interpolate
543 s = Float_(smoothing=s, Error=HeightError, low=4)
545 xs, ys, hs = self._xyhs3(knots, **name_wrap)
546 try:
547 self._ev = spi.SmoothSphereBivariateSpline(ys, xs, hs,
548 eps=EPS, s=s).ev
549 except Exception as x:
550 raise _SciPyIssue(x, self._ev_name)
552 def __call__(self, *llis, **wrap):
553 '''Interpolate the height for one or several locations.
555 @see: L{Here<_HeightBase.__call__>} for further details.
556 '''
557 return self._evalls(llis, **wrap)
560class _HeightIDW(_HeightNamed):
561 '''(INTERNAL) Base class for U{Inverse Distance Weighting
562 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) height
563 interpolators.
565 @see: U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
566 geostatistics/Inverse-Distance-Weighting/index.html>},
567 U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
568 shepard_interp_2d/shepard_interp_2d.html>} and other C{_HeightIDW*}
569 classes.
570 '''
571 _beta = 0 # fidw inverse power
572 _func = None # formy function
573 _knots = () # knots list or tuple
574 _kwds = {} # func_ options
576 def __init__(self, knots, beta=2, **name__kwds):
577 '''New C{_HeightIDW*} interpolator.
579 @arg knots: The points with known height (C{LatLon}s).
580 @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
581 @kwarg name__kwds: Optional C{B{name}=NN} for this height interpolator
582 (C{str}) and any keyword arguments for the distance function,
583 retrievable with property C{kwds}.
585 @raise HeightError: Insufficient number of B{C{knots}} or an invalid
586 B{C{knot}} or B{C{beta}}.
587 '''
588 name, kwds = _name2__(**name__kwds)
589 if name:
590 self.name = name
592 n, self._knots = len2(knots)
593 if n < self.kmin:
594 raise _InsufficientError(self.kmin, knots=n)
595 self.beta = beta
596 self._kwds = kwds or {}
598 def __call__(self, *llis, **wrap):
599 '''Interpolate the height for one or several locations.
601 @arg llis: One or more locations (C{LatLon}s), all positional.
602 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
603 locations (C{bool}).
605 @return: A single interpolated height (C{float}) or a list
606 or tuple of interpolated heights (C{float}s).
608 @raise HeightError: Insufficient number of B{C{llis}}, an
609 invalid B{C{lli}} or L{pygeodesy.fidw}
610 issue.
611 '''
612 def _xy2(wrap=False):
613 _w = _Wrap._latlonop(wrap)
614 try: # like _xyhs above, but degrees
615 for i, ll in enumerate(llis):
616 yield _w(ll.lon, ll.lat)
617 except Exception as x:
618 i = Fmt.INDEX(llis=i)
619 raise HeightError(i, ll, cause=x)
621 _as, llis = _as_llis2(llis)
622 return _as(map(self._hIDW, *zip(*_xy2(**wrap))))
624 @property_RO
625 def adjust(self):
626 '''Get the C{adjust} setting (C{bool}) or C{None}.
627 '''
628 return _xkwds_get(self._kwds, adjust=None)
630 @property
631 def beta(self):
632 '''Get the inverse distance power (C{int}).
633 '''
634 return self._beta
636 @beta.setter # PYCHOK setter!
637 def beta(self, beta):
638 '''Set the inverse distance power (C{int} 1, 2, or 3).
640 @raise HeightError: Invalid B{C{beta}}.
641 '''
642 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
644 @property_RO
645 def datum(self):
646 '''Get the C{datum} setting or the default (L{Datum}).
647 '''
648 return _xkwds_get(self._kwds, datum=self._datum)
650 def _datum_setter(self, datum):
651 '''(INTERNAL) Set the default C{datum}.
652 '''
653 d = datum or _xattr(self._knots[0], datum=None)
654 if d and d is not self._datum:
655 self._datum = _ellipsoidal_datum(d, name=self.name)
657 def _distances(self, x, y):
658 '''(INTERNAL) Yield distances to C{(x, y)}.
659 '''
660 _f, kwds = self._func, self._kwds
661 if not callable(_f): # PYCHOK no cover
662 self._notOverloaded(distance_function=_f)
663 try:
664 for i, k in enumerate(self._knots):
665 yield _f(y, x, k.lat, k.lon, **kwds)
666 except Exception as e:
667 i = Fmt.INDEX(knots=i)
668 raise HeightError(i, k, cause=e)
670 def _distancesTo(self, _To):
671 '''(INTERNAL) Yield distances C{_To}.
672 '''
673 try:
674 for i, k in enumerate(self._knots):
675 yield _To(k)
676 except Exception as e:
677 i = Fmt.INDEX(knots=i)
678 raise HeightError(i, k, cause=e)
680 def height(self, lats, lons, **wrap):
681 '''Interpolate the height for one or several lat-/longitudes.
683 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
684 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
685 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
686 and B{C{lons}} (C{bool}).
688 @return: A single interpolated height (C{float}) or a list of
689 interpolated heights (C{float}s).
691 @raise HeightError: Insufficient or non-matching number of B{C{lats}}
692 and B{C{lons}} or a L{pygeodesy.fidw} issue.
693 '''
694 lls = self._as_lls(lats, lons)
695 return self(lls, **wrap) # __call__(ll) or __call__(lls)
697 @Property_RO
698 def _heights(self):
699 '''(INTERNAL) Get the knots' heights.
700 '''
701 return tuple(_xattr(k, height=0) for k in self.knots)
703 def _hIDW(self, x, y):
704 '''(INTERNAL) Return the IDW-interpolated height at
705 location (x, y), both C{degrees} or C{radians}.
706 '''
707 ds, hs = self._distances(x, y), self._heights
708 try:
709 return _MODS.fmath.fidw(hs, ds, beta=self.beta)
710 except (TypeError, ValueError) as e:
711 raise HeightError(x=x, y=y, cause=e)
713 @property_RO
714 def hypot(self):
715 '''Get the C{hypot} setting (C{callable}) or C{None}.
716 '''
717 return _xkwds_get(self._kwds, hypot=None)
719 @property_RO
720 def knots(self):
721 '''Get the B{C{knots}} (C{list} or C{tuple}).
722 '''
723 return self._knots
725 @property_RO
726 def kwds(self):
727 '''Get the optional keyword arguments (C{dict}).
728 '''
729 return self._kwds
731 @property_RO
732 def limit(self):
733 '''Get the C{limit} setting (C{degrees}) or C{None}.
734 '''
735 return _xkwds_get(self._kwds, limit=None)
737 @property_RO
738 def radius(self):
739 '''Get the C{radius} setting (C{bool}) or C{None}.
740 '''
741 return _xkwds_get(self._kwds, radius=None)
743 @property_RO
744 def scaled(self):
745 '''Get the C{scaled} setting (C{bool}) or C{None}.
746 '''
747 return _xkwds_get(self._kwds, scaled=None)
749 @property_RO
750 def wrap(self):
751 '''Get the C{wrap} setting or the default (C{bool}) or C{None}.
752 '''
753 return _xkwds_get(self._kwds, wrap=self._wrap)
756class HeightIDWcosineAndoyerLambert(_HeightIDW):
757 '''Height interpolator using U{Inverse Distance Weighting
758 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
759 and the I{angular} distance in C{radians} from function
760 L{pygeodesy.cosineAndoyerLambert_}.
761 '''
762 def __init__(self, knots, beta=2, **name__datum_wrap):
763 '''New L{HeightIDWcosineAndoyerLambert} interpolator.
765 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height
766 interpolator (C{str}) and any keyword arguments for
767 function L{pygeodesy.cosineAndoyerLambert}.
769 @see: L{Here<_HeightIDW.__init__>} for further details.
770 '''
771 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
772 self._func = _formy.cosineAndoyerLambert
774 if _FOR_DOCS:
775 __call__ = _HeightIDW.__call__
776 height = _HeightIDW.height
779class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):
780 '''Height interpolator using U{Inverse Distance Weighting
781 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
782 and the I{angular} distance in C{radians} from function
783 L{pygeodesy.cosineForsytheAndoyerLambert_}.
784 '''
785 def __init__(self, knots, beta=2, **name__datum_wrap):
786 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
788 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height
789 interpolator (C{str}) and any keyword arguments for
790 function L{pygeodesy.cosineForsytheAndoyerLambert}.
792 @see: L{Here<_HeightIDW.__init__>} for further details.
793 '''
794 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
795 self._func = _formy.cosineForsytheAndoyerLambert
797 if _FOR_DOCS:
798 __call__ = _HeightIDW.__call__
799 height = _HeightIDW.height
802class HeightIDWcosineLaw(_HeightIDW):
803 '''Height interpolator using U{Inverse Distance Weighting
804 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
805 I{angular} distance in C{radians} from function L{pygeodesy.cosineLaw_}.
807 @note: See note at function L{pygeodesy.vincentys_}.
808 '''
809 def __init__(self, knots, beta=2, **name__radius_wrap):
810 '''New L{HeightIDWcosineLaw} interpolator.
812 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this height
813 interpolator (C{str}) and any keyword arguments for
814 function L{pygeodesy.cosineLaw}.
816 @see: L{Here<_HeightIDW.__init__>} for further details.
817 '''
818 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
819 self._func = _formy.cosineLaw
821 if _FOR_DOCS:
822 __call__ = _HeightIDW.__call__
823 height = _HeightIDW.height
826class HeightIDWdistanceTo(_HeightIDW):
827 '''Height interpolator using U{Inverse Distance Weighting
828 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
829 and the distance from the points' C{LatLon.distanceTo} method,
830 conventionally in C{meter}.
831 '''
832 def __init__(self, knots, beta=2, **name__distanceTo_kwds):
833 '''New L{HeightIDWdistanceTo} interpolator.
835 @kwarg name__distanceTo_kwds: Optional C{B{name}=NN} for this
836 height interpolator (C{str}) and keyword arguments
837 for B{C{knots}}' method C{LatLon.distanceTo}.
839 @see: L{Here<_HeightIDW.__init__>} for further details.
841 @note: All B{C{points}} I{must} be instances of the same
842 ellipsoidal or spherical C{LatLon} class, I{not
843 checked}.
844 '''
845 _HeightIDW.__init__(self, knots, beta=beta, **name__distanceTo_kwds)
846 ks0 = _distanceTo(HeightError, knots=self._knots)[0]
847 # use knots[0] class and datum to create compatible points
848 # in ._as_lls instead of class LatLon_ and datum None
849 self._datum = ks0.datum
850 self._LLis = ks0.classof # type(ks0)
852 def _distances(self, x, y):
853 '''(INTERNAL) Yield distances to C{(x, y)}.
854 '''
855 kwds, ll = self._kwds, self._LLis(y, x)
857 def _To(k):
858 return k.distanceTo(ll, **kwds)
860 return self._distancesTo(_To)
862 if _FOR_DOCS:
863 __call__ = _HeightIDW.__call__
864 height = _HeightIDW.height
867class HeightIDWequirectangular(_HeightIDW):
868 '''Height interpolator using U{Inverse Distance Weighting
869 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
870 the I{angular} distance in C{radians squared} like function
871 L{pygeodesy.equirectangular4}.
872 '''
873 def __init__(self, knots, beta=2, **name__adjust_limit_wrap): # XXX beta=1
874 '''New L{HeightIDWequirectangular} interpolator.
876 @kwarg name__adjust_limit_wrap: Optional C{B{name}=NN} for this
877 height interpolator (C{str}) and keyword arguments
878 for function L{pygeodesy.equirectangular4}.
880 @see: L{Here<_HeightIDW.__init__>} for further details.
881 '''
882 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_limit_wrap)
884 def _distances(self, x, y):
885 '''(INTERNAL) Yield distances to C{(x, y)}.
886 '''
887 _f, kwds = _formy.equirectangular4, self._kwds
889 def _To(k):
890 return _f(y, x, k.lat, k.lon, **kwds).distance2
892 return self._distancesTo(_To)
894 if _FOR_DOCS:
895 __call__ = _HeightIDW.__call__
896 height = _HeightIDW.height
899class HeightIDWeuclidean(_HeightIDW):
900 '''Height interpolator using U{Inverse Distance Weighting
901 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
902 I{angular} distance in C{radians} from function L{pygeodesy.euclidean_}.
903 '''
904 def __init__(self, knots, beta=2, **name__adjust_radius_wrap):
905 '''New L{HeightIDWeuclidean} interpolator.
907 @kwarg name__adjust_radius_wrap: Optional C{B{name}=NN} for this
908 height interpolator (C{str}) and keyword arguments
909 for function function L{pygeodesy.euclidean}.
911 @see: L{Here<_HeightIDW.__init__>} for further details.
912 '''
913 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_radius_wrap)
914 self._func = _formy.euclidean
916 if _FOR_DOCS:
917 __call__ = _HeightIDW.__call__
918 height = _HeightIDW.height
921class HeightIDWexact(_HeightIDW):
922 '''Height interpolator using U{Inverse Distance Weighting
923 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
924 and the I{angular} distance in C{degrees} from method
925 L{GeodesicExact.Inverse}.
926 '''
927 def __init__(self, knots, beta=2, datum=None, **name__wrap):
928 '''New L{HeightIDWexact} interpolator.
930 @kwarg datum: Datum to override the default C{Datums.WGS84} and
931 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
932 L{Ellipsoid2} or L{a_f2Tuple}).
933 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator
934 (C{str}) and a keyword argument for method C{Inverse1} of
935 class L{geodesicx.GeodesicExact}.
937 @raise TypeError: Invalid B{C{datum}}.
939 @see: L{Here<_HeightIDW.__init__>} for further details.
940 '''
941 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)
942 self._datum_setter(datum)
943 self._func = self.datum.ellipsoid.geodesicx.Inverse1
945 if _FOR_DOCS:
946 __call__ = _HeightIDW.__call__
947 height = _HeightIDW.height
950class HeightIDWflatLocal(_HeightIDW):
951 '''Height interpolator using U{Inverse Distance Weighting
952 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
953 and the I{angular} distance in C{radians squared} like function
954 L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}.
955 '''
956 def __init__(self, knots, beta=2, **name__datum_hypot_scaled_wrap):
957 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
959 @kwarg name__datum_hypot_scaled_wrap: Optional C{B{name}=NN}
960 for this height interpolator (C{str}) and any
961 keyword arguments for L{pygeodesy.flatLocal}.
963 @see: L{HeightIDW<_HeightIDW.__init__>} for further details.
964 '''
965 _HeightIDW.__init__(self, knots, beta=beta,
966 **name__datum_hypot_scaled_wrap)
967 self._func = _formy.flatLocal
969 if _FOR_DOCS:
970 __call__ = _HeightIDW.__call__
971 height = _HeightIDW.height
974class HeightIDWflatPolar(_HeightIDW):
975 '''Height interpolator using U{Inverse Distance Weighting
976 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
977 and the I{angular} distance in C{radians} from function
978 L{pygeodesy.flatPolar_}.
979 '''
980 def __init__(self, knots, beta=2, **name__radius_wrap):
981 '''New L{HeightIDWflatPolar} interpolator.
983 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
984 height interpolator (C{str}) and any keyword
985 arguments for function L{pygeodesy.flatPolar}.
987 @see: L{Here<_HeightIDW.__init__>} for further details.
988 '''
989 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
990 self._func = _formy.flatPolar
992 if _FOR_DOCS:
993 __call__ = _HeightIDW.__call__
994 height = _HeightIDW.height
997class HeightIDWhaversine(_HeightIDW):
998 '''Height interpolator using U{Inverse Distance Weighting
999 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1000 I{angular} distance in C{radians} from function L{pygeodesy.haversine_}.
1002 @note: See note at function L{pygeodesy.vincentys_}.
1003 '''
1004 def __init__(self, knots, beta=2, **name__radius_wrap):
1005 '''New L{HeightIDWhaversine} interpolator.
1007 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1008 height interpolator (C{str}) and any keyword
1009 arguments for function L{pygeodesy.haversine}.
1011 @see: L{Here<_HeightIDW.__init__>} for further details.
1012 '''
1013 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1014 self._func = _formy.haversine
1016 if _FOR_DOCS:
1017 __call__ = _HeightIDW.__call__
1018 height = _HeightIDW.height
1021class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny
1022 if _FOR_DOCS:
1023 __doc__ = HeightIDWflatLocal.__doc__
1024 __init__ = HeightIDWflatLocal.__init__
1025 __call__ = HeightIDWflatLocal.__call__
1026 height = HeightIDWflatLocal.height
1029class HeightIDWkarney(_HeightIDW):
1030 '''Height interpolator using U{Inverse Distance Weighting
1031 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1032 I{angular} distance in C{degrees} from I{Karney}'s U{geographiclib
1033 <https://PyPI.org/project/geographiclib>} method U{Geodesic.Inverse
1034 <https://GeographicLib.SourceForge.io/Python/doc/code.html#
1035 geographiclib.geodesic.Geodesic.Inverse>}.
1036 '''
1037 def __init__(self, knots, beta=2, datum=None, **name__wrap):
1038 '''New L{HeightIDWkarney} interpolator.
1040 @kwarg datum: Datum to override the default C{Datums.WGS84} and
1041 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
1042 L{Ellipsoid2} or L{a_f2Tuple}).
1043 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator
1044 (C{str}) and a keyword argument for method C{Inverse1} of
1045 class L{geodesicw.Geodesic}.
1047 @raise ImportError: Package U{geographiclib
1048 <https://PyPI.org/project/geographiclib>} missing.
1050 @raise TypeError: Invalid B{C{datum}}.
1052 @see: L{Here<_HeightIDW.__init__>} for further details.
1053 '''
1054 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)
1055 self._datum_setter(datum)
1056 self._func = self.datum.ellipsoid.geodesic.Inverse1
1058 if _FOR_DOCS:
1059 __call__ = _HeightIDW.__call__
1060 height = _HeightIDW.height
1063class HeightIDWthomas(_HeightIDW):
1064 '''Height interpolator using U{Inverse Distance Weighting
1065 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1066 I{angular} distance in C{radians} from function L{pygeodesy.thomas_}.
1067 '''
1068 def __init__(self, knots, beta=2, **name__datum_wrap):
1069 '''New L{HeightIDWthomas} interpolator.
1071 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this
1072 height interpolator (C{str}) and any keyword
1073 arguments for function L{pygeodesy.thomas}.
1075 @see: L{Here<_HeightIDW.__init__>} for further details.
1076 '''
1077 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
1078 self._func = _formy.thomas
1080 if _FOR_DOCS:
1081 __call__ = _HeightIDW.__call__
1082 height = _HeightIDW.height
1085class HeightIDWvincentys(_HeightIDW):
1086 '''Height interpolator using U{Inverse Distance Weighting
1087 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1088 I{angular} distance in C{radians} from function L{pygeodesy.vincentys_}.
1090 @note: See note at function L{pygeodesy.vincentys_}.
1091 '''
1092 def __init__(self, knots, beta=2, **name__radius_wrap):
1093 '''New L{HeightIDWvincentys} interpolator.
1095 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1096 height interpolator (C{str}) and any keyword
1097 arguments for function L{pygeodesy.vincentys}.
1099 @see: L{Here<_HeightIDW.__init__>} for further details.
1100 '''
1101 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1102 self._func = _formy.vincentys
1104 if _FOR_DOCS:
1105 __call__ = _HeightIDW.__call__
1106 height = _HeightIDW.height
1109__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightNamed)
1111# **) MIT License
1112#
1113# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1114#
1115# Permission is hereby granted, free of charge, to any person obtaining a
1116# copy of this software and associated documentation files (the "Software"),
1117# to deal in the Software without restriction, including without limitation
1118# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1119# and/or sell copies of the Software, and to permit persons to whom the
1120# Software is furnished to do so, subject to the following conditions:
1121#
1122# The above copyright notice and this permission notice shall be included
1123# in all copies or substantial portions of the Software.
1124#
1125# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1126# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1127# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1128# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1129# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1130# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1131# OTHER DEALINGS IN THE SOFTWARE.