Coverage for pygeodesy/heights.py: 96%
324 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-06 16:50 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-05-06 16:50 -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
76# from pygeodesy.fmath import fidw # _MODS
77# from pygeodesy.formy import cosineAndoyerLambert, cosineForsytheAndoyerLambert, \
78# cosineLaw, equirectangular_, euclidean, flatLocal, \
79# flatPolar, haversine, thomas, vincentys # _MODS
80from pygeodesy.interns import NN, _COMMASPACE_, _cubic_, _insufficient_, _knots_, \
81 _linear_, _NOTEQUAL_, _PLUS_, _scipy_, _SPACE_, \
82 _STAR_, _version2
83from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, _FOR_DOCS
84from pygeodesy.named import _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.04.07'
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(lls, wrap=False, name=_llis_):
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 _0, _90, _180 = _0_0, _90_0, _180_0
186 _m, _r, _w = max, radians, _Wrap._latlonop(wrap)
187 try:
188 for i, ll in enumerate(lls):
189 y, x = _w(ll.lat, ll.lon)
190 yield _m(_0, _r(x + _180)), \
191 _m(_0, _r(y + _90)), ll.height
192 except Exception as e:
193 i = Fmt.INDEX(name, i)
194 raise HeightError(i, ll, cause=e)
197class _HeightBase(_Named): # in .geoids
198 '''(INTERNAL) Interpolator base class.
199 '''
200 _datum = _WGS84 # default
201 _kmin = 2 # min number of knots
202 _LLis = LatLon_ # ._height class
203 _np_sp = None # (numpy, scipy)
204 _wrap = None # wrap knots and llis
206 @property_RO
207 def datum(self):
208 '''Get the C{datum} setting or the default (L{Datum}).
209 '''
210 return self._datum
212 @property_RO
213 def kmin(self):
214 '''Get the minimum number of knots (C{int}).
215 '''
216 return self._kmin
218 @property_RO
219 def wrap(self):
220 '''Get the C{wrap} setting (C{bool}) or C{None}.
221 '''
222 return self._wrap
225class _HeightsBase(_HeightBase): # in .geoids
226 '''(INTERNAL) Interpolator base class.
227 '''
228 _np_sp = None # (numpy, scipy)
230 def __call__(self, *llis, **wrap): # PYCHOK no cover
231 '''Interpolate the height for one or several locations. I{Must be overloaded}.
233 @arg llis: One or more locations (C{LatLon}s), all positional.
234 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
235 locations (C{bool}), overriding the B{C{knots}}
236 setting.
238 @return: A single interpolated height (C{float}) or a list
239 or tuple of interpolated heights (C{float}s).
241 @raise HeightError: Insufficient number of B{C{llis}} or
242 an invalid B{C{lli}}.
244 @raise SciPyError: A C{scipy} issue.
246 @raise SciPyWarning: A C{scipy} warning as exception.
247 '''
248 self._notOverloaded(callername='__call__', *llis, **wrap)
250 def _as_xyllis4(self, llis, **wrap):
251 # convert lli C{LatLon}s to tuples or C{NumPy} arrays of
252 # C{SciPy} sphericals and determine the return type
253 atype = self.numpy.array
254 wrap = _xkwds(wrap, wrap=self._wrap)
255 _as, llis = _as_llis2(llis)
256 xis, yis, _ = zip(*_xyhs(llis, **wrap)) # PYCHOK yield
257 return _as, atype(xis), atype(yis), llis
259 def _ev(self, *args): # PYCHOK no cover
260 '''(INTERNAL) I{Must be overloaded}.'''
261 self._notOverloaded(*args)
263 def _eval(self, llis, **wrap): # XXX single arg, not *args
264 _as, xis, yis, _ = self._as_xyllis4(llis, **wrap)
265 try: # SciPy .ev signature: y first, then x!
266 return _as(self._ev(yis, xis))
267 except Exception as x:
268 raise _SciPyIssue(x)
270 def height(self, lats, lons, **wrap): # PYCHOK no cover
271 '''Interpolate the height for one or several lat-/longitudes. I{Must be overloaded}.
273 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
274 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
275 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
276 and B{C{lons}} locationts (C{bool}), overriding
277 the B{C{knots}} setting.
279 @return: A single interpolated height (C{float}) or a list of
280 interpolated heights (C{float}s).
282 @raise HeightError: Insufficient or non-matching number of
283 B{C{lats}} and B{C{lons}}.
285 @raise SciPyError: A C{scipy} issue.
287 @raise SciPyWarning: A C{scipy} warning as exception.
288 '''
289 self._notOverloaded(lats, lons, **wrap)
291 def _np_sp2(self, throwarnings=False):
292 '''(INTERNAL) Import C{numpy} and C{scipy}, once.
293 '''
294 t = _HeightsBase._np_sp
295 if not t:
296 # raise SciPyWarnings, but not if
297 # scipy has already been imported
298 if throwarnings: # PYCHOK no cover
299 import sys
300 if _scipy_ not in sys.modules:
301 import warnings
302 warnings.filterwarnings(_error_)
304 sp = _xscipy(self.__class__, 1, 2)
305 np = _xnumpy(self.__class__, 1, 9)
307 _HeightsBase._np_sp = t = np, sp
308 return t
310 @Property_RO
311 def numpy(self):
312 '''Get the C{numpy} module or C{None}.
313 '''
314 np, _ = self._np_sp2()
315 return np
317 @Property_RO
318 def scipy(self):
319 '''Get the C{scipy} module or C{None}.
320 '''
321 _, sp = self._np_sp2()
322 return sp
324 @Property_RO
325 def scipy_interpolate(self):
326 '''Get the C{scipy.interpolate} module or C{None}.
327 '''
328 _ = self.scipy
329 import scipy.interpolate as spi # scipy 1.2.2
330 return spi
332 def _scipy_version(self, **n):
333 '''Get the C{scipy} version as 2- or 3-tuple C{(major, minor, micro)}.
334 '''
335 return _version2(self.scipy.version.version, **n)
337 def _xyhs3(self, knots, **wrap):
338 # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals
339 xs, ys, hs = zip(*_xyhs(knots, name=_knots_, **wrap)) # PYCHOK yield
340 n = len(hs)
341 if n < self.kmin:
342 raise _insufficientError(self.kmin, knots=n)
343 return map1(self.numpy.array, xs, ys, hs)
346class HeightCubic(_HeightsBase):
347 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
348 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
349 C{kind='cubic'}.
350 '''
351 _interp2d = None
352 _kind = _cubic_
353 _kmin = 16
355 def __init__(self, knots, name=NN, **wrap):
356 '''New L{HeightCubic} interpolator.
358 @arg knots: The points with known height (C{LatLon}s).
359 @kwarg name: Optional name for this height interpolator (C{str}).
360 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{knots}}
361 and B{C{llis}} locations (C{bool}).
363 @raise HeightError: Insufficient number of B{C{knots}} or
364 invalid B{C{knot}}.
366 @raise ImportError: Package C{numpy} or C{scipy} not found
367 or not installed.
369 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
371 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
372 as exception.
373 '''
374 spi = self.scipy_interpolate
376 xs, ys, hs = self._xyhs3(knots, **wrap)
377 try: # SciPy.interpolate.interp2d kind 'linear' or 'cubic'
378 self._interp2d = spi.interp2d(xs, ys, hs, kind=self._kind)
379 except Exception as x:
380 raise _SciPyIssue(x)
382 if name:
383 self.name = name
385 def __call__(self, *llis, **wrap):
386 '''Interpolate the height for one or several locations.
388 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
390 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
391 as exception.
393 @see: L{Here<_HeightsBase.__call__>} for other details.
394 '''
395 return _HeightsBase._eval(self, llis, **wrap)
397 def _ev(self, yis, xis): # PYCHOK expected
398 # to make SciPy .interp2d signature(x, y), single (x, y)
399 # match SciPy .ev signature(ys, xs), flipped multiples
400 return map(self._interp2d, xis, yis)
402 def height(self, lats, lons, **wrap):
403 '''Interpolate the height for one or several lat-/longitudes.
405 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
407 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
408 as exception.
410 @see: L{Here<_HeightsBase.height>} for other details.
411 '''
412 return _height_called(self, lats, lons, **wrap)
415class HeightLinear(HeightCubic):
416 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
417 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
418 C{kind='linear'}.
419 '''
420 _kind = _linear_
421 _kmin = 2
423 def __init__(self, knots, name=NN, **wrap):
424 '''New L{HeightLinear} interpolator.
426 @see: L{Here<HeightCubic.__init__>} for details.
427 '''
428 HeightCubic.__init__(self, knots, name=name, **wrap)
430 if _FOR_DOCS:
431 __call__ = HeightCubic.__call__
432 height = HeightCubic.height
435class HeightLSQBiSpline(_HeightsBase):
436 '''Height interpolator using C{SciPy} U{LSQSphereBivariateSpline
437 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
438 interpolate.LSQSphereBivariateSpline.html>}.
439 '''
440 _kmin = 16 # k = 3, always
442 def __init__(self, knots, weight=None, name=NN, **wrap):
443 '''New L{HeightLSQBiSpline} interpolator.
445 @arg knots: The points with known height (C{LatLon}s).
446 @kwarg weight: Optional weight or weights for each B{C{knot}}
447 (C{scalar} or C{scalar}s).
448 @kwarg name: Optional name for this height interpolator (C{str}).
449 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{knots}}
450 and B{C{llis}} locations (C{bool}).
452 @raise HeightError: Insufficient number of B{C{knots}} or
453 an invalid B{C{knot}} or B{C{weight}}.
455 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
457 @raise ImportError: Package C{numpy} or C{scipy} not found
458 or not installed.
460 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
462 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
463 warning as exception.
464 '''
465 np = self.numpy
466 spi = self.scipy_interpolate
468 xs, ys, hs = self._xyhs3(knots, **wrap)
469 n = len(hs)
471 w = weight
472 if isscalar(w):
473 w = float(w)
474 if w <= 0:
475 raise HeightError(weight=w)
476 w = [w] * n
477 elif w is not None:
478 m, w = len2(w)
479 if m != n:
480 raise LenError(HeightLSQBiSpline, weight=m, knots=n)
481 w = map2(float, w)
482 m = min(w)
483 if m <= 0: # PYCHOK no cover
484 w = Fmt.INDEX(weight=w.index(m))
485 raise HeightError(w, m)
486 try:
487 T = 1.0e-4 # like SciPy example
488 ps = np.array(_ordedup(xs, T, PI2 - T))
489 ts = np.array(_ordedup(ys, T, PI - T))
490 self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs,
491 ts, ps, eps=EPS, w=w).ev
492 except Exception as x:
493 raise _SciPyIssue(x)
495 if name:
496 self.name = name
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 other 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 other 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=NN, **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: Optional name for this height interpolator (C{str}).
536 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{knots}}
537 and any called B{C{llis}} and height B{C{lats}}
538 and B{C{lons}} locations (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
544 or not 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, **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 if name:
563 self.name = name
565 def __call__(self, *llis, **wrap):
566 '''Interpolate the height for one or several locations.
568 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
570 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
571 warning as exception.
573 @see: L{Here<_HeightsBase.__call__>} for other details.
574 '''
575 return _HeightsBase._eval(self, llis, **wrap)
577 def height(self, lats, lons, **wrap):
578 '''Interpolate the height for one or several lat-/longitudes.
580 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
582 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
583 warning as exception.
585 @see: L{Here<_HeightsBase.height>} for other details.
586 '''
587 return _height_called(self, lats, lons, **wrap)
590class _HeightIDW(_HeightBase):
591 '''(INTERNAL) Base class for U{Inverse Distance Weighting
592 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) height
593 interpolators.
595 @see: U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
596 geostatistics/Inverse-Distance-Weighting/index.html>},
597 U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
598 shepard_interp_2d/shepard_interp_2d.html>} and other C{_HeightIDW*}
599 classes.
600 '''
601 _beta = 0 # fidw inverse power
602 _func = None # formy function
603 _knots = () # knots list or tuple
604 _kwds = {} # func_ options
606 def __init__(self, knots, beta=2, name=NN, **kwds):
607 '''New C{_HeightIDW*} interpolator.
609 @arg knots: The points with known height (C{LatLon}s).
610 @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
611 @kwarg name: Optional name for this height interpolator (C{str}).
612 @kwarg kwds: Optional keyword argument for distance function,
613 retrievable with property C{kwds}.
615 @raise HeightError: Insufficient number of B{C{knots}} or
616 an invalid B{C{knot}} or B{C{beta}}.
617 '''
618 if name:
619 self.name = name
620 n, self._knots = len2(knots)
621 if n < self.kmin:
622 raise _insufficientError(self.kmin, knots=n)
623 self.beta = beta
624 self._kwds = kwds or {}
626 def __call__(self, *llis, **wrap):
627 '''Interpolate the height for one or several locations.
629 @arg llis: One or more locations (C{LatLon}s), all positional.
630 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
631 locations (C{bool}).
633 @return: A single interpolated height (C{float}) or a list
634 or tuple of interpolated heights (C{float}s).
636 @raise HeightError: Insufficient number of B{C{llis}}, an
637 invalid B{C{lli}} or L{pygeodesy.fidw}
638 issue.
639 '''
640 def _xy2(lls, wrap=False):
641 _w = _Wrap._latlonop(wrap)
642 try: # like _xyhs above, but degrees
643 for i, ll in enumerate(lls):
644 yield _w(ll.lon, ll.lat)
645 except Exception as e:
646 i = Fmt.INDEX(llis=i)
647 raise HeightError(i, ll, cause=e)
649 _as, llis = _as_llis2(llis)
650 return _as(map(self._hIDW, *zip(*_xy2(llis, **wrap))))
652 @property_RO
653 def adjust(self):
654 '''Get the C{adjust} setting (C{bool}) or C{None}.
655 '''
656 return _xkwds_get(self._kwds, adjust=None)
658 @property
659 def beta(self):
660 '''Get the inverse distance power (C{int}).
661 '''
662 return self._beta
664 @beta.setter # PYCHOK setter!
665 def beta(self, beta):
666 '''Set the inverse distance power (C{int} 1, 2, or 3).
668 @raise HeightError: Invalid B{C{beta}}.
669 '''
670 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
672 @property_RO
673 def datum(self):
674 '''Get the C{datum} setting or the default (L{Datum}).
675 '''
676 return _xkwds_get(self._kwds, datum=self._datum)
678 def _datum_setter(self, datum):
679 '''(INTERNAL) Set the default C{datum}.
680 '''
681 d = datum or _xattr(self._knots[0], datum=None)
682 if d and d is not self._datum:
683 self._datum = _ellipsoidal_datum(d, name=self.name)
685 def _distances(self, x, y):
686 '''(INTERNAL) Yield distances to C{(x, y)}.
687 '''
688 _f, kwds = self._func, self._kwds
689 if not callable(_f): # PYCHOK no cover
690 self._notOverloaded(distance_function=_f)
691 try:
692 for i, k in enumerate(self._knots):
693 yield _f(y, x, k.lat, k.lon, **kwds)
694 except Exception as e:
695 i = Fmt.INDEX(knots=i)
696 raise HeightError(i, k, cause=e)
698 @property_RO
699 def _fmath(self):
700 '''(INTERNAL) Get module C{fmath}, I{once}.
701 '''
702 _HeightIDW._fmath = f = _MODS.fmath # overwrite property_RO
703 return f
705 @property_RO
706 def _formy(self):
707 '''(INTERNAL) Get module C{formy}, I{once}.
708 '''
709 _HeightIDW._formy = f = _MODS.formy # overwrite property_RO
710 return f
712 def height(self, lats, lons, **wrap):
713 '''Interpolate the height for one or several lat-/longitudes.
715 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
716 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
717 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
718 and B{C{lons}} (C{bool}).
720 @return: A single interpolated height (C{float}) or a list of
721 interpolated heights (C{float}s).
723 @raise HeightError: Insufficient or non-matching number of B{C{lats}}
724 and B{C{lons}} or a L{pygeodesy.fidw} issue.
725 '''
726 return _height_called(self, lats, lons, **wrap)
728 @Property_RO
729 def _heights(self):
730 '''(INTERNAL) Get the knots' heights.
731 '''
732 return tuple(_xattr(k, height=0) for k in self.knots)
734 def _hIDW(self, x, y):
735 '''(INTERNAL) Return the IDW-interpolated height at
736 location (x, y), both C{degrees} or C{radians}.
737 '''
738 ds, hs = self._distances(x, y), self._heights
739 try:
740 return self._fmath.fidw(hs, ds, beta=self.beta)
741 except (TypeError, ValueError) as e:
742 raise HeightError(x=x, y=y, cause=e)
744 @property_RO
745 def hypot(self):
746 '''Get the C{hypot} setting (C{callable}) or C{None}.
747 '''
748 return _xkwds_get(self._kwds, hypot=None)
750 @property_RO
751 def knots(self):
752 '''Get the B{C{knots}} (C{list} or C{tuple}).
753 '''
754 return self._knots
756 @property_RO
757 def kwds(self):
758 '''Get the optional keyword arguments (C{dict}).
759 '''
760 return self._kwds
762 @property_RO
763 def limit(self):
764 '''Get the C{limit} setting (C{degrees}) or C{None}.
765 '''
766 return _xkwds_get(self._kwds, limit=None)
768 @property_RO
769 def radius(self):
770 '''Get the C{radius} setting (C{bool}) or C{None}.
771 '''
772 return _xkwds_get(self._kwds, radius=None)
774 @property_RO
775 def scaled(self):
776 '''Get the C{scaled} setting (C{bool}) or C{None}.
777 '''
778 return _xkwds_get(self._kwds, scaled=None)
780 @property_RO
781 def wrap(self):
782 '''Get the C{wrap} setting or the default (C{bool}) or C{None}.
783 '''
784 return _xkwds_get(self._kwds, wrap=self._wrap)
787class HeightIDWcosineAndoyerLambert(_HeightIDW):
788 '''Height interpolator using U{Inverse Distance Weighting
789 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
790 and the I{angular} distance in C{radians} from function
791 L{pygeodesy.cosineAndoyerLambert_}.
792 '''
793 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
794 '''New L{HeightIDWcosineAndoyerLambert} interpolator.
796 @kwarg datum_wrap: Optional keyword arguments for function
797 L{pygeodesy.cosineAndoyerLambert}.
799 @see: L{Here<_HeightIDW.__init__>} for other details.
800 '''
801 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)
802 self._func = self._formy.cosineAndoyerLambert
804 if _FOR_DOCS:
805 __call__ = _HeightIDW.__call__
806 height = _HeightIDW.height
809class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):
810 '''Height interpolator using U{Inverse Distance Weighting
811 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
812 and the I{angular} distance in C{radians} from function
813 L{pygeodesy.cosineForsytheAndoyerLambert_}.
814 '''
815 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
816 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
818 @kwarg datum_wrap: Optional keyword arguments for function
819 L{pygeodesy.cosineForsytheAndoyerLambert}.
821 @see: L{Here<_HeightIDW.__init__>} for other details.
822 '''
823 _HeightIDW.__init__(self, knots, beta=beta, name=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=NN, **radius_wrap):
839 '''New L{HeightIDWcosineLaw} interpolator.
841 @kwarg radius_wrap: Optional keyword arguments for function
842 L{pygeodesy.cosineLaw}.
844 @see: L{Here<_HeightIDW.__init__>} for other details.
845 '''
846 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
847 self._func = self._formy.cosineLaw
849 if _FOR_DOCS:
850 __call__ = _HeightIDW.__call__
851 height = _HeightIDW.height
854class HeightIDWdistanceTo(_HeightIDW):
855 '''Height interpolator using U{Inverse Distance Weighting
856 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
857 and the distance from the points' C{LatLon.distanceTo} method,
858 conventionally in C{meter}.
859 '''
860 def __init__(self, knots, beta=2, name=NN, **distanceTo_kwds):
861 '''New L{HeightIDWdistanceTo} interpolator.
863 @kwarg distanceTo_kwds: Optional keyword arguments for the
864 B{C{knots}}' C{LatLon.distanceTo}
865 method.
867 @see: L{Here<_HeightIDW.__init__>} for other details.
869 @note: All B{C{points}} I{must} be instances of the same
870 ellipsoidal or spherical C{LatLon} class, I{not
871 checked}.
872 '''
873 _HeightIDW.__init__(self, knots, beta=beta, name=name,
874 **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.equirectangular_}.
902 '''
903 def __init__(self, knots, beta=2, name=NN, **adjust_limit_wrap): # XXX beta=1
904 '''New L{HeightIDWequirectangular} interpolator.
906 @kwarg adjust_limit_wrap: Optional keyword arguments for
907 function L{pygeodesy.equirectangular_}.
909 @see: L{Here<_HeightIDW.__init__>} for other details.
910 '''
911 _HeightIDW.__init__(self, knots, beta=beta, name=name,
912 **adjust_limit_wrap)
914 def _distances(self, x, y):
915 '''(INTERNAL) Yield distances to C{(x, y)}.
916 '''
917 _f, kwds = self._formy.equirectangular_, 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=NN, **adjust_radius_wrap):
936 '''New L{HeightIDWeuclidean} interpolator.
938 @kwarg adjust_radius_wrap: Optional keyword arguments for
939 function L{pygeodesy.euclidean}.
941 @see: L{Here<_HeightIDW.__init__>} for other details.
942 '''
943 _HeightIDW.__init__(self, knots, beta=beta, name=name,
944 **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, name=NN, datum=None, **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 wrap: Optional keyword argument for method C{Inverse1}
965 of class L{geodesicx.GeodesicExact}.
967 @raise TypeError: Invalid B{C{datum}}.
969 @see: L{Here<_HeightIDW.__init__>} for other details.
970 '''
971 _HeightIDW.__init__(self, knots, beta=beta, name=name, **wrap)
972 self._datum_setter(datum)
973 self._func = self.datum.ellipsoid.geodesicx.Inverse1
975 if _FOR_DOCS:
976 __call__ = _HeightIDW.__call__
977 height = _HeightIDW.height
980class HeightIDWflatLocal(_HeightIDW):
981 '''Height interpolator using U{Inverse Distance Weighting
982 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
983 and the I{angular} distance in C{radians squared} like function
984 L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}.
985 '''
986 def __init__(self, knots, beta=2, name=NN, **datum_hypot_scaled_wrap):
987 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
989 @kwarg datum_hypot_scaled_wrap: Optional keyword arguments for
990 function L{pygeodesy.flatLocal}.
992 @see: L{HeightIDW<_HeightIDW.__init__>} for other details.
993 '''
994 _HeightIDW.__init__(self, knots, beta=beta, name=name,
995 **datum_hypot_scaled_wrap)
996 self._func = self._formy.flatLocal
998 if _FOR_DOCS:
999 __call__ = _HeightIDW.__call__
1000 height = _HeightIDW.height
1003class HeightIDWflatPolar(_HeightIDW):
1004 '''Height interpolator using U{Inverse Distance Weighting
1005 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
1006 and the I{angular} distance in C{radians} from function
1007 L{pygeodesy.flatPolar_}.
1008 '''
1009 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
1010 '''New L{HeightIDWflatPolar} interpolator.
1012 @kwarg radius_wrap: Optional keyword arguments for function
1013 L{pygeodesy.flatPolar}.
1015 @see: L{Here<_HeightIDW.__init__>} for other details.
1016 '''
1017 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1018 self._func = self._formy.flatPolar
1020 if _FOR_DOCS:
1021 __call__ = _HeightIDW.__call__
1022 height = _HeightIDW.height
1025class HeightIDWhaversine(_HeightIDW):
1026 '''Height interpolator using U{Inverse Distance Weighting
1027 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1028 I{angular} distance in C{radians} from function L{pygeodesy.haversine_}.
1030 @note: See note at function L{pygeodesy.vincentys_}.
1031 '''
1032 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
1033 '''New L{HeightIDWhaversine} interpolator.
1035 @kwarg radius_wrap: Optional keyword arguments for function
1036 L{pygeodesy.haversine}.
1038 @see: L{Here<_HeightIDW.__init__>} for other details.
1039 '''
1040 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1041 self._func = self._formy.haversine
1043 if _FOR_DOCS:
1044 __call__ = _HeightIDW.__call__
1045 height = _HeightIDW.height
1048class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny
1049 if _FOR_DOCS:
1050 __doc__ = HeightIDWflatLocal.__doc__
1051 __init__ = HeightIDWflatLocal.__init__
1052 __call__ = HeightIDWflatLocal.__call__
1053 height = HeightIDWflatLocal.height
1056class HeightIDWkarney(_HeightIDW):
1057 '''Height interpolator using U{Inverse Distance Weighting
1058 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1059 I{angular} distance in C{degrees} from I{Karney}'s U{geographiclib
1060 <https://PyPI.org/project/geographiclib>} method U{Geodesic.Inverse
1061 <https://geographiclib.sourceforge.io/Python/doc/code.html#
1062 geographiclib.geodesic.Geodesic.Inverse>}.
1063 '''
1064 def __init__(self, knots, beta=2, name=NN, datum=None, **wrap):
1065 '''New L{HeightIDWkarney} interpolator.
1067 @kwarg datum: Datum to override the default C{Datums.WGS84} and
1068 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
1069 L{Ellipsoid2} or L{a_f2Tuple}).
1070 @kwarg wrap: Optional keyword argument for method C{Inverse1}
1071 of class L{geodesicw.Geodesic}.
1073 @raise ImportError: Package U{geographiclib
1074 <https://PyPI.org/project/geographiclib>} missing.
1076 @raise TypeError: Invalid B{C{datum}}.
1078 @see: L{Here<_HeightIDW.__init__>} for other details.
1079 '''
1080 _HeightIDW.__init__(self, knots, beta=beta, name=name, **wrap)
1081 self._datum_setter(datum)
1082 self._func = self.datum.ellipsoid.geodesic.Inverse1
1084 if _FOR_DOCS:
1085 __call__ = _HeightIDW.__call__
1086 height = _HeightIDW.height
1089class HeightIDWthomas(_HeightIDW):
1090 '''Height interpolator using U{Inverse Distance Weighting
1091 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1092 I{angular} distance in C{radians} from function L{pygeodesy.thomas_}.
1093 '''
1094 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
1095 '''New L{HeightIDWthomas} interpolator.
1097 @kwarg datum_wrap: Optional keyword argument for function
1098 L{pygeodesy.thomas}.
1100 @see: L{Here<_HeightIDW.__init__>} for other details.
1101 '''
1102 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)
1103 self._func = self._formy.thomas
1105 if _FOR_DOCS:
1106 __call__ = _HeightIDW.__call__
1107 height = _HeightIDW.height
1110class HeightIDWvincentys(_HeightIDW):
1111 '''Height interpolator using U{Inverse Distance Weighting
1112 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1113 I{angular} distance in C{radians} from function L{pygeodesy.vincentys_}.
1115 @note: See note at function L{pygeodesy.vincentys_}.
1116 '''
1117 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
1118 '''New L{HeightIDWvincentys} interpolator.
1120 @kwarg radius_wrap: Optional keyword arguments for function
1121 L{pygeodesy.vincentys}.
1123 @see: L{Here<_HeightIDW.__init__>} for other details.
1124 '''
1125 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1126 self._func = self._formy.vincentys
1128 if _FOR_DOCS:
1129 __call__ = _HeightIDW.__call__
1130 height = _HeightIDW.height
1133__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightsBase)
1135# **) MIT License
1136#
1137# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1138#
1139# Permission is hereby granted, free of charge, to any person obtaining a
1140# copy of this software and associated documentation files (the "Software"),
1141# to deal in the Software without restriction, including without limitation
1142# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1143# and/or sell copies of the Software, and to permit persons to whom the
1144# Software is furnished to do so, subject to the following conditions:
1145#
1146# The above copyright notice and this permission notice shall be included
1147# in all copies or substantial portions of the Software.
1148#
1149# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1150# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1151# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1152# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1153# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1154# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1155# OTHER DEALINGS IN THE SOFTWARE.