Coverage for pygeodesy/heights.py: 96%
315 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
2# -*- coding: utf-8 -*-
4u'''Height interpolations at C{LatLon} points from 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, 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
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.06.03'
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 _ordedup(ts, lo=EPS, hi=PI2-EPS):
173 # clip, order and remove duplicates
174 # p, ks = 0, []
175 # for k in sorted(max(lo, min(hi, t)) for t in ts):
176 # if k > p:
177 # ks.append(k)
178 # p = k
179 # return ks
180 return sorted(set(max(lo, min(hi, t)) for t in ts)) # list
183def _xyhs(wrap=False, **name_lls):
184 # map (lat, lon, h) to (x, y, h) in radians, offset
185 # x as 0 <= lon <= PI2 and y as 0 <= lat <= PI
186 name, lls = _xkwds_item2(name_lls)
187 _0, _90, _180 = _0_0, _90_0, _180_0
188 _m, _r, _w = max, radians, _Wrap._latlonop(wrap)
189 try:
190 for i, ll in enumerate(lls):
191 y, x = _w(ll.lat, ll.lon)
192 yield _m(_0, _r(x + _180)), \
193 _m(_0, _r(y + _90)), ll.height
194 except Exception as e:
195 i = Fmt.INDEX(name, i)
196 raise HeightError(i, ll, cause=e)
199class _HeightBase(_Named): # in .geoids
200 '''(INTERNAL) Interpolator base class.
201 '''
202 _datum = _WGS84 # default
203 _kmin = 2 # min number of knots
204 _LLis = LatLon_ # ._height class
205 _np_sp = None # (numpy, scipy)
206 _wrap = None # wrap knots and llis
208 @property_RO
209 def datum(self):
210 '''Get the C{datum} setting or the default (L{Datum}).
211 '''
212 return self._datum
214 @property_RO
215 def kmin(self):
216 '''Get the minimum number of knots (C{int}).
217 '''
218 return self._kmin
220 @property_RO
221 def wrap(self):
222 '''Get the C{wrap} setting (C{bool}) or C{None}.
223 '''
224 return self._wrap
227class _HeightsBase(_HeightBase): # in .geoids
228 '''(INTERNAL) Interpolator base class.
229 '''
230 _np_sp = None # (numpy, scipy)
232 def __call__(self, *llis, **wrap): # PYCHOK no cover
233 '''Interpolate the height for one or several locations. I{Must be overloaded}.
235 @arg llis: One or more locations (C{LatLon}s), all positional.
236 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
237 locations (C{bool}), overriding the B{C{knots}}
238 setting.
240 @return: A single interpolated height (C{float}) or a list
241 or tuple of interpolated heights (C{float}s).
243 @raise HeightError: Insufficient number of B{C{llis}} or
244 an invalid B{C{lli}}.
246 @raise SciPyError: A C{scipy} issue.
248 @raise SciPyWarning: A C{scipy} warning as exception.
249 '''
250 self._notOverloaded(callername='__call__', *llis, **wrap)
252 def _as_xyllis4(self, llis, **wrap):
253 # convert lli C{LatLon}s to tuples or C{NumPy} arrays of
254 # C{SciPy} sphericals and determine the return type
255 atype = self.numpy.array
256 wrap = _xkwds(wrap, wrap=self._wrap)
257 _as, llis = _as_llis2(llis)
258 xis, yis, _ = zip(*_xyhs(llis=llis, **wrap)) # PYCHOK yield
259 return _as, atype(xis), atype(yis), llis
261 def _ev(self, *args): # PYCHOK no cover
262 '''(INTERNAL) I{Must be overloaded}.'''
263 self._notOverloaded(*args)
265 def _eval(self, llis, **wrap): # XXX single arg, not *args
266 _as, xis, yis, _ = self._as_xyllis4(llis, **wrap)
267 try: # SciPy .ev signature: y first, then x!
268 return _as(self._ev(yis, xis))
269 except Exception as x:
270 raise _SciPyIssue(x)
272 def height(self, lats, lons, **wrap): # PYCHOK no cover
273 '''Interpolate the height for one or several lat-/longitudes. I{Must be overloaded}.
275 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
276 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
277 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
278 and B{C{lons}} locationts (C{bool}), overriding
279 the B{C{knots}} setting.
281 @return: A single interpolated height (C{float}) or a list of
282 interpolated heights (C{float}s).
284 @raise HeightError: Insufficient or non-matching number of
285 B{C{lats}} and B{C{lons}}.
287 @raise SciPyError: A C{scipy} issue.
289 @raise SciPyWarning: A C{scipy} warning as exception.
290 '''
291 self._notOverloaded(lats, lons, **wrap)
293 def _np_sp2(self, throwarnings=False):
294 '''(INTERNAL) Import C{numpy} and C{scipy}, once.
295 '''
296 t = _HeightsBase._np_sp
297 if not t:
298 # raise SciPyWarnings, but not if
299 # scipy has already been imported
300 if throwarnings: # PYCHOK no cover
301 import sys
302 if _scipy_ not in sys.modules:
303 import warnings
304 warnings.filterwarnings(_error_)
306 sp = _xscipy(self.__class__, 1, 2)
307 np = _xnumpy(self.__class__, 1, 9)
309 _HeightsBase._np_sp = t = np, sp
310 return t
312 @Property_RO
313 def numpy(self):
314 '''Get the C{numpy} module or C{None}.
315 '''
316 np, _ = self._np_sp2()
317 return np
319 @Property_RO
320 def scipy(self):
321 '''Get the C{scipy} module or C{None}.
322 '''
323 _, sp = self._np_sp2()
324 return sp
326 @Property_RO
327 def scipy_interpolate(self):
328 '''Get the C{scipy.interpolate} module or C{None}.
329 '''
330 _ = self.scipy
331 import scipy.interpolate as spi # scipy 1.2.2
332 return spi
334 def _scipy_version(self, **n):
335 '''Get the C{scipy} version as 2- or 3-tuple C{(major, minor, micro)}.
336 '''
337 return _MODS.internals._version2(self.scipy.version.version, **n)
339 def _xyhs3(self, knots, wrap=False, **name):
340 # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals
341 xs, ys, hs = zip(*_xyhs(knots=knots, wrap=wrap)) # PYCHOK yield
342 n = len(hs)
343 if n < self.kmin:
344 raise _insufficientError(self.kmin, knots=n)
345 if name:
346 self.name = name
347 return map1(self.numpy.array, xs, ys, hs)
350class HeightCubic(_HeightsBase):
351 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
352 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
353 C{kind='cubic'}.
354 '''
355 _interp2d = None
356 _kind = _cubic_
357 _kmin = 16
359 def __init__(self, knots, **name_wrap):
360 '''New L{HeightCubic} interpolator.
362 @arg knots: The points with known height (C{LatLon}s).
363 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
364 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
365 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
366 C{True} (C{bool}).
368 @raise HeightError: Insufficient number of B{C{knots}} or
369 invalid B{C{knot}}.
371 @raise ImportError: Package C{numpy} or C{scipy} not found
372 or not installed.
374 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
376 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
377 as exception.
378 '''
379 spi = self.scipy_interpolate
381 xs, ys, hs = self._xyhs3(knots, **name_wrap)
382 try: # SciPy.interpolate.interp2d kind 'linear' or 'cubic'
383 self._interp2d = spi.interp2d(xs, ys, hs, kind=self._kind)
384 except Exception as x:
385 raise _SciPyIssue(x)
387 def __call__(self, *llis, **wrap):
388 '''Interpolate the height for one or several locations.
390 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
392 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
393 as exception.
395 @see: L{Here<_HeightsBase.__call__>} for further details.
396 '''
397 return _HeightsBase._eval(self, llis, **wrap)
399 def _ev(self, yis, xis): # PYCHOK expected
400 # to make SciPy .interp2d signature(x, y), single (x, y)
401 # match SciPy .ev signature(ys, xs), flipped multiples
402 return map(self._interp2d, xis, yis)
404 def height(self, lats, lons, **wrap):
405 '''Interpolate the height for one or several lat-/longitudes.
407 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
409 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
410 as exception.
412 @see: L{Here<_HeightsBase.height>} for further details.
413 '''
414 return _height_called(self, lats, lons, **wrap)
417class HeightLinear(HeightCubic):
418 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
419 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
420 C{kind='linear'}.
421 '''
422 _kind = _linear_
423 _kmin = 2
425 def __init__(self, knots, **name_wrap):
426 '''New L{HeightLinear} interpolator.
428 @see: L{Here<HeightCubic.__init__>} for all details.
429 '''
430 HeightCubic.__init__(self, knots, **name_wrap)
432 if _FOR_DOCS:
433 __call__ = HeightCubic.__call__
434 height = HeightCubic.height
437class HeightLSQBiSpline(_HeightsBase):
438 '''Height interpolator using C{SciPy} U{LSQSphereBivariateSpline
439 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
440 interpolate.LSQSphereBivariateSpline.html>}.
441 '''
442 _kmin = 16 # k = 3, always
444 def __init__(self, knots, weight=None, **name_wrap):
445 '''New L{HeightLSQBiSpline} interpolator.
447 @arg knots: The points with known height (C{LatLon}s).
448 @kwarg weight: Optional weight or weights for each B{C{knot}}
449 (C{scalar} or C{scalar}s).
450 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
451 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
452 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
453 C{True} (C{bool}).
455 @raise HeightError: Insufficient number of B{C{knots}} or
456 an invalid B{C{knot}} or B{C{weight}}.
458 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
460 @raise ImportError: Package C{numpy} or C{scipy} not found or
461 not installed.
463 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
465 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
466 warning as exception.
467 '''
468 np = self.numpy
469 spi = self.scipy_interpolate
471 xs, ys, hs = self._xyhs3(knots, **name_wrap)
472 n = len(hs)
474 w = weight
475 if isscalar(w):
476 w = float(w)
477 if w <= 0:
478 raise HeightError(weight=w)
479 w = [w] * n
480 elif w is not None:
481 m, w = len2(w)
482 if m != n:
483 raise LenError(HeightLSQBiSpline, weight=m, knots=n)
484 w = map2(float, w)
485 m = min(w)
486 if m <= 0: # PYCHOK no cover
487 w = Fmt.INDEX(weight=w.index(m))
488 raise HeightError(w, m)
489 try:
490 T = 1.0e-4 # like SciPy example
491 ps = np.array(_ordedup(xs, T, PI2 - T))
492 ts = np.array(_ordedup(ys, T, PI - T))
493 self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs,
494 ts, ps, eps=EPS, w=w).ev
495 except Exception as x:
496 raise _SciPyIssue(x)
498 def __call__(self, *llis, **wrap):
499 '''Interpolate the height for one or several locations.
501 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
503 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
504 warning as exception.
506 @see: L{Here<_HeightsBase.__call__>} for further details.
507 '''
508 return _HeightsBase._eval(self, llis, **wrap)
510 def height(self, lats, lons, **wrap):
511 '''Interpolate the height for one or several lat-/longitudes.
513 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
515 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
516 warning as exception.
518 @see: L{Here<_HeightsBase.height>} for further details.
519 '''
520 return _height_called(self, lats, lons, **wrap)
523class HeightSmoothBiSpline(_HeightsBase):
524 '''Height interpolator using C{SciPy} U{SmoothSphereBivariateSpline
525 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
526 interpolate.SmoothSphereBivariateSpline.html>}.
527 '''
528 _kmin = 16 # k = 3, always
530 def __init__(self, knots, s=4, **name_wrap):
531 '''New L{HeightSmoothBiSpline} interpolator.
533 @arg knots: The points with known height (C{LatLon}s).
534 @kwarg s: The spline smoothing factor (C{scalar}), default C{4}.
535 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
536 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
537 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
538 C{True} (C{bool}).
540 @raise HeightError: Insufficient number of B{C{knots}} or
541 an invalid B{C{knot}} or B{C{s}}.
543 @raise ImportError: Package C{numpy} or C{scipy} not found or not
544 installed.
546 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
548 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
549 warning as exception.
550 '''
551 spi = self.scipy_interpolate
553 s = Float_(s, name=_smoothing_, Error=HeightError, low=4)
555 xs, ys, hs = self._xyhs3(knots, **name_wrap)
556 try:
557 self._ev = spi.SmoothSphereBivariateSpline(ys, xs, hs,
558 eps=EPS, s=s).ev
559 except Exception as x:
560 raise _SciPyIssue(x)
562 def __call__(self, *llis, **wrap):
563 '''Interpolate the height for one or several locations.
565 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
567 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
568 warning as exception.
570 @see: L{Here<_HeightsBase.__call__>} for further details.
571 '''
572 return _HeightsBase._eval(self, llis, **wrap)
574 def height(self, lats, lons, **wrap):
575 '''Interpolate the height for one or several lat-/longitudes.
577 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
579 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
580 warning as exception.
582 @see: L{Here<_HeightsBase.height>} for further details.
583 '''
584 return _height_called(self, lats, lons, **wrap)
587class _HeightIDW(_HeightBase):
588 '''(INTERNAL) Base class for U{Inverse Distance Weighting
589 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) height
590 interpolators.
592 @see: U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
593 geostatistics/Inverse-Distance-Weighting/index.html>},
594 U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
595 shepard_interp_2d/shepard_interp_2d.html>} and other C{_HeightIDW*}
596 classes.
597 '''
598 _beta = 0 # fidw inverse power
599 _func = None # formy function
600 _knots = () # knots list or tuple
601 _kwds = {} # func_ options
603 def __init__(self, knots, beta=2, **name__kwds):
604 '''New C{_HeightIDW*} interpolator.
606 @arg knots: The points with known height (C{LatLon}s).
607 @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
608 @kwarg name__kwds: Optional C{B{name}=NN} for this height interpolator
609 (C{str}) and any keyword arguments for the distance function,
610 retrievable with property C{kwds}.
612 @raise HeightError: Insufficient number of B{C{knots}} or an invalid
613 B{C{knot}} or B{C{beta}}.
614 '''
615 name, kwds = _name2__(**name__kwds)
616 if name:
617 self.name = name
619 n, self._knots = len2(knots)
620 if n < self.kmin:
621 raise _insufficientError(self.kmin, knots=n)
622 self.beta = beta
623 self._kwds = kwds or {}
625 def __call__(self, *llis, **wrap):
626 '''Interpolate the height for one or several locations.
628 @arg llis: One or more locations (C{LatLon}s), all positional.
629 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
630 locations (C{bool}).
632 @return: A single interpolated height (C{float}) or a list
633 or tuple of interpolated heights (C{float}s).
635 @raise HeightError: Insufficient number of B{C{llis}}, an
636 invalid B{C{lli}} or L{pygeodesy.fidw}
637 issue.
638 '''
639 def _xy2(wrap=False):
640 _w = _Wrap._latlonop(wrap)
641 try: # like _xyhs above, but degrees
642 for i, ll in enumerate(llis):
643 yield _w(ll.lon, ll.lat)
644 except Exception as x:
645 i = Fmt.INDEX(llis=i)
646 raise HeightError(i, ll, cause=x)
648 _as, llis = _as_llis2(llis)
649 return _as(map(self._hIDW, *zip(*_xy2(**wrap))))
651 @property_RO
652 def adjust(self):
653 '''Get the C{adjust} setting (C{bool}) or C{None}.
654 '''
655 return _xkwds_get(self._kwds, adjust=None)
657 @property
658 def beta(self):
659 '''Get the inverse distance power (C{int}).
660 '''
661 return self._beta
663 @beta.setter # PYCHOK setter!
664 def beta(self, beta):
665 '''Set the inverse distance power (C{int} 1, 2, or 3).
667 @raise HeightError: Invalid B{C{beta}}.
668 '''
669 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
671 @property_RO
672 def datum(self):
673 '''Get the C{datum} setting or the default (L{Datum}).
674 '''
675 return _xkwds_get(self._kwds, datum=self._datum)
677 def _datum_setter(self, datum):
678 '''(INTERNAL) Set the default C{datum}.
679 '''
680 d = datum or _xattr(self._knots[0], datum=None)
681 if d and d is not self._datum:
682 self._datum = _ellipsoidal_datum(d, name=self.name)
684 def _distances(self, x, y):
685 '''(INTERNAL) Yield distances to C{(x, y)}.
686 '''
687 _f, kwds = self._func, self._kwds
688 if not callable(_f): # PYCHOK no cover
689 self._notOverloaded(distance_function=_f)
690 try:
691 for i, k in enumerate(self._knots):
692 yield _f(y, x, k.lat, k.lon, **kwds)
693 except Exception as e:
694 i = Fmt.INDEX(knots=i)
695 raise HeightError(i, k, cause=e)
697 def height(self, lats, lons, **wrap):
698 '''Interpolate the height for one or several lat-/longitudes.
700 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
701 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
702 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
703 and B{C{lons}} (C{bool}).
705 @return: A single interpolated height (C{float}) or a list of
706 interpolated heights (C{float}s).
708 @raise HeightError: Insufficient or non-matching number of B{C{lats}}
709 and B{C{lons}} or a L{pygeodesy.fidw} issue.
710 '''
711 return _height_called(self, lats, lons, **wrap)
713 @Property_RO
714 def _heights(self):
715 '''(INTERNAL) Get the knots' heights.
716 '''
717 return tuple(_xattr(k, height=0) for k in self.knots)
719 def _hIDW(self, x, y):
720 '''(INTERNAL) Return the IDW-interpolated height at
721 location (x, y), both C{degrees} or C{radians}.
722 '''
723 ds, hs = self._distances(x, y), self._heights
724 try:
725 return _MODS.fmath.fidw(hs, ds, beta=self.beta)
726 except (TypeError, ValueError) as e:
727 raise HeightError(x=x, y=y, cause=e)
729 @property_RO
730 def hypot(self):
731 '''Get the C{hypot} setting (C{callable}) or C{None}.
732 '''
733 return _xkwds_get(self._kwds, hypot=None)
735 @property_RO
736 def knots(self):
737 '''Get the B{C{knots}} (C{list} or C{tuple}).
738 '''
739 return self._knots
741 @property_RO
742 def kwds(self):
743 '''Get the optional keyword arguments (C{dict}).
744 '''
745 return self._kwds
747 @property_RO
748 def limit(self):
749 '''Get the C{limit} setting (C{degrees}) or C{None}.
750 '''
751 return _xkwds_get(self._kwds, limit=None)
753 @property_RO
754 def radius(self):
755 '''Get the C{radius} setting (C{bool}) or C{None}.
756 '''
757 return _xkwds_get(self._kwds, radius=None)
759 @property_RO
760 def scaled(self):
761 '''Get the C{scaled} setting (C{bool}) or C{None}.
762 '''
763 return _xkwds_get(self._kwds, scaled=None)
765 @property_RO
766 def wrap(self):
767 '''Get the C{wrap} setting or the default (C{bool}) or C{None}.
768 '''
769 return _xkwds_get(self._kwds, wrap=self._wrap)
772class HeightIDWcosineAndoyerLambert(_HeightIDW):
773 '''Height interpolator using U{Inverse Distance Weighting
774 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
775 and the I{angular} distance in C{radians} from function
776 L{pygeodesy.cosineAndoyerLambert_}.
777 '''
778 def __init__(self, knots, beta=2, **name__datum_wrap):
779 '''New L{HeightIDWcosineAndoyerLambert} interpolator.
781 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height
782 interpolator (C{str}) and any keyword arguments for
783 function L{pygeodesy.cosineAndoyerLambert}.
785 @see: L{Here<_HeightIDW.__init__>} for further details.
786 '''
787 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
788 self._func = _formy.cosineAndoyerLambert
790 if _FOR_DOCS:
791 __call__ = _HeightIDW.__call__
792 height = _HeightIDW.height
795class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):
796 '''Height interpolator using U{Inverse Distance Weighting
797 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
798 and the I{angular} distance in C{radians} from function
799 L{pygeodesy.cosineForsytheAndoyerLambert_}.
800 '''
801 def __init__(self, knots, beta=2, **name__datum_wrap):
802 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
804 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height
805 interpolator (C{str}) and any keyword arguments for
806 function L{pygeodesy.cosineForsytheAndoyerLambert}.
808 @see: L{Here<_HeightIDW.__init__>} for further details.
809 '''
810 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
811 self._func = _formy.cosineForsytheAndoyerLambert
813 if _FOR_DOCS:
814 __call__ = _HeightIDW.__call__
815 height = _HeightIDW.height
818class HeightIDWcosineLaw(_HeightIDW):
819 '''Height interpolator using U{Inverse Distance Weighting
820 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
821 I{angular} distance in C{radians} from function L{pygeodesy.cosineLaw_}.
823 @note: See note at function L{pygeodesy.vincentys_}.
824 '''
825 def __init__(self, knots, beta=2, **name__radius_wrap):
826 '''New L{HeightIDWcosineLaw} interpolator.
828 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this height
829 interpolator (C{str}) and any keyword arguments for
830 function L{pygeodesy.cosineLaw}.
832 @see: L{Here<_HeightIDW.__init__>} for further details.
833 '''
834 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
835 self._func = _formy.cosineLaw
837 if _FOR_DOCS:
838 __call__ = _HeightIDW.__call__
839 height = _HeightIDW.height
842class HeightIDWdistanceTo(_HeightIDW):
843 '''Height interpolator using U{Inverse Distance Weighting
844 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
845 and the distance from the points' C{LatLon.distanceTo} method,
846 conventionally in C{meter}.
847 '''
848 def __init__(self, knots, beta=2, **name__distanceTo_kwds):
849 '''New L{HeightIDWdistanceTo} interpolator.
851 @kwarg name__distanceTo_kwds: Optional C{B{name}=NN} for this
852 height interpolator (C{str}) and keyword arguments
853 for B{C{knots}}' method C{LatLon.distanceTo}.
855 @see: L{Here<_HeightIDW.__init__>} for further details.
857 @note: All B{C{points}} I{must} be instances of the same
858 ellipsoidal or spherical C{LatLon} class, I{not
859 checked}.
860 '''
861 _HeightIDW.__init__(self, knots, beta=beta, **name__distanceTo_kwds)
862 ks0 = _distanceTo(HeightError, knots=self._knots)[0]
863 # use knots[0] class and datum to create compatible points
864 # in _height_called instead of class LatLon_ and datum None
865 self._datum = ks0.datum
866 self._LLis = ks0.classof # type(ks0)
868 def _distances(self, x, y):
869 '''(INTERNAL) Yield distances to C{(x, y)}.
870 '''
871 kwds, ll = self._kwds, self._LLis(y, x)
872 try:
873 for i, k in enumerate(self._knots):
874 yield k.distanceTo(ll, **kwds)
875 except Exception as e:
876 i = Fmt.INDEX(knots=i)
877 raise HeightError(i, k, cause=e)
879 if _FOR_DOCS:
880 __call__ = _HeightIDW.__call__
881 height = _HeightIDW.height
884class HeightIDWequirectangular(_HeightIDW):
885 '''Height interpolator using U{Inverse Distance Weighting
886 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
887 the I{angular} distance in C{radians squared} like function
888 L{pygeodesy.equirectangular4}.
889 '''
890 def __init__(self, knots, beta=2, **name__adjust_limit_wrap): # XXX beta=1
891 '''New L{HeightIDWequirectangular} interpolator.
893 @kwarg name__adjust_limit_wrap: Optional C{B{name}=NN} for this
894 height interpolator (C{str}) and keyword arguments
895 for function L{pygeodesy.equirectangular4}.
897 @see: L{Here<_HeightIDW.__init__>} for further details.
898 '''
899 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_limit_wrap)
901 def _distances(self, x, y):
902 '''(INTERNAL) Yield distances to C{(x, y)}.
903 '''
904 _f, kwds = _formy.equirectangular4, self._kwds
905 try:
906 for i, k in enumerate(self._knots):
907 yield _f(y, x, k.lat, k.lon, **kwds).distance2
908 except Exception as e:
909 i = Fmt.INDEX(knots=i)
910 raise HeightError(i, k, cause=e)
912 if _FOR_DOCS:
913 __call__ = _HeightIDW.__call__
914 height = _HeightIDW.height
917class HeightIDWeuclidean(_HeightIDW):
918 '''Height interpolator using U{Inverse Distance Weighting
919 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
920 I{angular} distance in C{radians} from function L{pygeodesy.euclidean_}.
921 '''
922 def __init__(self, knots, beta=2, **name__adjust_radius_wrap):
923 '''New L{HeightIDWeuclidean} interpolator.
925 @kwarg name__adjust_radius_wrap: Optional C{B{name}=NN} for this
926 height interpolator (C{str}) and keyword arguments
927 for function function L{pygeodesy.euclidean}.
929 @see: L{Here<_HeightIDW.__init__>} for further details.
930 '''
931 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_radius_wrap)
932 self._func = _formy.euclidean
934 if _FOR_DOCS:
935 __call__ = _HeightIDW.__call__
936 height = _HeightIDW.height
939class HeightIDWexact(_HeightIDW):
940 '''Height interpolator using U{Inverse Distance Weighting
941 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
942 and the I{angular} distance in C{degrees} from method
943 L{GeodesicExact.Inverse}.
944 '''
945 def __init__(self, knots, beta=2, datum=None, **name__wrap):
946 '''New L{HeightIDWexact} interpolator.
948 @kwarg datum: Datum to override the default C{Datums.WGS84} and
949 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
950 L{Ellipsoid2} or L{a_f2Tuple}).
951 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator
952 (C{str}) and a keyword argument for method C{Inverse1} of
953 class L{geodesicx.GeodesicExact}.
955 @raise TypeError: Invalid B{C{datum}}.
957 @see: L{Here<_HeightIDW.__init__>} for further details.
958 '''
959 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)
960 self._datum_setter(datum)
961 self._func = self.datum.ellipsoid.geodesicx.Inverse1
963 if _FOR_DOCS:
964 __call__ = _HeightIDW.__call__
965 height = _HeightIDW.height
968class HeightIDWflatLocal(_HeightIDW):
969 '''Height interpolator using U{Inverse Distance Weighting
970 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
971 and the I{angular} distance in C{radians squared} like function
972 L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}.
973 '''
974 def __init__(self, knots, beta=2, **name__datum_hypot_scaled_wrap):
975 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
977 @kwarg name__datum_hypot_scaled_wrap: Optional C{B{name}=NN}
978 for this height interpolator (C{str}) and any
979 keyword arguments for L{pygeodesy.flatLocal}.
981 @see: L{HeightIDW<_HeightIDW.__init__>} for further details.
982 '''
983 _HeightIDW.__init__(self, knots, beta=beta,
984 **name__datum_hypot_scaled_wrap)
985 self._func = _formy.flatLocal
987 if _FOR_DOCS:
988 __call__ = _HeightIDW.__call__
989 height = _HeightIDW.height
992class HeightIDWflatPolar(_HeightIDW):
993 '''Height interpolator using U{Inverse Distance Weighting
994 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
995 and the I{angular} distance in C{radians} from function
996 L{pygeodesy.flatPolar_}.
997 '''
998 def __init__(self, knots, beta=2, **name__radius_wrap):
999 '''New L{HeightIDWflatPolar} interpolator.
1001 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1002 height interpolator (C{str}) and any keyword
1003 arguments for function L{pygeodesy.flatPolar}.
1005 @see: L{Here<_HeightIDW.__init__>} for further details.
1006 '''
1007 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1008 self._func = _formy.flatPolar
1010 if _FOR_DOCS:
1011 __call__ = _HeightIDW.__call__
1012 height = _HeightIDW.height
1015class HeightIDWhaversine(_HeightIDW):
1016 '''Height interpolator using U{Inverse Distance Weighting
1017 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1018 I{angular} distance in C{radians} from function L{pygeodesy.haversine_}.
1020 @note: See note at function L{pygeodesy.vincentys_}.
1021 '''
1022 def __init__(self, knots, beta=2, **name__radius_wrap):
1023 '''New L{HeightIDWhaversine} interpolator.
1025 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1026 height interpolator (C{str}) and any keyword
1027 arguments for function L{pygeodesy.haversine}.
1029 @see: L{Here<_HeightIDW.__init__>} for further details.
1030 '''
1031 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1032 self._func = _formy.haversine
1034 if _FOR_DOCS:
1035 __call__ = _HeightIDW.__call__
1036 height = _HeightIDW.height
1039class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny
1040 if _FOR_DOCS:
1041 __doc__ = HeightIDWflatLocal.__doc__
1042 __init__ = HeightIDWflatLocal.__init__
1043 __call__ = HeightIDWflatLocal.__call__
1044 height = HeightIDWflatLocal.height
1047class HeightIDWkarney(_HeightIDW):
1048 '''Height interpolator using U{Inverse Distance Weighting
1049 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1050 I{angular} distance in C{degrees} from I{Karney}'s U{geographiclib
1051 <https://PyPI.org/project/geographiclib>} method U{Geodesic.Inverse
1052 <https://geographiclib.sourceforge.io/Python/doc/code.html#
1053 geographiclib.geodesic.Geodesic.Inverse>}.
1054 '''
1055 def __init__(self, knots, beta=2, datum=None, **name__wrap):
1056 '''New L{HeightIDWkarney} interpolator.
1058 @kwarg datum: Datum to override the default C{Datums.WGS84} and
1059 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
1060 L{Ellipsoid2} or L{a_f2Tuple}).
1061 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator
1062 (C{str}) and a keyword argument for method C{Inverse1} of
1063 class L{geodesicw.Geodesic}.
1065 @raise ImportError: Package U{geographiclib
1066 <https://PyPI.org/project/geographiclib>} missing.
1068 @raise TypeError: Invalid B{C{datum}}.
1070 @see: L{Here<_HeightIDW.__init__>} for further details.
1071 '''
1072 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)
1073 self._datum_setter(datum)
1074 self._func = self.datum.ellipsoid.geodesic.Inverse1
1076 if _FOR_DOCS:
1077 __call__ = _HeightIDW.__call__
1078 height = _HeightIDW.height
1081class HeightIDWthomas(_HeightIDW):
1082 '''Height interpolator using U{Inverse Distance Weighting
1083 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1084 I{angular} distance in C{radians} from function L{pygeodesy.thomas_}.
1085 '''
1086 def __init__(self, knots, beta=2, **name__datum_wrap):
1087 '''New L{HeightIDWthomas} interpolator.
1089 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this
1090 height interpolator (C{str}) and any keyword
1091 arguments for function L{pygeodesy.thomas}.
1093 @see: L{Here<_HeightIDW.__init__>} for further details.
1094 '''
1095 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
1096 self._func = _formy.thomas
1098 if _FOR_DOCS:
1099 __call__ = _HeightIDW.__call__
1100 height = _HeightIDW.height
1103class HeightIDWvincentys(_HeightIDW):
1104 '''Height interpolator using U{Inverse Distance Weighting
1105 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1106 I{angular} distance in C{radians} from function L{pygeodesy.vincentys_}.
1108 @note: See note at function L{pygeodesy.vincentys_}.
1109 '''
1110 def __init__(self, knots, beta=2, **name__radius_wrap):
1111 '''New L{HeightIDWvincentys} interpolator.
1113 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1114 height interpolator (C{str}) and any keyword
1115 arguments for function L{pygeodesy.vincentys}.
1117 @see: L{Here<_HeightIDW.__init__>} for further details.
1118 '''
1119 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1120 self._func = _formy.vincentys
1122 if _FOR_DOCS:
1123 __call__ = _HeightIDW.__call__
1124 height = _HeightIDW.height
1127__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightsBase)
1129# **) MIT License
1130#
1131# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1132#
1133# Permission is hereby granted, free of charge, to any person obtaining a
1134# copy of this software and associated documentation files (the "Software"),
1135# to deal in the Software without restriction, including without limitation
1136# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1137# and/or sell copies of the Software, and to permit persons to whom the
1138# Software is furnished to do so, subject to the following conditions:
1139#
1140# The above copyright notice and this permission notice shall be included
1141# in all copies or substantial portions of the Software.
1142#
1143# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1144# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1145# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1146# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1147# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1148# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1149# OTHER DEALINGS IN THE SOFTWARE.