Coverage for pygeodesy/heights.py: 96%
304 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-08-02 18:24 -0400
« prev ^ index » next coverage.py v7.6.0, created at 2024-08-02 18:24 -0400
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_, _cubic_, _insufficient_, _linear_, \
82 _NOTEQUAL_, _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.07.25'
96_error_ = 'error'
97_formy = _MODS.into(formy=__name__)
98_llis_ = 'llis'
99_smoothing_ = 'smoothing'
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 _height_called(inst, lats, lons, Error=HeightError, **wrap): # in .geoids
152 LLis, d = inst._LLis, inst.datum
153 if _isDegrees(lats) and _isDegrees(lons):
154 llis = LLis(lats, lons, datum=d)
155 else:
156 n, lats = len2(lats)
157 m, lons = len2(lons)
158 if n != m:
159 # format a LenError, but raise an Error
160 e = LenError(inst.__class__, lats=n, lons=m, txt=None)
161 raise e if Error is LenError else Error(str(e))
162 llis = [LLis(*t, datum=d) for t in zip(lats, lons)]
163 return inst(llis, **wrap) # __call__(lli) or __call__(llis)
166def _insufficientError(need, Error=HeightError, **name_value): # PYCHOK no cover
167 # create an insufficient Error instance
168 t = _COMMASPACE_(_insufficient_, str(need) + _PLUS_)
169 return Error(txt=t, **name_value)
172def _orderedup(ts, lo=EPS, hi=PI2-EPS):
173 # clip, order and remove duplicates
174 # p = 0
175 # for k in sorted(max(lo, min(hi, t)) for t in ts):
176 # if k > p:
177 # yield k
178 # p = k
179 return sorted(set(max(lo, min(hi, t)) for t in ts)) # list
182def _xyhs(wrap=False, **name_lls):
183 # map (lat, lon, h) to (x, y, h) in radians, offset
184 # x as 0 <= lon <= PI2 and y as 0 <= lat <= PI
185 name, lls = _xkwds_item2(name_lls)
186 _0, _90, _180 = _0_0, _90_0, _180_0
187 _m, _r, _w = max, radians, _Wrap._latlonop(wrap)
188 try:
189 for i, ll in enumerate(lls):
190 y, x = _w(ll.lat, ll.lon)
191 yield _m(_0, _r(x + _180)), \
192 _m(_0, _r(y + _90)), ll.height
193 except Exception as e:
194 i = Fmt.INDEX(name, i)
195 raise HeightError(i, ll, cause=e)
198class _HeightBase(_Named): # in .geoids
199 '''(INTERNAL) Interpolator base class.
200 '''
201 _datum = _WGS84 # default
202 _kmin = 2 # min number of knots
203 _LLis = LatLon_ # ._height class
204 _np_sp = None # (numpy, scipy)
205 _wrap = None # wrap knots and llis
207 @property_RO
208 def datum(self):
209 '''Get the C{datum} setting or the default (L{Datum}).
210 '''
211 return self._datum
213 @property_RO
214 def kmin(self):
215 '''Get the minimum number of knots (C{int}).
216 '''
217 return self._kmin
219 @property_RO
220 def wrap(self):
221 '''Get the C{wrap} setting (C{bool}) or C{None}.
222 '''
223 return self._wrap
226class _HeightsBase(_HeightBase): # in .geoids
227 '''(INTERNAL) Interpolator base class.
228 '''
229 def __call__(self, *llis, **wrap): # PYCHOK no cover
230 '''Interpolate the height for one or several locations. I{Must be overloaded}.
232 @arg llis: One or more locations (C{LatLon}s), all positional.
233 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
234 locations (C{bool}), overriding the B{C{knots}}
235 setting.
237 @return: A single interpolated height (C{float}) or a list
238 or tuple of interpolated heights (C{float}s).
240 @raise HeightError: Insufficient number of B{C{llis}} or
241 an invalid B{C{lli}}.
243 @raise SciPyError: A C{scipy} issue.
245 @raise SciPyWarning: A C{scipy} warning as exception.
246 '''
247 self._notOverloaded(callername='__call__', *llis, **wrap)
249 def _as_xyllis4(self, llis, **wrap):
250 # convert lli C{LatLon}s to tuples or C{NumPy} arrays of
251 # C{SciPy} sphericals and determine the return type
252 atype = self.numpy.array
253 wrap = _xkwds(wrap, wrap=self._wrap)
254 _as, llis = _as_llis2(llis)
255 xis, yis, _ = zip(*_xyhs(llis=llis, **wrap)) # PYCHOK yield
256 return _as, atype(xis), atype(yis), llis
258 def _ev(self, *args): # PYCHOK no cover
259 '''(INTERNAL) I{Must be overloaded}.'''
260 self._notOverloaded(*args)
262 def _eval(self, llis, **wrap): # XXX single arg, not *args
263 _as, xis, yis, _ = self._as_xyllis4(llis, **wrap)
264 try: # SciPy .ev signature: y first, then x!
265 return _as(self._ev(yis, xis))
266 except Exception as x:
267 raise _SciPyIssue(x)
269 def height(self, lats, lons, **wrap): # PYCHOK no cover
270 '''Interpolate the height for one or several lat-/longitudes. I{Must be overloaded}.
272 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
273 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
274 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
275 and B{C{lons}} locationts (C{bool}), overriding
276 the B{C{knots}} setting.
278 @return: A single interpolated height (C{float}) or a list of
279 interpolated heights (C{float}s).
281 @raise HeightError: Insufficient or non-matching number of
282 B{C{lats}} and B{C{lons}}.
284 @raise SciPyError: A C{scipy} issue.
286 @raise SciPyWarning: A C{scipy} warning as exception.
287 '''
288 self._notOverloaded(lats, lons, **wrap)
290 def _np_sp2(self, throwarnings=False): # PYCHOK no cover
291 '''(INTERNAL) Import C{numpy} and C{scipy}, once.
292 '''
293 # raise SciPyWarnings, but not if
294 # scipy has already been imported
295 if throwarnings: # PYCHOK no cover
296 import sys
297 if _scipy_ not in sys.modules:
298 import warnings
299 warnings.filterwarnings(_error_)
300 return self.numpy, self.scipy
302 @property_ROver
303 def numpy(self):
304 '''Get the C{numpy} module or C{None}.
305 '''
306 return _xnumpy(self.__class__, 1, 9) # overwrite property_ROver
308 @property_ROver
309 def scipy(self):
310 '''Get the C{scipy} module or C{None}.
311 '''
312 return _xscipy(self.__class__, 1, 2) # overwrite property_ROver
314 @property_ROver
315 def scipy_interpolate(self):
316 '''Get the C{scipy.interpolate} module or C{None}.
317 '''
318 _ = self.scipy
319 import scipy.interpolate as spi # scipy 1.2.2
320 return spi # overwrite property_ROver
322 def _scipy_version(self, **n):
323 '''Get the C{scipy} version as 2- or 3-tuple C{(major, minor, micro)}.
324 '''
325 return _MODS.internals._version2(self.scipy.version.version, **n)
327 def _xyhs3(self, knots, wrap=False, **name):
328 # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals
329 xs, ys, hs = zip(*_xyhs(knots=knots, wrap=wrap)) # PYCHOK yield
330 n = len(hs)
331 if n < self.kmin:
332 raise _insufficientError(self.kmin, knots=n)
333 if name:
334 self.name = name
335 return map1(self.numpy.array, xs, ys, hs)
338class HeightCubic(_HeightsBase):
339 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
340 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
341 C{kind='cubic'}.
342 '''
343 _interp2d = None
344 _kind = _cubic_
345 _kmin = 16
347 def __init__(self, knots, **name_wrap):
348 '''New L{HeightCubic} interpolator.
350 @arg knots: The points with known height (C{LatLon}s).
351 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
352 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
353 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
354 C{True} (C{bool}).
356 @raise HeightError: Insufficient number of B{C{knots}} or
357 invalid B{C{knot}}.
359 @raise ImportError: Package C{numpy} or C{scipy} not found
360 or not installed.
362 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
364 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
365 as exception.
366 '''
367 spi = self.scipy_interpolate
369 xs, ys, hs = self._xyhs3(knots, **name_wrap)
370 try: # SciPy.interpolate.interp2d kind 'linear' or 'cubic'
371 self._interp2d = spi.interp2d(xs, ys, hs, kind=self._kind)
372 except Exception as x:
373 raise _SciPyIssue(x)
375 def __call__(self, *llis, **wrap):
376 '''Interpolate the height for one or several locations.
378 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
380 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
381 as exception.
383 @see: L{Here<_HeightsBase.__call__>} for further details.
384 '''
385 return _HeightsBase._eval(self, llis, **wrap)
387 def _ev(self, yis, xis): # PYCHOK expected
388 # to make SciPy .interp2d signature(x, y), single (x, y)
389 # match SciPy .ev signature(ys, xs), flipped multiples
390 return map(self._interp2d, xis, yis)
392 def height(self, lats, lons, **wrap):
393 '''Interpolate the height for one or several lat-/longitudes.
395 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
397 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
398 as exception.
400 @see: L{Here<_HeightsBase.height>} for further details.
401 '''
402 return _height_called(self, lats, lons, **wrap)
405class HeightLinear(HeightCubic):
406 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
407 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
408 C{kind='linear'}.
409 '''
410 _kind = _linear_
411 _kmin = 2
413 def __init__(self, knots, **name_wrap):
414 '''New L{HeightLinear} interpolator.
416 @see: L{Here<HeightCubic.__init__>} for all details.
417 '''
418 HeightCubic.__init__(self, knots, **name_wrap)
420 if _FOR_DOCS:
421 __call__ = HeightCubic.__call__
422 height = HeightCubic.height
425class HeightLSQBiSpline(_HeightsBase):
426 '''Height interpolator using C{SciPy} U{LSQSphereBivariateSpline
427 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
428 interpolate.LSQSphereBivariateSpline.html>}.
429 '''
430 _kmin = 16 # k = 3, always
432 def __init__(self, knots, weight=None, low=1e-4, **name_wrap):
433 '''New L{HeightLSQBiSpline} interpolator.
435 @arg knots: The points with known height (C{LatLon}s).
436 @kwarg weight: Optional weight or weights for each B{C{knot}}
437 (C{scalar} or C{scalar}s).
438 @kwarg low: Optional lower bound for I{ordered knots} (C{radians}).
439 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
440 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
441 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
442 C{True} (C{bool}).
444 @raise HeightError: Insufficient number of B{C{knots}} or an invalid
445 B{C{knot}}, B{C{weight}} or B{C{eps}}.
447 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
449 @raise ImportError: Package C{numpy} or C{scipy} not found or
450 not installed.
452 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
454 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
455 warning as exception.
456 '''
457 np = self.numpy
458 spi = self.scipy_interpolate
460 xs, ys, hs = self._xyhs3(knots, **name_wrap)
461 n = len(hs)
463 w = weight
464 if isscalar(w):
465 w = float(w)
466 if w <= 0:
467 raise HeightError(weight=w)
468 w = [w] * n
469 elif w is not None:
470 m, w = len2(w)
471 if m != n:
472 raise LenError(HeightLSQBiSpline, weight=m, knots=n)
473 w = map2(float, w)
474 m = min(w)
475 if m <= 0: # PYCHOK no cover
476 w = Fmt.INDEX(weight=w.index(m))
477 raise HeightError(w, m)
478 try:
479 if not EPS < low < (PI_2 - EPS): # 1e-4 like SciPy example
480 raise HeightError(low=low)
481 ps = np.array(_orderedup(xs, low, PI2 - low))
482 ts = np.array(_orderedup(ys, low, PI - low))
483 self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs,
484 ts, ps, eps=EPS, w=w).ev
485 except Exception as x:
486 raise _SciPyIssue(x)
488 def __call__(self, *llis, **wrap):
489 '''Interpolate the height for one or several locations.
491 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
493 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
494 warning as exception.
496 @see: L{Here<_HeightsBase.__call__>} for further details.
497 '''
498 return _HeightsBase._eval(self, llis, **wrap)
500 def height(self, lats, lons, **wrap):
501 '''Interpolate the height for one or several lat-/longitudes.
503 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
505 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
506 warning as exception.
508 @see: L{Here<_HeightsBase.height>} for further details.
509 '''
510 return _height_called(self, lats, lons, **wrap)
513class HeightSmoothBiSpline(_HeightsBase):
514 '''Height interpolator using C{SciPy} U{SmoothSphereBivariateSpline
515 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
516 interpolate.SmoothSphereBivariateSpline.html>}.
517 '''
518 _kmin = 16 # k = 3, always
520 def __init__(self, knots, s=4, **name_wrap):
521 '''New L{HeightSmoothBiSpline} interpolator.
523 @arg knots: The points with known height (C{LatLon}s).
524 @kwarg s: The spline smoothing factor (C{scalar}), default C{4}.
525 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
526 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
527 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
528 C{True} (C{bool}).
530 @raise HeightError: Insufficient number of B{C{knots}} or
531 an invalid B{C{knot}} or B{C{s}}.
533 @raise ImportError: Package C{numpy} or C{scipy} not found or not
534 installed.
536 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
538 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
539 warning as exception.
540 '''
541 spi = self.scipy_interpolate
543 s = Float_(s, name=_smoothing_, 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)
552 def __call__(self, *llis, **wrap):
553 '''Interpolate the height for one or several locations.
555 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
557 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
558 warning as exception.
560 @see: L{Here<_HeightsBase.__call__>} for further details.
561 '''
562 return _HeightsBase._eval(self, llis, **wrap)
564 def height(self, lats, lons, **wrap):
565 '''Interpolate the height for one or several lat-/longitudes.
567 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
569 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
570 warning as exception.
572 @see: L{Here<_HeightsBase.height>} for further details.
573 '''
574 return _height_called(self, lats, lons, **wrap)
577class _HeightIDW(_HeightBase):
578 '''(INTERNAL) Base class for U{Inverse Distance Weighting
579 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) height
580 interpolators.
582 @see: U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
583 geostatistics/Inverse-Distance-Weighting/index.html>},
584 U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
585 shepard_interp_2d/shepard_interp_2d.html>} and other C{_HeightIDW*}
586 classes.
587 '''
588 _beta = 0 # fidw inverse power
589 _func = None # formy function
590 _knots = () # knots list or tuple
591 _kwds = {} # func_ options
593 def __init__(self, knots, beta=2, **name__kwds):
594 '''New C{_HeightIDW*} interpolator.
596 @arg knots: The points with known height (C{LatLon}s).
597 @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
598 @kwarg name__kwds: Optional C{B{name}=NN} for this height interpolator
599 (C{str}) and any keyword arguments for the distance function,
600 retrievable with property C{kwds}.
602 @raise HeightError: Insufficient number of B{C{knots}} or an invalid
603 B{C{knot}} or B{C{beta}}.
604 '''
605 name, kwds = _name2__(**name__kwds)
606 if name:
607 self.name = name
609 n, self._knots = len2(knots)
610 if n < self.kmin:
611 raise _insufficientError(self.kmin, knots=n)
612 self.beta = beta
613 self._kwds = kwds or {}
615 def __call__(self, *llis, **wrap):
616 '''Interpolate the height for one or several locations.
618 @arg llis: One or more locations (C{LatLon}s), all positional.
619 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
620 locations (C{bool}).
622 @return: A single interpolated height (C{float}) or a list
623 or tuple of interpolated heights (C{float}s).
625 @raise HeightError: Insufficient number of B{C{llis}}, an
626 invalid B{C{lli}} or L{pygeodesy.fidw}
627 issue.
628 '''
629 def _xy2(wrap=False):
630 _w = _Wrap._latlonop(wrap)
631 try: # like _xyhs above, but degrees
632 for i, ll in enumerate(llis):
633 yield _w(ll.lon, ll.lat)
634 except Exception as x:
635 i = Fmt.INDEX(llis=i)
636 raise HeightError(i, ll, cause=x)
638 _as, llis = _as_llis2(llis)
639 return _as(map(self._hIDW, *zip(*_xy2(**wrap))))
641 @property_RO
642 def adjust(self):
643 '''Get the C{adjust} setting (C{bool}) or C{None}.
644 '''
645 return _xkwds_get(self._kwds, adjust=None)
647 @property
648 def beta(self):
649 '''Get the inverse distance power (C{int}).
650 '''
651 return self._beta
653 @beta.setter # PYCHOK setter!
654 def beta(self, beta):
655 '''Set the inverse distance power (C{int} 1, 2, or 3).
657 @raise HeightError: Invalid B{C{beta}}.
658 '''
659 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
661 @property_RO
662 def datum(self):
663 '''Get the C{datum} setting or the default (L{Datum}).
664 '''
665 return _xkwds_get(self._kwds, datum=self._datum)
667 def _datum_setter(self, datum):
668 '''(INTERNAL) Set the default C{datum}.
669 '''
670 d = datum or _xattr(self._knots[0], datum=None)
671 if d and d is not self._datum:
672 self._datum = _ellipsoidal_datum(d, name=self.name)
674 def _distances(self, x, y):
675 '''(INTERNAL) Yield distances to C{(x, y)}.
676 '''
677 _f, kwds = self._func, self._kwds
678 if not callable(_f): # PYCHOK no cover
679 self._notOverloaded(distance_function=_f)
680 try:
681 for i, k in enumerate(self._knots):
682 yield _f(y, x, k.lat, k.lon, **kwds)
683 except Exception as e:
684 i = Fmt.INDEX(knots=i)
685 raise HeightError(i, k, cause=e)
687 def height(self, lats, lons, **wrap):
688 '''Interpolate the height for one or several lat-/longitudes.
690 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
691 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
692 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
693 and B{C{lons}} (C{bool}).
695 @return: A single interpolated height (C{float}) or a list of
696 interpolated heights (C{float}s).
698 @raise HeightError: Insufficient or non-matching number of B{C{lats}}
699 and B{C{lons}} or a L{pygeodesy.fidw} issue.
700 '''
701 return _height_called(self, lats, lons, **wrap)
703 @Property_RO
704 def _heights(self):
705 '''(INTERNAL) Get the knots' heights.
706 '''
707 return tuple(_xattr(k, height=0) for k in self.knots)
709 def _hIDW(self, x, y):
710 '''(INTERNAL) Return the IDW-interpolated height at
711 location (x, y), both C{degrees} or C{radians}.
712 '''
713 ds, hs = self._distances(x, y), self._heights
714 try:
715 return _MODS.fmath.fidw(hs, ds, beta=self.beta)
716 except (TypeError, ValueError) as e:
717 raise HeightError(x=x, y=y, cause=e)
719 @property_RO
720 def hypot(self):
721 '''Get the C{hypot} setting (C{callable}) or C{None}.
722 '''
723 return _xkwds_get(self._kwds, hypot=None)
725 @property_RO
726 def knots(self):
727 '''Get the B{C{knots}} (C{list} or C{tuple}).
728 '''
729 return self._knots
731 @property_RO
732 def kwds(self):
733 '''Get the optional keyword arguments (C{dict}).
734 '''
735 return self._kwds
737 @property_RO
738 def limit(self):
739 '''Get the C{limit} setting (C{degrees}) or C{None}.
740 '''
741 return _xkwds_get(self._kwds, limit=None)
743 @property_RO
744 def radius(self):
745 '''Get the C{radius} setting (C{bool}) or C{None}.
746 '''
747 return _xkwds_get(self._kwds, radius=None)
749 @property_RO
750 def scaled(self):
751 '''Get the C{scaled} setting (C{bool}) or C{None}.
752 '''
753 return _xkwds_get(self._kwds, scaled=None)
755 @property_RO
756 def wrap(self):
757 '''Get the C{wrap} setting or the default (C{bool}) or C{None}.
758 '''
759 return _xkwds_get(self._kwds, wrap=self._wrap)
762class HeightIDWcosineAndoyerLambert(_HeightIDW):
763 '''Height interpolator using U{Inverse Distance Weighting
764 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
765 and the I{angular} distance in C{radians} from function
766 L{pygeodesy.cosineAndoyerLambert_}.
767 '''
768 def __init__(self, knots, beta=2, **name__datum_wrap):
769 '''New L{HeightIDWcosineAndoyerLambert} interpolator.
771 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height
772 interpolator (C{str}) and any keyword arguments for
773 function L{pygeodesy.cosineAndoyerLambert}.
775 @see: L{Here<_HeightIDW.__init__>} for further details.
776 '''
777 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
778 self._func = _formy.cosineAndoyerLambert
780 if _FOR_DOCS:
781 __call__ = _HeightIDW.__call__
782 height = _HeightIDW.height
785class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):
786 '''Height interpolator using U{Inverse Distance Weighting
787 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
788 and the I{angular} distance in C{radians} from function
789 L{pygeodesy.cosineForsytheAndoyerLambert_}.
790 '''
791 def __init__(self, knots, beta=2, **name__datum_wrap):
792 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
794 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height
795 interpolator (C{str}) and any keyword arguments for
796 function L{pygeodesy.cosineForsytheAndoyerLambert}.
798 @see: L{Here<_HeightIDW.__init__>} for further details.
799 '''
800 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
801 self._func = _formy.cosineForsytheAndoyerLambert
803 if _FOR_DOCS:
804 __call__ = _HeightIDW.__call__
805 height = _HeightIDW.height
808class HeightIDWcosineLaw(_HeightIDW):
809 '''Height interpolator using U{Inverse Distance Weighting
810 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
811 I{angular} distance in C{radians} from function L{pygeodesy.cosineLaw_}.
813 @note: See note at function L{pygeodesy.vincentys_}.
814 '''
815 def __init__(self, knots, beta=2, **name__radius_wrap):
816 '''New L{HeightIDWcosineLaw} interpolator.
818 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this height
819 interpolator (C{str}) and any keyword arguments for
820 function L{pygeodesy.cosineLaw}.
822 @see: L{Here<_HeightIDW.__init__>} for further details.
823 '''
824 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
825 self._func = _formy.cosineLaw
827 if _FOR_DOCS:
828 __call__ = _HeightIDW.__call__
829 height = _HeightIDW.height
832class HeightIDWdistanceTo(_HeightIDW):
833 '''Height interpolator using U{Inverse Distance Weighting
834 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
835 and the distance from the points' C{LatLon.distanceTo} method,
836 conventionally in C{meter}.
837 '''
838 def __init__(self, knots, beta=2, **name__distanceTo_kwds):
839 '''New L{HeightIDWdistanceTo} interpolator.
841 @kwarg name__distanceTo_kwds: Optional C{B{name}=NN} for this
842 height interpolator (C{str}) and keyword arguments
843 for B{C{knots}}' method C{LatLon.distanceTo}.
845 @see: L{Here<_HeightIDW.__init__>} for further details.
847 @note: All B{C{points}} I{must} be instances of the same
848 ellipsoidal or spherical C{LatLon} class, I{not
849 checked}.
850 '''
851 _HeightIDW.__init__(self, knots, beta=beta, **name__distanceTo_kwds)
852 ks0 = _distanceTo(HeightError, knots=self._knots)[0]
853 # use knots[0] class and datum to create compatible points
854 # in _height_called instead of class LatLon_ and datum None
855 self._datum = ks0.datum
856 self._LLis = ks0.classof # type(ks0)
858 def _distances(self, x, y):
859 '''(INTERNAL) Yield distances to C{(x, y)}.
860 '''
861 kwds, ll = self._kwds, self._LLis(y, x)
862 try:
863 for i, k in enumerate(self._knots):
864 yield k.distanceTo(ll, **kwds)
865 except Exception as e:
866 i = Fmt.INDEX(knots=i)
867 raise HeightError(i, k, cause=e)
869 if _FOR_DOCS:
870 __call__ = _HeightIDW.__call__
871 height = _HeightIDW.height
874class HeightIDWequirectangular(_HeightIDW):
875 '''Height interpolator using U{Inverse Distance Weighting
876 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
877 the I{angular} distance in C{radians squared} like function
878 L{pygeodesy.equirectangular4}.
879 '''
880 def __init__(self, knots, beta=2, **name__adjust_limit_wrap): # XXX beta=1
881 '''New L{HeightIDWequirectangular} interpolator.
883 @kwarg name__adjust_limit_wrap: Optional C{B{name}=NN} for this
884 height interpolator (C{str}) and keyword arguments
885 for function L{pygeodesy.equirectangular4}.
887 @see: L{Here<_HeightIDW.__init__>} for further details.
888 '''
889 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_limit_wrap)
891 def _distances(self, x, y):
892 '''(INTERNAL) Yield distances to C{(x, y)}.
893 '''
894 _f, kwds = _formy.equirectangular4, self._kwds
895 try:
896 for i, k in enumerate(self._knots):
897 yield _f(y, x, k.lat, k.lon, **kwds).distance2
898 except Exception as e:
899 i = Fmt.INDEX(knots=i)
900 raise HeightError(i, k, cause=e)
902 if _FOR_DOCS:
903 __call__ = _HeightIDW.__call__
904 height = _HeightIDW.height
907class HeightIDWeuclidean(_HeightIDW):
908 '''Height interpolator using U{Inverse Distance Weighting
909 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
910 I{angular} distance in C{radians} from function L{pygeodesy.euclidean_}.
911 '''
912 def __init__(self, knots, beta=2, **name__adjust_radius_wrap):
913 '''New L{HeightIDWeuclidean} interpolator.
915 @kwarg name__adjust_radius_wrap: Optional C{B{name}=NN} for this
916 height interpolator (C{str}) and keyword arguments
917 for function function L{pygeodesy.euclidean}.
919 @see: L{Here<_HeightIDW.__init__>} for further details.
920 '''
921 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_radius_wrap)
922 self._func = _formy.euclidean
924 if _FOR_DOCS:
925 __call__ = _HeightIDW.__call__
926 height = _HeightIDW.height
929class HeightIDWexact(_HeightIDW):
930 '''Height interpolator using U{Inverse Distance Weighting
931 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
932 and the I{angular} distance in C{degrees} from method
933 L{GeodesicExact.Inverse}.
934 '''
935 def __init__(self, knots, beta=2, datum=None, **name__wrap):
936 '''New L{HeightIDWexact} interpolator.
938 @kwarg datum: Datum to override the default C{Datums.WGS84} and
939 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
940 L{Ellipsoid2} or L{a_f2Tuple}).
941 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator
942 (C{str}) and a keyword argument for method C{Inverse1} of
943 class L{geodesicx.GeodesicExact}.
945 @raise TypeError: Invalid B{C{datum}}.
947 @see: L{Here<_HeightIDW.__init__>} for further details.
948 '''
949 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)
950 self._datum_setter(datum)
951 self._func = self.datum.ellipsoid.geodesicx.Inverse1
953 if _FOR_DOCS:
954 __call__ = _HeightIDW.__call__
955 height = _HeightIDW.height
958class HeightIDWflatLocal(_HeightIDW):
959 '''Height interpolator using U{Inverse Distance Weighting
960 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
961 and the I{angular} distance in C{radians squared} like function
962 L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}.
963 '''
964 def __init__(self, knots, beta=2, **name__datum_hypot_scaled_wrap):
965 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
967 @kwarg name__datum_hypot_scaled_wrap: Optional C{B{name}=NN}
968 for this height interpolator (C{str}) and any
969 keyword arguments for L{pygeodesy.flatLocal}.
971 @see: L{HeightIDW<_HeightIDW.__init__>} for further details.
972 '''
973 _HeightIDW.__init__(self, knots, beta=beta,
974 **name__datum_hypot_scaled_wrap)
975 self._func = _formy.flatLocal
977 if _FOR_DOCS:
978 __call__ = _HeightIDW.__call__
979 height = _HeightIDW.height
982class HeightIDWflatPolar(_HeightIDW):
983 '''Height interpolator using U{Inverse Distance Weighting
984 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
985 and the I{angular} distance in C{radians} from function
986 L{pygeodesy.flatPolar_}.
987 '''
988 def __init__(self, knots, beta=2, **name__radius_wrap):
989 '''New L{HeightIDWflatPolar} interpolator.
991 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
992 height interpolator (C{str}) and any keyword
993 arguments for function L{pygeodesy.flatPolar}.
995 @see: L{Here<_HeightIDW.__init__>} for further details.
996 '''
997 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
998 self._func = _formy.flatPolar
1000 if _FOR_DOCS:
1001 __call__ = _HeightIDW.__call__
1002 height = _HeightIDW.height
1005class HeightIDWhaversine(_HeightIDW):
1006 '''Height interpolator using U{Inverse Distance Weighting
1007 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1008 I{angular} distance in C{radians} from function L{pygeodesy.haversine_}.
1010 @note: See note at function L{pygeodesy.vincentys_}.
1011 '''
1012 def __init__(self, knots, beta=2, **name__radius_wrap):
1013 '''New L{HeightIDWhaversine} interpolator.
1015 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1016 height interpolator (C{str}) and any keyword
1017 arguments for function L{pygeodesy.haversine}.
1019 @see: L{Here<_HeightIDW.__init__>} for further details.
1020 '''
1021 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1022 self._func = _formy.haversine
1024 if _FOR_DOCS:
1025 __call__ = _HeightIDW.__call__
1026 height = _HeightIDW.height
1029class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny
1030 if _FOR_DOCS:
1031 __doc__ = HeightIDWflatLocal.__doc__
1032 __init__ = HeightIDWflatLocal.__init__
1033 __call__ = HeightIDWflatLocal.__call__
1034 height = HeightIDWflatLocal.height
1037class HeightIDWkarney(_HeightIDW):
1038 '''Height interpolator using U{Inverse Distance Weighting
1039 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1040 I{angular} distance in C{degrees} from I{Karney}'s U{geographiclib
1041 <https://PyPI.org/project/geographiclib>} method U{Geodesic.Inverse
1042 <https://geographiclib.sourceforge.io/Python/doc/code.html#
1043 geographiclib.geodesic.Geodesic.Inverse>}.
1044 '''
1045 def __init__(self, knots, beta=2, datum=None, **name__wrap):
1046 '''New L{HeightIDWkarney} interpolator.
1048 @kwarg datum: Datum to override the default C{Datums.WGS84} and
1049 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
1050 L{Ellipsoid2} or L{a_f2Tuple}).
1051 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator
1052 (C{str}) and a keyword argument for method C{Inverse1} of
1053 class L{geodesicw.Geodesic}.
1055 @raise ImportError: Package U{geographiclib
1056 <https://PyPI.org/project/geographiclib>} missing.
1058 @raise TypeError: Invalid B{C{datum}}.
1060 @see: L{Here<_HeightIDW.__init__>} for further details.
1061 '''
1062 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)
1063 self._datum_setter(datum)
1064 self._func = self.datum.ellipsoid.geodesic.Inverse1
1066 if _FOR_DOCS:
1067 __call__ = _HeightIDW.__call__
1068 height = _HeightIDW.height
1071class HeightIDWthomas(_HeightIDW):
1072 '''Height interpolator using U{Inverse Distance Weighting
1073 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1074 I{angular} distance in C{radians} from function L{pygeodesy.thomas_}.
1075 '''
1076 def __init__(self, knots, beta=2, **name__datum_wrap):
1077 '''New L{HeightIDWthomas} interpolator.
1079 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this
1080 height interpolator (C{str}) and any keyword
1081 arguments for function L{pygeodesy.thomas}.
1083 @see: L{Here<_HeightIDW.__init__>} for further details.
1084 '''
1085 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
1086 self._func = _formy.thomas
1088 if _FOR_DOCS:
1089 __call__ = _HeightIDW.__call__
1090 height = _HeightIDW.height
1093class HeightIDWvincentys(_HeightIDW):
1094 '''Height interpolator using U{Inverse Distance Weighting
1095 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1096 I{angular} distance in C{radians} from function L{pygeodesy.vincentys_}.
1098 @note: See note at function L{pygeodesy.vincentys_}.
1099 '''
1100 def __init__(self, knots, beta=2, **name__radius_wrap):
1101 '''New L{HeightIDWvincentys} interpolator.
1103 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1104 height interpolator (C{str}) and any keyword
1105 arguments for function L{pygeodesy.vincentys}.
1107 @see: L{Here<_HeightIDW.__init__>} for further details.
1108 '''
1109 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1110 self._func = _formy.vincentys
1112 if _FOR_DOCS:
1113 __call__ = _HeightIDW.__call__
1114 height = _HeightIDW.height
1117__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightsBase)
1119# **) MIT License
1120#
1121# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1122#
1123# Permission is hereby granted, free of charge, to any person obtaining a
1124# copy of this software and associated documentation files (the "Software"),
1125# to deal in the Software without restriction, including without limitation
1126# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1127# and/or sell copies of the Software, and to permit persons to whom the
1128# Software is furnished to do so, subject to the following conditions:
1129#
1130# The above copyright notice and this permission notice shall be included
1131# in all copies or substantial portions of the Software.
1132#
1133# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1134# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1135# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1136# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1137# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1138# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1139# OTHER DEALINGS IN THE SOFTWARE.