Coverage for pygeodesy/heights.py: 96%
322 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-01 11:43 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-01 11:43 -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
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.05.25'
96_error_ = 'error'
97_llis_ = 'llis'
98_smoothing_ = 'smoothing'
101class HeightError(PointsError):
102 '''Height interpolator C{Height...} or interpolation issue.
103 '''
104 pass
107def _alist(ais):
108 # return list of floats, not numpy.float64s
109 return list(map(float, ais))
112def _ascalar(ais): # in .geoids
113 # return single float, not numpy.float64
114 ais = list(ais) # np.array, etc. to list
115 if len(ais) != 1:
116 n = Fmt.PAREN(len=repr(ais))
117 t = _SPACE_(len(ais), _NOTEQUAL_, 1)
118 raise _AssertionError(n, txt=t)
119 return float(ais[0]) # remove np.<type>
122def _atuple(ais):
123 # return tuple of floats, not numpy.float64s
124 return tuple(map(float, ais))
127def _as_llis2(llis, m=1, Error=HeightError): # in .geoids
128 # determine return type and convert lli C{LatLon}s to list
129 if not isinstance(llis, tuple): # llis are *args
130 n = Fmt.PAREN(type_=_STAR_(NN, _llis_))
131 raise _AssertionError(n, txt=repr(llis))
133 n = len(llis)
134 if n == 1: # convert single lli to 1-item list
135 llis = llis[0]
136 try:
137 n, llis = len2(llis)
138 _as = _alist # return list of interpolated heights
139 except TypeError: # single lli
140 n, llis = 1, [llis]
141 _as = _ascalar # return single interpolated heights
142 else: # of 0, 2 or more llis
143 _as = _atuple # return tuple of interpolated heights
145 if n < m:
146 raise _insufficientError(m, Error=Error, llis=n)
147 return _as, llis
150def _height_called(inst, lats, lons, Error=HeightError, **wrap): # in .geoids
151 LLis, d = inst._LLis, inst.datum
152 if _isDegrees(lats) and _isDegrees(lons):
153 llis = LLis(lats, lons, datum=d)
154 else:
155 n, lats = len2(lats)
156 m, lons = len2(lons)
157 if n != m:
158 # format a LenError, but raise an Error
159 e = LenError(inst.__class__, lats=n, lons=m, txt=None)
160 raise e if Error is LenError else Error(str(e))
161 llis = [LLis(*t, datum=d) for t in zip(lats, lons)]
162 return inst(llis, **wrap) # __call__(lli) or __call__(llis)
165def _insufficientError(need, Error=HeightError, **name_value): # PYCHOK no cover
166 # create an insufficient Error instance
167 t = _COMMASPACE_(_insufficient_, str(need) + _PLUS_)
168 return Error(txt=t, **name_value)
171def _ordedup(ts, lo=EPS, hi=PI2-EPS):
172 # clip, order and remove duplicates
173 # p, ks = 0, []
174 # for k in sorted(max(lo, min(hi, t)) for t in ts):
175 # if k > p:
176 # ks.append(k)
177 # p = k
178 # return ks
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 _np_sp = None # (numpy, scipy)
231 def __call__(self, *llis, **wrap): # PYCHOK no cover
232 '''Interpolate the height for one or several locations. I{Must be overloaded}.
234 @arg llis: One or more locations (C{LatLon}s), all positional.
235 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
236 locations (C{bool}), overriding the B{C{knots}}
237 setting.
239 @return: A single interpolated height (C{float}) or a list
240 or tuple of interpolated heights (C{float}s).
242 @raise HeightError: Insufficient number of B{C{llis}} or
243 an invalid B{C{lli}}.
245 @raise SciPyError: A C{scipy} issue.
247 @raise SciPyWarning: A C{scipy} warning as exception.
248 '''
249 self._notOverloaded(callername='__call__', *llis, **wrap)
251 def _as_xyllis4(self, llis, **wrap):
252 # convert lli C{LatLon}s to tuples or C{NumPy} arrays of
253 # C{SciPy} sphericals and determine the return type
254 atype = self.numpy.array
255 wrap = _xkwds(wrap, wrap=self._wrap)
256 _as, llis = _as_llis2(llis)
257 xis, yis, _ = zip(*_xyhs(llis=llis, **wrap)) # PYCHOK yield
258 return _as, atype(xis), atype(yis), llis
260 def _ev(self, *args): # PYCHOK no cover
261 '''(INTERNAL) I{Must be overloaded}.'''
262 self._notOverloaded(*args)
264 def _eval(self, llis, **wrap): # XXX single arg, not *args
265 _as, xis, yis, _ = self._as_xyllis4(llis, **wrap)
266 try: # SciPy .ev signature: y first, then x!
267 return _as(self._ev(yis, xis))
268 except Exception as x:
269 raise _SciPyIssue(x)
271 def height(self, lats, lons, **wrap): # PYCHOK no cover
272 '''Interpolate the height for one or several lat-/longitudes. I{Must be overloaded}.
274 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
275 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
276 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
277 and B{C{lons}} locationts (C{bool}), overriding
278 the B{C{knots}} setting.
280 @return: A single interpolated height (C{float}) or a list of
281 interpolated heights (C{float}s).
283 @raise HeightError: Insufficient or non-matching number of
284 B{C{lats}} and B{C{lons}}.
286 @raise SciPyError: A C{scipy} issue.
288 @raise SciPyWarning: A C{scipy} warning as exception.
289 '''
290 self._notOverloaded(lats, lons, **wrap)
292 def _np_sp2(self, throwarnings=False):
293 '''(INTERNAL) Import C{numpy} and C{scipy}, once.
294 '''
295 t = _HeightsBase._np_sp
296 if not t:
297 # raise SciPyWarnings, but not if
298 # scipy has already been imported
299 if throwarnings: # PYCHOK no cover
300 import sys
301 if _scipy_ not in sys.modules:
302 import warnings
303 warnings.filterwarnings(_error_)
305 sp = _xscipy(self.__class__, 1, 2)
306 np = _xnumpy(self.__class__, 1, 9)
308 _HeightsBase._np_sp = t = np, sp
309 return t
311 @Property_RO
312 def numpy(self):
313 '''Get the C{numpy} module or C{None}.
314 '''
315 np, _ = self._np_sp2()
316 return np
318 @Property_RO
319 def scipy(self):
320 '''Get the C{scipy} module or C{None}.
321 '''
322 _, sp = self._np_sp2()
323 return sp
325 @Property_RO
326 def scipy_interpolate(self):
327 '''Get the C{scipy.interpolate} module or C{None}.
328 '''
329 _ = self.scipy
330 import scipy.interpolate as spi # scipy 1.2.2
331 return spi
333 def _scipy_version(self, **n):
334 '''Get the C{scipy} version as 2- or 3-tuple C{(major, minor, micro)}.
335 '''
336 return _MODS.internals._version2(self.scipy.version.version, **n)
338 def _xyhs3(self, knots, wrap=False, **name):
339 # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals
340 xs, ys, hs = zip(*_xyhs(knots=knots, wrap=wrap)) # PYCHOK yield
341 n = len(hs)
342 if n < self.kmin:
343 raise _insufficientError(self.kmin, knots=n)
344 if name:
345 self.name = name
346 return map1(self.numpy.array, xs, ys, hs)
349class HeightCubic(_HeightsBase):
350 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
351 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
352 C{kind='cubic'}.
353 '''
354 _interp2d = None
355 _kind = _cubic_
356 _kmin = 16
358 def __init__(self, knots, **name_wrap):
359 '''New L{HeightCubic} interpolator.
361 @arg knots: The points with known height (C{LatLon}s).
362 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
363 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
364 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
365 C{True} (C{bool}).
367 @raise HeightError: Insufficient number of B{C{knots}} or
368 invalid B{C{knot}}.
370 @raise ImportError: Package C{numpy} or C{scipy} not found
371 or not installed.
373 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
375 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
376 as exception.
377 '''
378 spi = self.scipy_interpolate
380 xs, ys, hs = self._xyhs3(knots, **name_wrap)
381 try: # SciPy.interpolate.interp2d kind 'linear' or 'cubic'
382 self._interp2d = spi.interp2d(xs, ys, hs, kind=self._kind)
383 except Exception as x:
384 raise _SciPyIssue(x)
386 def __call__(self, *llis, **wrap):
387 '''Interpolate the height for one or several locations.
389 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
391 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
392 as exception.
394 @see: L{Here<_HeightsBase.__call__>} for further details.
395 '''
396 return _HeightsBase._eval(self, llis, **wrap)
398 def _ev(self, yis, xis): # PYCHOK expected
399 # to make SciPy .interp2d signature(x, y), single (x, y)
400 # match SciPy .ev signature(ys, xs), flipped multiples
401 return map(self._interp2d, xis, yis)
403 def height(self, lats, lons, **wrap):
404 '''Interpolate the height for one or several lat-/longitudes.
406 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
408 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
409 as exception.
411 @see: L{Here<_HeightsBase.height>} for further details.
412 '''
413 return _height_called(self, lats, lons, **wrap)
416class HeightLinear(HeightCubic):
417 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
418 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
419 C{kind='linear'}.
420 '''
421 _kind = _linear_
422 _kmin = 2
424 def __init__(self, knots, **name_wrap):
425 '''New L{HeightLinear} interpolator.
427 @see: L{Here<HeightCubic.__init__>} for all details.
428 '''
429 HeightCubic.__init__(self, knots, **name_wrap)
431 if _FOR_DOCS:
432 __call__ = HeightCubic.__call__
433 height = HeightCubic.height
436class HeightLSQBiSpline(_HeightsBase):
437 '''Height interpolator using C{SciPy} U{LSQSphereBivariateSpline
438 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
439 interpolate.LSQSphereBivariateSpline.html>}.
440 '''
441 _kmin = 16 # k = 3, always
443 def __init__(self, knots, weight=None, **name_wrap):
444 '''New L{HeightLSQBiSpline} interpolator.
446 @arg knots: The points with known height (C{LatLon}s).
447 @kwarg weight: Optional weight or weights for each B{C{knot}}
448 (C{scalar} or C{scalar}s).
449 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
450 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
451 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
452 C{True} (C{bool}).
454 @raise HeightError: Insufficient number of B{C{knots}} or
455 an invalid B{C{knot}} or B{C{weight}}.
457 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
459 @raise ImportError: Package C{numpy} or C{scipy} not found or
460 not installed.
462 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
464 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
465 warning as exception.
466 '''
467 np = self.numpy
468 spi = self.scipy_interpolate
470 xs, ys, hs = self._xyhs3(knots, **name_wrap)
471 n = len(hs)
473 w = weight
474 if isscalar(w):
475 w = float(w)
476 if w <= 0:
477 raise HeightError(weight=w)
478 w = [w] * n
479 elif w is not None:
480 m, w = len2(w)
481 if m != n:
482 raise LenError(HeightLSQBiSpline, weight=m, knots=n)
483 w = map2(float, w)
484 m = min(w)
485 if m <= 0: # PYCHOK no cover
486 w = Fmt.INDEX(weight=w.index(m))
487 raise HeightError(w, m)
488 try:
489 T = 1.0e-4 # like SciPy example
490 ps = np.array(_ordedup(xs, T, PI2 - T))
491 ts = np.array(_ordedup(ys, T, PI - T))
492 self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs,
493 ts, ps, eps=EPS, w=w).ev
494 except Exception as x:
495 raise _SciPyIssue(x)
497 def __call__(self, *llis, **wrap):
498 '''Interpolate the height for one or several locations.
500 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
502 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
503 warning as exception.
505 @see: L{Here<_HeightsBase.__call__>} for further details.
506 '''
507 return _HeightsBase._eval(self, llis, **wrap)
509 def height(self, lats, lons, **wrap):
510 '''Interpolate the height for one or several lat-/longitudes.
512 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
514 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
515 warning as exception.
517 @see: L{Here<_HeightsBase.height>} for further details.
518 '''
519 return _height_called(self, lats, lons, **wrap)
522class HeightSmoothBiSpline(_HeightsBase):
523 '''Height interpolator using C{SciPy} U{SmoothSphereBivariateSpline
524 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
525 interpolate.SmoothSphereBivariateSpline.html>}.
526 '''
527 _kmin = 16 # k = 3, always
529 def __init__(self, knots, s=4, **name_wrap):
530 '''New L{HeightSmoothBiSpline} interpolator.
532 @arg knots: The points with known height (C{LatLon}s).
533 @kwarg s: The spline smoothing factor (C{scalar}), default C{4}.
534 @kwarg name_wrap: Optional C{B{name}=NN} for this height interpolator
535 (C{str}) and keyword argument C{b{wrap}=False} to wrap or
536 I{normalize} all B{C{knots}} and B{C{llis}} locations iff
537 C{True} (C{bool}).
539 @raise HeightError: Insufficient number of B{C{knots}} or
540 an invalid B{C{knot}} or B{C{s}}.
542 @raise ImportError: Package C{numpy} or C{scipy} not found or not
543 installed.
545 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
547 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
548 warning as exception.
549 '''
550 spi = self.scipy_interpolate
552 s = Float_(s, name=_smoothing_, Error=HeightError, low=4)
554 xs, ys, hs = self._xyhs3(knots, **name_wrap)
555 try:
556 self._ev = spi.SmoothSphereBivariateSpline(ys, xs, hs,
557 eps=EPS, s=s).ev
558 except Exception as x:
559 raise _SciPyIssue(x)
561 def __call__(self, *llis, **wrap):
562 '''Interpolate the height for one or several locations.
564 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
566 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
567 warning as exception.
569 @see: L{Here<_HeightsBase.__call__>} for further details.
570 '''
571 return _HeightsBase._eval(self, llis, **wrap)
573 def height(self, lats, lons, **wrap):
574 '''Interpolate the height for one or several lat-/longitudes.
576 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
578 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
579 warning as exception.
581 @see: L{Here<_HeightsBase.height>} for further details.
582 '''
583 return _height_called(self, lats, lons, **wrap)
586class _HeightIDW(_HeightBase):
587 '''(INTERNAL) Base class for U{Inverse Distance Weighting
588 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) height
589 interpolators.
591 @see: U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
592 geostatistics/Inverse-Distance-Weighting/index.html>},
593 U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
594 shepard_interp_2d/shepard_interp_2d.html>} and other C{_HeightIDW*}
595 classes.
596 '''
597 _beta = 0 # fidw inverse power
598 _func = None # formy function
599 _knots = () # knots list or tuple
600 _kwds = {} # func_ options
602 def __init__(self, knots, beta=2, **name__kwds):
603 '''New C{_HeightIDW*} interpolator.
605 @arg knots: The points with known height (C{LatLon}s).
606 @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
607 @kwarg name__kwds: Optional C{B{name}=NN} for this height interpolator
608 (C{str}) and any keyword arguments for the distance function,
609 retrievable with property C{kwds}.
611 @raise HeightError: Insufficient number of B{C{knots}} or an invalid
612 B{C{knot}} or B{C{beta}}.
613 '''
614 name, kwds = _name2__(**name__kwds)
615 if name:
616 self.name = name
618 n, self._knots = len2(knots)
619 if n < self.kmin:
620 raise _insufficientError(self.kmin, knots=n)
621 self.beta = beta
622 self._kwds = kwds or {}
624 def __call__(self, *llis, **wrap):
625 '''Interpolate the height for one or several locations.
627 @arg llis: One or more locations (C{LatLon}s), all positional.
628 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
629 locations (C{bool}).
631 @return: A single interpolated height (C{float}) or a list
632 or tuple of interpolated heights (C{float}s).
634 @raise HeightError: Insufficient number of B{C{llis}}, an
635 invalid B{C{lli}} or L{pygeodesy.fidw}
636 issue.
637 '''
638 def _xy2(wrap=False):
639 _w = _Wrap._latlonop(wrap)
640 try: # like _xyhs above, but degrees
641 for i, ll in enumerate(llis):
642 yield _w(ll.lon, ll.lat)
643 except Exception as x:
644 i = Fmt.INDEX(llis=i)
645 raise HeightError(i, ll, cause=x)
647 _as, llis = _as_llis2(llis)
648 return _as(map(self._hIDW, *zip(*_xy2(**wrap))))
650 @property_RO
651 def adjust(self):
652 '''Get the C{adjust} setting (C{bool}) or C{None}.
653 '''
654 return _xkwds_get(self._kwds, adjust=None)
656 @property
657 def beta(self):
658 '''Get the inverse distance power (C{int}).
659 '''
660 return self._beta
662 @beta.setter # PYCHOK setter!
663 def beta(self, beta):
664 '''Set the inverse distance power (C{int} 1, 2, or 3).
666 @raise HeightError: Invalid B{C{beta}}.
667 '''
668 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
670 @property_RO
671 def datum(self):
672 '''Get the C{datum} setting or the default (L{Datum}).
673 '''
674 return _xkwds_get(self._kwds, datum=self._datum)
676 def _datum_setter(self, datum):
677 '''(INTERNAL) Set the default C{datum}.
678 '''
679 d = datum or _xattr(self._knots[0], datum=None)
680 if d and d is not self._datum:
681 self._datum = _ellipsoidal_datum(d, name=self.name)
683 def _distances(self, x, y):
684 '''(INTERNAL) Yield distances to C{(x, y)}.
685 '''
686 _f, kwds = self._func, self._kwds
687 if not callable(_f): # PYCHOK no cover
688 self._notOverloaded(distance_function=_f)
689 try:
690 for i, k in enumerate(self._knots):
691 yield _f(y, x, k.lat, k.lon, **kwds)
692 except Exception as e:
693 i = Fmt.INDEX(knots=i)
694 raise HeightError(i, k, cause=e)
696 @property_RO
697 def _fmath(self):
698 '''(INTERNAL) Get module C{fmath}, I{once}.
699 '''
700 _HeightIDW._fmath = f = _MODS.fmath # overwrite property_RO
701 return f
703 @property_RO
704 def _formy(self):
705 '''(INTERNAL) Get module C{formy}, I{once}.
706 '''
707 _HeightIDW._formy = f = _MODS.formy # overwrite property_RO
708 return f
710 def height(self, lats, lons, **wrap):
711 '''Interpolate the height for one or several lat-/longitudes.
713 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
714 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
715 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
716 and B{C{lons}} (C{bool}).
718 @return: A single interpolated height (C{float}) or a list of
719 interpolated heights (C{float}s).
721 @raise HeightError: Insufficient or non-matching number of B{C{lats}}
722 and B{C{lons}} or a L{pygeodesy.fidw} issue.
723 '''
724 return _height_called(self, lats, lons, **wrap)
726 @Property_RO
727 def _heights(self):
728 '''(INTERNAL) Get the knots' heights.
729 '''
730 return tuple(_xattr(k, height=0) for k in self.knots)
732 def _hIDW(self, x, y):
733 '''(INTERNAL) Return the IDW-interpolated height at
734 location (x, y), both C{degrees} or C{radians}.
735 '''
736 ds, hs = self._distances(x, y), self._heights
737 try:
738 return self._fmath.fidw(hs, ds, beta=self.beta)
739 except (TypeError, ValueError) as e:
740 raise HeightError(x=x, y=y, cause=e)
742 @property_RO
743 def hypot(self):
744 '''Get the C{hypot} setting (C{callable}) or C{None}.
745 '''
746 return _xkwds_get(self._kwds, hypot=None)
748 @property_RO
749 def knots(self):
750 '''Get the B{C{knots}} (C{list} or C{tuple}).
751 '''
752 return self._knots
754 @property_RO
755 def kwds(self):
756 '''Get the optional keyword arguments (C{dict}).
757 '''
758 return self._kwds
760 @property_RO
761 def limit(self):
762 '''Get the C{limit} setting (C{degrees}) or C{None}.
763 '''
764 return _xkwds_get(self._kwds, limit=None)
766 @property_RO
767 def radius(self):
768 '''Get the C{radius} setting (C{bool}) or C{None}.
769 '''
770 return _xkwds_get(self._kwds, radius=None)
772 @property_RO
773 def scaled(self):
774 '''Get the C{scaled} setting (C{bool}) or C{None}.
775 '''
776 return _xkwds_get(self._kwds, scaled=None)
778 @property_RO
779 def wrap(self):
780 '''Get the C{wrap} setting or the default (C{bool}) or C{None}.
781 '''
782 return _xkwds_get(self._kwds, wrap=self._wrap)
785class HeightIDWcosineAndoyerLambert(_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.cosineAndoyerLambert_}.
790 '''
791 def __init__(self, knots, beta=2, **name__datum_wrap):
792 '''New L{HeightIDWcosineAndoyerLambert} 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.cosineAndoyerLambert}.
798 @see: L{Here<_HeightIDW.__init__>} for further details.
799 '''
800 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
801 self._func = self._formy.cosineAndoyerLambert
803 if _FOR_DOCS:
804 __call__ = _HeightIDW.__call__
805 height = _HeightIDW.height
808class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):
809 '''Height interpolator using U{Inverse Distance Weighting
810 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
811 and the I{angular} distance in C{radians} from function
812 L{pygeodesy.cosineForsytheAndoyerLambert_}.
813 '''
814 def __init__(self, knots, beta=2, **name__datum_wrap):
815 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
817 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this height
818 interpolator (C{str}) and any keyword arguments for
819 function L{pygeodesy.cosineForsytheAndoyerLambert}.
821 @see: L{Here<_HeightIDW.__init__>} for further details.
822 '''
823 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
824 self._func = self._formy.cosineForsytheAndoyerLambert
826 if _FOR_DOCS:
827 __call__ = _HeightIDW.__call__
828 height = _HeightIDW.height
831class HeightIDWcosineLaw(_HeightIDW):
832 '''Height interpolator using U{Inverse Distance Weighting
833 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
834 I{angular} distance in C{radians} from function L{pygeodesy.cosineLaw_}.
836 @note: See note at function L{pygeodesy.vincentys_}.
837 '''
838 def __init__(self, knots, beta=2, **name__radius_wrap):
839 '''New L{HeightIDWcosineLaw} interpolator.
841 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this height
842 interpolator (C{str}) and any keyword arguments for
843 function L{pygeodesy.cosineLaw}.
845 @see: L{Here<_HeightIDW.__init__>} for further details.
846 '''
847 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
848 self._func = self._formy.cosineLaw
850 if _FOR_DOCS:
851 __call__ = _HeightIDW.__call__
852 height = _HeightIDW.height
855class HeightIDWdistanceTo(_HeightIDW):
856 '''Height interpolator using U{Inverse Distance Weighting
857 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
858 and the distance from the points' C{LatLon.distanceTo} method,
859 conventionally in C{meter}.
860 '''
861 def __init__(self, knots, beta=2, **name__distanceTo_kwds):
862 '''New L{HeightIDWdistanceTo} interpolator.
864 @kwarg name__distanceTo_kwds: Optional C{B{name}=NN} for this
865 height interpolator (C{str}) and keyword arguments
866 for B{C{knots}}' method C{LatLon.distanceTo}.
868 @see: L{Here<_HeightIDW.__init__>} for further details.
870 @note: All B{C{points}} I{must} be instances of the same
871 ellipsoidal or spherical C{LatLon} class, I{not
872 checked}.
873 '''
874 _HeightIDW.__init__(self, knots, beta=beta, **name__distanceTo_kwds)
875 ks0 = _distanceTo(HeightError, knots=self._knots)[0]
876 # use knots[0] class and datum to create compatible points
877 # in _height_called instead of class LatLon_ and datum None
878 self._datum = ks0.datum
879 self._LLis = ks0.classof # type(ks0)
881 def _distances(self, x, y):
882 '''(INTERNAL) Yield distances to C{(x, y)}.
883 '''
884 kwds, ll = self._kwds, self._LLis(y, x)
885 try:
886 for i, k in enumerate(self._knots):
887 yield k.distanceTo(ll, **kwds)
888 except Exception as e:
889 i = Fmt.INDEX(knots=i)
890 raise HeightError(i, k, cause=e)
892 if _FOR_DOCS:
893 __call__ = _HeightIDW.__call__
894 height = _HeightIDW.height
897class HeightIDWequirectangular(_HeightIDW):
898 '''Height interpolator using U{Inverse Distance Weighting
899 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
900 the I{angular} distance in C{radians squared} like function
901 L{pygeodesy.equirectangular4}.
902 '''
903 def __init__(self, knots, beta=2, **name__adjust_limit_wrap): # XXX beta=1
904 '''New L{HeightIDWequirectangular} interpolator.
906 @kwarg name__adjust_limit_wrap: Optional C{B{name}=NN} for this
907 height interpolator (C{str}) and keyword arguments
908 for function L{pygeodesy.equirectangular4}.
910 @see: L{Here<_HeightIDW.__init__>} for further details.
911 '''
912 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_limit_wrap)
914 def _distances(self, x, y):
915 '''(INTERNAL) Yield distances to C{(x, y)}.
916 '''
917 _f, kwds = self._formy.equirectangular4, self._kwds
918 try:
919 for i, k in enumerate(self._knots):
920 yield _f(y, x, k.lat, k.lon, **kwds).distance2
921 except Exception as e:
922 i = Fmt.INDEX(knots=i)
923 raise HeightError(i, k, cause=e)
925 if _FOR_DOCS:
926 __call__ = _HeightIDW.__call__
927 height = _HeightIDW.height
930class HeightIDWeuclidean(_HeightIDW):
931 '''Height interpolator using U{Inverse Distance Weighting
932 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
933 I{angular} distance in C{radians} from function L{pygeodesy.euclidean_}.
934 '''
935 def __init__(self, knots, beta=2, **name__adjust_radius_wrap):
936 '''New L{HeightIDWeuclidean} interpolator.
938 @kwarg name__adjust_radius_wrap: Optional C{B{name}=NN} for this
939 height interpolator (C{str}) and keyword arguments
940 for function function L{pygeodesy.euclidean}.
942 @see: L{Here<_HeightIDW.__init__>} for further details.
943 '''
944 _HeightIDW.__init__(self, knots, beta=beta, **name__adjust_radius_wrap)
945 self._func = self._formy.euclidean
947 if _FOR_DOCS:
948 __call__ = _HeightIDW.__call__
949 height = _HeightIDW.height
952class HeightIDWexact(_HeightIDW):
953 '''Height interpolator using U{Inverse Distance Weighting
954 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
955 and the I{angular} distance in C{degrees} from method
956 L{GeodesicExact.Inverse}.
957 '''
958 def __init__(self, knots, beta=2, datum=None, **name__wrap):
959 '''New L{HeightIDWexact} interpolator.
961 @kwarg datum: Datum to override the default C{Datums.WGS84} and
962 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
963 L{Ellipsoid2} or L{a_f2Tuple}).
964 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator
965 (C{str}) and a keyword argument for method C{Inverse1} of
966 class L{geodesicx.GeodesicExact}.
968 @raise TypeError: Invalid B{C{datum}}.
970 @see: L{Here<_HeightIDW.__init__>} for further details.
971 '''
972 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)
973 self._datum_setter(datum)
974 self._func = self.datum.ellipsoid.geodesicx.Inverse1
976 if _FOR_DOCS:
977 __call__ = _HeightIDW.__call__
978 height = _HeightIDW.height
981class HeightIDWflatLocal(_HeightIDW):
982 '''Height interpolator using U{Inverse Distance Weighting
983 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
984 and the I{angular} distance in C{radians squared} like function
985 L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}.
986 '''
987 def __init__(self, knots, beta=2, **name__datum_hypot_scaled_wrap):
988 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
990 @kwarg name__datum_hypot_scaled_wrap: Optional C{B{name}=NN}
991 for this height interpolator (C{str}) and any
992 keyword arguments for L{pygeodesy.flatLocal}.
994 @see: L{HeightIDW<_HeightIDW.__init__>} for further details.
995 '''
996 _HeightIDW.__init__(self, knots, beta=beta,
997 **name__datum_hypot_scaled_wrap)
998 self._func = self._formy.flatLocal
1000 if _FOR_DOCS:
1001 __call__ = _HeightIDW.__call__
1002 height = _HeightIDW.height
1005class HeightIDWflatPolar(_HeightIDW):
1006 '''Height interpolator using U{Inverse Distance Weighting
1007 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
1008 and the I{angular} distance in C{radians} from function
1009 L{pygeodesy.flatPolar_}.
1010 '''
1011 def __init__(self, knots, beta=2, **name__radius_wrap):
1012 '''New L{HeightIDWflatPolar} interpolator.
1014 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1015 height interpolator (C{str}) and any keyword
1016 arguments for function L{pygeodesy.flatPolar}.
1018 @see: L{Here<_HeightIDW.__init__>} for further details.
1019 '''
1020 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1021 self._func = self._formy.flatPolar
1023 if _FOR_DOCS:
1024 __call__ = _HeightIDW.__call__
1025 height = _HeightIDW.height
1028class HeightIDWhaversine(_HeightIDW):
1029 '''Height interpolator using U{Inverse Distance Weighting
1030 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1031 I{angular} distance in C{radians} from function L{pygeodesy.haversine_}.
1033 @note: See note at function L{pygeodesy.vincentys_}.
1034 '''
1035 def __init__(self, knots, beta=2, **name__radius_wrap):
1036 '''New L{HeightIDWhaversine} interpolator.
1038 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1039 height interpolator (C{str}) and any keyword
1040 arguments for function L{pygeodesy.haversine}.
1042 @see: L{Here<_HeightIDW.__init__>} for further details.
1043 '''
1044 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1045 self._func = self._formy.haversine
1047 if _FOR_DOCS:
1048 __call__ = _HeightIDW.__call__
1049 height = _HeightIDW.height
1052class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny
1053 if _FOR_DOCS:
1054 __doc__ = HeightIDWflatLocal.__doc__
1055 __init__ = HeightIDWflatLocal.__init__
1056 __call__ = HeightIDWflatLocal.__call__
1057 height = HeightIDWflatLocal.height
1060class HeightIDWkarney(_HeightIDW):
1061 '''Height interpolator using U{Inverse Distance Weighting
1062 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1063 I{angular} distance in C{degrees} from I{Karney}'s U{geographiclib
1064 <https://PyPI.org/project/geographiclib>} method U{Geodesic.Inverse
1065 <https://geographiclib.sourceforge.io/Python/doc/code.html#
1066 geographiclib.geodesic.Geodesic.Inverse>}.
1067 '''
1068 def __init__(self, knots, beta=2, datum=None, **name__wrap):
1069 '''New L{HeightIDWkarney} interpolator.
1071 @kwarg datum: Datum to override the default C{Datums.WGS84} and
1072 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
1073 L{Ellipsoid2} or L{a_f2Tuple}).
1074 @kwarg name__wrap: Optional C{B{name}=NN} for this height interpolator
1075 (C{str}) and a keyword argument for method C{Inverse1} of
1076 class L{geodesicw.Geodesic}.
1078 @raise ImportError: Package U{geographiclib
1079 <https://PyPI.org/project/geographiclib>} missing.
1081 @raise TypeError: Invalid B{C{datum}}.
1083 @see: L{Here<_HeightIDW.__init__>} for further details.
1084 '''
1085 _HeightIDW.__init__(self, knots, beta=beta, **name__wrap)
1086 self._datum_setter(datum)
1087 self._func = self.datum.ellipsoid.geodesic.Inverse1
1089 if _FOR_DOCS:
1090 __call__ = _HeightIDW.__call__
1091 height = _HeightIDW.height
1094class HeightIDWthomas(_HeightIDW):
1095 '''Height interpolator using U{Inverse Distance Weighting
1096 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1097 I{angular} distance in C{radians} from function L{pygeodesy.thomas_}.
1098 '''
1099 def __init__(self, knots, beta=2, **name__datum_wrap):
1100 '''New L{HeightIDWthomas} interpolator.
1102 @kwarg name__datum_wrap: Optional C{B{name}=NN} for this
1103 height interpolator (C{str}) and any keyword
1104 arguments for function L{pygeodesy.thomas}.
1106 @see: L{Here<_HeightIDW.__init__>} for further details.
1107 '''
1108 _HeightIDW.__init__(self, knots, beta=beta, **name__datum_wrap)
1109 self._func = self._formy.thomas
1111 if _FOR_DOCS:
1112 __call__ = _HeightIDW.__call__
1113 height = _HeightIDW.height
1116class HeightIDWvincentys(_HeightIDW):
1117 '''Height interpolator using U{Inverse Distance Weighting
1118 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1119 I{angular} distance in C{radians} from function L{pygeodesy.vincentys_}.
1121 @note: See note at function L{pygeodesy.vincentys_}.
1122 '''
1123 def __init__(self, knots, beta=2, **name__radius_wrap):
1124 '''New L{HeightIDWvincentys} interpolator.
1126 @kwarg name__radius_wrap: Optional C{B{name}=NN} for this
1127 height interpolator (C{str}) and any keyword
1128 arguments for function L{pygeodesy.vincentys}.
1130 @see: L{Here<_HeightIDW.__init__>} for further details.
1131 '''
1132 _HeightIDW.__init__(self, knots, beta=beta, **name__radius_wrap)
1133 self._func = self._formy.vincentys
1135 if _FOR_DOCS:
1136 __call__ = _HeightIDW.__call__
1137 height = _HeightIDW.height
1140__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightsBase)
1142# **) MIT License
1143#
1144# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1145#
1146# Permission is hereby granted, free of charge, to any person obtaining a
1147# copy of this software and associated documentation files (the "Software"),
1148# to deal in the Software without restriction, including without limitation
1149# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1150# and/or sell copies of the Software, and to permit persons to whom the
1151# Software is furnished to do so, subject to the following conditions:
1152#
1153# The above copyright notice and this permission notice shall be included
1154# in all copies or substantial portions of the Software.
1155#
1156# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1157# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1158# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1159# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1160# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1161# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1162# OTHER DEALINGS IN THE SOFTWARE.