Coverage for pygeodesy/heights.py: 95%
310 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-06 15:27 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-06 15:27 -0400
2# -*- coding: utf-8 -*-
4u'''Height interpolations of C{LatLon} points.
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},
19called C{knots}. The C{knots} do not need to be ordered in any
20particular way.
222. Select one of the C{Height} classes for height interpolation
24C{>>> from pygeodesy import HeightCubic # or other Height... as HeightXyz}
263. Instantiate a height interpolator with the C{knots} and use keyword
27arguments to select different interpolation options
29C{>>> hinterpolator = HeightXyz(knots, **options)}
314. Get the interpolated height of other C{LatLon} location(s) with
33C{>>> ll = LatLon(1, 2, ...)}
34C{>>> h = hinterpolator(ll)}
36or
38C{>>> h0, h1, h2, ... = hinterpolator(ll0, ll1, ll2, ...)}
40or a list, tuple, generator, etc. of C{LatLon}s
42C{>>> hs = hinterpolator(lls)}
445. For separate lat- and longitudes invoke the C{.height} method
46C{>>> h = hinterpolator.height(lat, lon)}
48or as 2 lists, 2 tuples, etc.
50C{>>> hs = hinterpolator.height(lats, lons)}
52@note: Classes L{HeightCubic} and L{HeightLinear} require package U{numpy
53 <https://PyPI.org/project/numpy>}, classes L{HeightLSQBiSpline} and
54 L{HeightSmoothBiSpline} require package U{scipy<https://SciPy.org>}.
55 Classes L{HeightIDWkarney} and L{HeightIDWdistanceTo} -if used with
56 L{ellipsoidalKarney.LatLon} points- require I{Karney}'s U{geographiclib
57 <https://PyPI.org/project/geographiclib>} to be installed.
59@note: Errors from C{scipy} are raised as L{SciPyError}s. Warnings issued
60 by C{scipy} can be thrown as L{SciPyWarning} exceptions, provided
61 Python C{warnings} are filtered accordingly, see L{SciPyWarning}.
63@see: U{SciPy<https://docs.SciPy.org/doc/scipy/reference/interpolate.html>}
64 Interpolation.
65'''
66# make sure int/int division yields float quotient, see .basics
67from __future__ import division as _; del _ # PYCHOK semicolon
69from pygeodesy.basics import isscalar, len2, map1, map2, _xnumpy, _xscipy
70from pygeodesy.constants import EPS, PI, PI2, _0_0, _90_0, _180_0
71from pygeodesy.datums import _ellipsoidal_datum, _WGS84
72from pygeodesy.errors import _AssertionError, LenError, PointsError, \
73 _SciPyIssue, _xattr, _xkwds, _xkwds_get
74from pygeodesy.fmath import fidw, Float_, Int_
75from pygeodesy.formy import cosineAndoyerLambert, cosineForsytheAndoyerLambert, \
76 cosineLaw, equirectangular_, euclidean, flatLocal, \
77 flatPolar, haversine, thomas, vincentys, radians
78from pygeodesy.interns import NN, _COMMASPACE_, _cubic_, _insufficient_, _knots_, \
79 _linear_, _NOTEQUAL_, _PLUS_, _scipy_, _SPACE_, _STAR_
80from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS
81from pygeodesy.named import _Named, notOverloaded
82from pygeodesy.points import _distanceTo, LatLon_, Fmt, _Wrap
83from pygeodesy.props import Property_RO, property_RO
84# from pygeodesy.streprs import Fmt # from .points
85# from pygeodesy.units import Float_, Int_ # from .fmath
86# from pygeodesy.utily import _Wrap # from .points
88# from math import radians # from .formy
90__all__ = _ALL_LAZY.heights
91__version__ = '23.08.06'
93_error_ = 'error'
94_llis_ = 'llis'
95_smoothing_ = 'smoothing'
98class HeightError(PointsError):
99 '''Height interpolator C{Height...} or interpolation issue.
100 '''
101 pass
104def _alist(ais):
105 # return list of floats, not numpy.float64s
106 return list(map(float, ais))
109def _ascalar(ais): # in .geoids
110 # return single float, not numpy.float64
111 ais = list(ais) # np.array, etc. to list
112 if len(ais) != 1:
113 n = Fmt.PAREN(len=repr(ais))
114 t = _SPACE_(len(ais), _NOTEQUAL_, 1)
115 raise _AssertionError(n, txt=t)
116 return float(ais[0]) # remove np.<type>
119def _atuple(ais):
120 # return tuple of floats, not numpy.float64s
121 return tuple(map(float, ais))
124def _as_llis2(llis, m=1, Error=HeightError): # in .geoids
125 # determine return type and convert lli C{LatLon}s to list
126 if not isinstance(llis, tuple): # llis are *args
127 n = Fmt.PAREN(type_=_STAR_(NN, _llis_))
128 raise _AssertionError(n, txt=repr(llis))
130 n = len(llis)
131 if n == 1: # convert single lli to 1-item list
132 llis = llis[0]
133 try:
134 n, llis = len2(llis)
135 _as = _alist # return list of interpolated heights
136 except TypeError: # single lli
137 n, llis = 1, [llis]
138 _as = _ascalar # return single interpolated heights
139 else: # of 0, 2 or more llis
140 _as = _atuple # return tuple of interpolated heights
142 if n < m:
143 raise _insufficientError(m, Error=Error, llis=n)
144 return _as, llis
147def _height_called(inst, lats, lons, Error=HeightError, **wrap): # in .geoids
148 LLis, d = inst._LLis, inst.datum
149 if isscalar(lats) and isscalar(lons):
150 llis = LLis(lats, lons, datum=d)
151 else:
152 n, lats = len2(lats)
153 m, lons = len2(lons)
154 if n != m:
155 # format a LenError, but raise an Error
156 e = LenError(inst.__class__, lats=n, lons=m, txt=None)
157 raise e if Error is LenError else Error(str(e))
158 llis = [LLis(*t, datum=d) for t in zip(lats, lons)]
159 return inst(llis, **wrap) # __call__(lli) or __call__(llis)
162def _insufficientError(need, Error=HeightError, **name_value): # PYCHOK no cover
163 # create an insufficient Error instance
164 t = _COMMASPACE_(_insufficient_, str(need) + _PLUS_)
165 return Error(txt=t, **name_value)
168def _ordedup(ts, lo=EPS, hi=PI2-EPS):
169 # clip, order and remove duplicates
170 # p, ks = 0, []
171 # for k in sorted(max(lo, min(hi, t)) for t in ts):
172 # if k > p:
173 # ks.append(k)
174 # p = k
175 # return ks
176 return sorted(set(max(lo, min(hi, t)) for t in ts)) # list
179def _xyhs(lls, wrap=False, name=_llis_):
180 # map (lat, lon, h) to (x, y, h) in radians, offset
181 # x as 0 <= lon <= PI2 and y as 0 <= lat <= PI
182 _0, _90, _180 = _0_0, _90_0, _180_0
183 _m, _r, _w = max, radians, _Wrap._latlonop(wrap)
184 try:
185 for i, ll in enumerate(lls):
186 y, x = _w(ll.lat, ll.lon)
187 yield _m(_0, _r(x + _180)), \
188 _m(_0, _r(y + _90)), ll.height
189 except AttributeError as x:
190 i = Fmt.SQUARE(name, i)
191 raise HeightError(i, ll, cause=x)
194class _HeightBase(_Named): # in .geoids
195 '''(INTERNAL) Interpolator base class.
196 '''
197 _datum = _WGS84 # default
198 _kmin = 2 # min number of knots
199 _LLis = LatLon_ # ._height class
200 _np_sp = None # (numpy, scipy)
201 _wrap = None # wrap knots and llis
203 @property_RO
204 def datum(self):
205 '''Get the C{datum} setting or the default (L{Datum}).
206 '''
207 return self._datum
209 @property_RO
210 def kmin(self):
211 '''Get the minimum number of knots (C{int}).
212 '''
213 return self._kmin
215 @property_RO
216 def wrap(self):
217 '''Get the C{wrap} setting (C{bool}) or C{None}.
218 '''
219 return self._wrap
222class _HeightsBase(_HeightBase): # in .geoids
223 '''(INTERNAL) Interpolator base class.
224 '''
225 _np_sp = None # (numpy, scipy)
227 def __call__(self, *llis, **wrap): # PYCHOK no cover
228 '''Interpolate the height for one or several locations.
230 @arg llis: The location or locations (C{LatLon}, ... or
231 C{LatLon}s).
232 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
233 locations (C{bool}), overriding the B{C{knots}}
234 setting.
236 @return: A single interpolated height (C{float}) or a list
237 or tuple of interpolated heights (C{float}s).
239 @raise HeightError: Insufficient number of B{C{llis}} or
240 an invalid B{C{lli}}.
242 @raise SciPyError: A C{scipy} issue.
244 @raise SciPyWarning: A C{scipy} warning as exception.
245 '''
246 notOverloaded(self, callername='__call__', *llis, **wrap)
248 def _as_xyllis4(self, llis, **wrap):
249 # convert lli C{LatLon}s to tuples or C{NumPy} arrays of
250 # C{SciPy} sphericals and determine the return type
251 atype = self.numpy.array
252 wrap = _xkwds(wrap, wrap=self._wrap)
253 _as, llis = _as_llis2(llis)
254 xis, yis, _ = zip(*_xyhs(llis, **wrap)) # PYCHOK yield
255 return _as, atype(xis), atype(yis), llis
257 def _ev(self, *args): # PYCHOK no cover
258 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
259 '''
260 notOverloaded(self, *args)
262 def _eval(self, llis, **wrap): # XXX single arg, not *args
263 _as, xis, yis, _ = self._as_xyllis4(llis, **wrap)
264 try: # SciPy .ev signature: y first, then x!
265 return _as(self._ev(yis, xis))
266 except Exception as x:
267 raise _SciPyIssue(x)
269 def height(self, lats, lons, **wrap): # PYCHOK no cover
270 '''Interpolate the height for one or several lat-/longitudes.
272 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
273 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
274 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
275 and B{C{lons}} locationts (C{bool}), overriding
276 the B{C{knots}} setting.
278 @return: A single interpolated height (C{float}) or a list of
279 interpolated heights (C{float}s).
281 @raise HeightError: Insufficient or non-matching number of
282 B{C{lats}} and B{C{lons}}.
284 @raise SciPyError: A C{scipy} issue.
286 @raise SciPyWarning: A C{scipy} warning as exception.
287 '''
288 notOverloaded(self, lats, lons, **wrap)
290 def _np_sp2(self, throwarnings=False):
291 '''(INTERNAL) Import C{numpy} and C{scipy}, once.
292 '''
293 t = _HeightsBase._np_sp
294 if not t:
295 # raise SciPyWarnings, but not if
296 # scipy has been imported already
297 if throwarnings: # PYCHOK no cover
298 import sys
299 if _scipy_ not in sys.modules:
300 import warnings
301 warnings.filterwarnings(_error_)
303 sp = _xscipy(self.__class__, 1, 2)
304 np = _xnumpy(self.__class__, 1, 9)
306 _HeightsBase._np_sp = t = np, sp
307 return t
309 @Property_RO
310 def numpy(self):
311 '''Get the C{numpy} module or C{None}.
312 '''
313 np, _ = self._np_sp2()
314 return np
316 @Property_RO
317 def scipy(self):
318 '''Get the C{scipy} module or C{None}.
319 '''
320 _, sp = self._np_sp2()
321 return sp
323 @Property_RO
324 def scipy_interpolate(self):
325 '''Get the C{scipy.interpolate} module or C{None}.
326 '''
327 _ = self.scipy
328 import scipy.interpolate as spi
329 return spi
331 def _xyhs3(self, knots, **wrap):
332 # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals
333 xs, ys, hs = zip(*_xyhs(knots, name=_knots_, **wrap)) # PYCHOK yield
334 n = len(hs)
335 if n < self.kmin:
336 raise _insufficientError(self.kmin, knots=n)
337 return map1(self.numpy.array, xs, ys, hs)
340class HeightCubic(_HeightsBase):
341 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
342 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
343 C{kind='cubic'}.
344 '''
345 _interp2d = None
346 _kind = _cubic_
347 _kmin = 16
349 def __init__(self, knots, name=NN, **wrap):
350 '''New L{HeightCubic} interpolator.
352 @arg knots: The points with known height (C{LatLon}s).
353 @kwarg name: Optional name for this height interpolator (C{str}).
354 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{knots}}
355 and B{C{llis}} locations (C{bool}).
357 @raise HeightError: Insufficient number of B{C{knots}} or
358 invalid B{C{knot}}.
360 @raise ImportError: Package C{numpy} or C{scipy} not found
361 or not installed.
363 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
365 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
366 as exception.
367 '''
368 spi = self.scipy_interpolate
370 xs, ys, hs = self._xyhs3(knots, **wrap)
371 try: # SciPy.interpolate.interp2d kind 'linear' or 'cubic'
372 self._interp2d = spi.interp2d(xs, ys, hs, kind=self._kind)
373 except Exception as x:
374 raise _SciPyIssue(x)
376 if name:
377 self.name = name
379 def __call__(self, *llis, **wrap):
380 '''Interpolate the height for one or several locations.
382 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
384 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
385 as exception.
387 @see: L{_HeightsBase.__call__} for details about B{C{llis}},
388 B{C{wrap}}, return values and other exceptions.
389 '''
390 return _HeightsBase._eval(self, llis, **wrap)
392 def _ev(self, yis, xis): # PYCHOK expected
393 # to make SciPy .interp2d signature(x, y), single (x, y)
394 # match SciPy .ev signature(ys, xs), flipped multiples
395 return map(self._interp2d, xis, yis)
397 def height(self, lats, lons, **wrap):
398 '''Interpolate the height for one or several lat-/longitudes.
400 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
402 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
403 as exception.
405 @see: L{_HeightsBase.height} for details about B{C{lats}},
406 B{C{lons}}, B{C{wrap}}, return values and other exceptions.
407 '''
408 return _height_called(self, lats, lons, **wrap)
411class HeightLinear(HeightCubic):
412 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
413 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
414 C{kind='linear'}.
415 '''
416 _kind = _linear_
417 _kmin = 2
419 def __init__(self, knots, name=NN, **wrap):
420 '''New L{HeightLinear} interpolator.
422 @see: L{HeightCubic.__init__} for details about B{C{knots}},
423 B{C{name}}, B{C{wrap}} and other exceptions.
424 '''
425 HeightCubic.__init__(self, knots, name=name, **wrap)
427 if _FOR_DOCS:
428 __call__ = HeightCubic.__call__
429 height = HeightCubic.height
432class HeightLSQBiSpline(_HeightsBase):
433 '''Height interpolator using C{SciPy} U{LSQSphereBivariateSpline
434 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
435 interpolate.LSQSphereBivariateSpline.html>}.
436 '''
437 _kmin = 16 # k = 3, always
439 def __init__(self, knots, weight=None, name=NN, **wrap):
440 '''New L{HeightLSQBiSpline} interpolator.
442 @arg knots: The points with known height (C{LatLon}s).
443 @kwarg weight: Optional weight or weights for each B{C{knot}}
444 (C{scalar} or C{scalar}s).
445 @kwarg name: Optional name for this height interpolator (C{str}).
446 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{knots}}
447 and B{C{llis}} locations (C{bool}).
449 @raise HeightError: Insufficient number of B{C{knots}} or
450 an invalid B{C{knot}} or B{C{weight}}.
452 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
454 @raise ImportError: Package C{numpy} or C{scipy} not found
455 or not installed.
457 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
459 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
460 warning as exception.
461 '''
462 np = self.numpy
463 spi = self.scipy_interpolate
465 xs, ys, hs = self._xyhs3(knots, **wrap)
466 n = len(hs)
468 w = weight
469 if isscalar(w):
470 w = float(w)
471 if w <= 0:
472 raise HeightError(weight=w)
473 w = [w] * n
474 elif w is not None:
475 m, w = len2(w)
476 if m != n:
477 raise LenError(HeightLSQBiSpline, weight=m, knots=n)
478 w = map2(float, w)
479 m = min(w)
480 if m <= 0: # PYCHOK no cover
481 w = Fmt.SQUARE(weight=w.find(m))
482 raise HeightError(w, m)
483 try:
484 T = 1.0e-4 # like SciPy example
485 ps = np.array(_ordedup(xs, T, PI2 - T))
486 ts = np.array(_ordedup(ys, T, PI - T))
487 self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs,
488 ts, ps, eps=EPS, w=w).ev
489 except Exception as x:
490 raise _SciPyIssue(x)
492 if name:
493 self.name = name
495 def __call__(self, *llis, **wrap):
496 '''Interpolate the height for one or several locations.
498 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
500 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
501 warning as exception.
503 @see: L{_HeightsBase.__call__} for details about B{C{llis}},
504 B{C{wrap}}, return values and other exceptions.
505 '''
506 return _HeightsBase._eval(self, llis, **wrap)
508 def height(self, lats, lons, **wrap):
509 '''Interpolate the height for one or several lat-/longitudes.
511 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
513 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
514 warning as exception.
516 @see: L{_HeightsBase.height} for details about B{C{lats}},
517 B{C{lons}}, B{C{wrap}}, return values and other exceptions.
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=NN, **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: Optional name for this height interpolator (C{str}).
535 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{knots}}
536 and any called B{C{llis}} and height B{C{lats}}
537 and B{C{lons}} locations (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
543 or not 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, **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 if name:
562 self.name = name
564 def __call__(self, *llis, **wrap):
565 '''Interpolate the height for one or several locations.
567 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
569 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
570 warning as exception.
572 @see: L{_HeightsBase.__call__} for details about B{C{llis}},
573 B{C{wrap}}, return values and other exceptions.
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{_HeightsBase.height} for details about B{C{lats}},
586 B{C{lons}}, B{C{wrap}}, return values and other exceptions.
587 '''
588 return _height_called(self, lats, lons, **wrap)
591class _HeightIDW(_HeightBase):
592 '''(INTERNAL) Base class for U{Inverse Distance Weighting
593 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) height
594 interpolators.
596 @see: U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
597 geostatistics/Inverse-Distance-Weighting/index.html>},
598 U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
599 shepard_interp_2d/shepard_interp_2d.html>} and other C{_HeightIDW*}
600 classes.
601 '''
602 _beta = 0 # fidw inverse power
603 _func = None # formy function
604 _ks = () # knots list or tuple
605 _kwds = {} # func_ options
607 def __init__(self, knots, beta=2, name=NN, **kwds):
608 '''New C{_HeightIDW*} interpolator.
610 @arg knots: The points with known height (C{LatLon}s).
611 @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
612 @kwarg name: Optional name for this height interpolator (C{str}).
613 @kwarg kwds: Optional keyword argument for distance function,
614 retrievable with property C{kwds}.
616 @raise HeightError: Insufficient number of B{C{knots}} or
617 an invalid B{C{knot}} or B{C{beta}}.
618 '''
619 if name:
620 self.name = name
621 n, self._ks = len2(knots)
622 if n < self.kmin:
623 raise _insufficientError(self.kmin, knots=n)
624 self.beta = beta
625 self._kwds = kwds or {}
627 def __call__(self, *llis, **wrap):
628 '''Interpolate the height for one or several locations.
630 @arg llis: The location or locations (C{LatLon}, ... or
631 C{LatLon}s).
632 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
633 locations (C{bool}).
635 @return: A single interpolated height (C{float}) or a list
636 or tuple of interpolated heights (C{float}s).
638 @raise HeightError: Insufficient number of B{C{llis}}, an
639 invalid B{C{lli}} or L{pygeodesy.fidw}
640 issue.
641 '''
642 def _xy2(lls, wrap=False):
643 _w = _Wrap._latlonop(wrap)
644 try: # like _xyhs above, but degrees
645 for i, ll in enumerate(lls):
646 yield _w(ll.lon, ll.lat)
647 except AttributeError as x:
648 i = Fmt.SQUARE(llis=i)
649 raise HeightError(i, ll, cause=x)
651 _as, llis = _as_llis2(llis)
652 return _as(map(self._hIDW, *zip(*_xy2(llis, **wrap))))
654 @property_RO
655 def adjust(self):
656 '''Get the C{adjust} setting (C{bool}) or C{None}.
657 '''
658 return _xkwds_get(self._kwds, adjust=None)
660 @property
661 def beta(self):
662 '''Get the inverse distance power (C{int}).
663 '''
664 return self._beta
666 @beta.setter # PYCHOK setter!
667 def beta(self, beta):
668 '''Set the inverse distance power (C{int} 1, 2, or 3).
670 @raise HeightError: Invalid B{C{beta}}.
671 '''
672 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
674 @property_RO
675 def datum(self):
676 '''Get the C{datum} setting or the default (L{Datum}).
677 '''
678 return _xkwds_get(self._kwds, datum=self._datum)
680 def _datum_setter(self, datum):
681 '''(INTERNAL) Set the default C{datum}.
682 '''
683 d = datum or _xattr(self._ks[0], datum=None)
684 if d and d is not self._datum:
685 self._datum = _ellipsoidal_datum(d, name=self.name)
687 def _distances(self, x, y):
688 '''(INTERNAL) Yield distances to C{(x, y)}.
689 '''
690 _f, kwds = self._func, self._kwds
691 if not callable(_f): # PYCHOK no cover
692 notOverloaded(self, distance_function=_f)
693 for _, k in enumerate(self._ks):
694 yield _f(y, x, k.lat, k.lon, **kwds)
696 def _hIDW(self, x, y):
697 '''(INTERNAL) Return the IDW-interpolate height at
698 location (x, y), both C{degrees} or C{radians}.
699 '''
700 try:
701 ds = self._distances(x, y)
702 hs = (k.height for k in self._ks)
703 return fidw(hs, ds, beta=self._beta)
704 except (TypeError, ValueError) as x:
705 raise HeightError(str(x), cause=x)
707 def height(self, lats, lons, **wrap):
708 '''Interpolate the height for one or several lat-/longitudes.
710 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
711 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
712 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
713 and B{C{lons}} (C{bool}).
715 @return: A single interpolated height (C{float}) or a list of
716 interpolated heights (C{float}s).
718 @raise HeightError: Insufficient or non-matching number of
719 B{C{lats}} and B{C{lons}} or L{pygeodesy.fidw}
720 issue.
721 '''
722 return _height_called(self, lats, lons, **wrap)
724 @property_RO
725 def hypot(self):
726 '''Get the C{hypot} setting (C{callable}) or C{None}.
727 '''
728 return _xkwds_get(self._kwds, hypot=None)
730 @property_RO
731 def knots(self):
732 '''Get the B{C{knots}} (C{list} or C{tuple}).
733 '''
734 return self._ks
736 @property_RO
737 def kwds(self):
738 '''Get the optional keyword arguments (C{dict}).
739 '''
740 return self._kwds
742 @property_RO
743 def limit(self):
744 '''Get the C{limit} setting (C{degrees}) or C{None}.
745 '''
746 return _xkwds_get(self._kwds, limit=None)
748 @property_RO
749 def radius(self):
750 '''Get the C{radius} setting (C{bool}) or C{None}.
751 '''
752 return _xkwds_get(self._kwds, radius=None)
754 @property_RO
755 def scaled(self):
756 '''Get the C{scaled} setting (C{bool}) or C{None}.
757 '''
758 return _xkwds_get(self._kwds, scaled=None)
760 @property_RO
761 def wrap(self):
762 '''Get the C{wrap} setting or the default (C{bool}) or C{None}.
763 '''
764 return _xkwds_get(self._kwds, wrap=self._wrap)
767class HeightIDWcosineAndoyerLambert(_HeightIDW):
768 '''Height interpolator using U{Inverse Distance Weighting
769 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
770 and the I{angular} distance in C{radians} from function
771 L{pygeodesy.cosineAndoyerLambert_}.
772 '''
773 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
774 '''New L{HeightIDWcosineAndoyerLambert} interpolator.
776 @kwarg datum_wrap: Optional keyword arguments for function
777 L{pygeodesy.cosineAndoyerLambert}.
779 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
780 B{C{beta}}, B{C{name}} and other exceptions.
781 '''
782 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)
783 self._func = cosineAndoyerLambert
785 if _FOR_DOCS:
786 __call__ = _HeightIDW.__call__
787 height = _HeightIDW.height
790class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):
791 '''Height interpolator using U{Inverse Distance Weighting
792 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
793 and the I{angular} distance in C{radians} from function
794 L{pygeodesy.cosineForsytheAndoyerLambert_}.
795 '''
796 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
797 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
799 @kwarg datum_wrap: Optional keyword arguments for function
800 L{pygeodesy.cosineAndoyerLambert}.
802 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
803 B{C{beta}}, B{C{name}} and other exceptions.
804 '''
805 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)
806 self._func = cosineForsytheAndoyerLambert
808 if _FOR_DOCS:
809 __call__ = _HeightIDW.__call__
810 height = _HeightIDW.height
813class HeightIDWcosineLaw(_HeightIDW):
814 '''Height interpolator using U{Inverse Distance Weighting
815 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
816 I{angular} distance in C{radians} from function L{pygeodesy.cosineLaw_}.
818 @note: See note at function L{pygeodesy.vincentys_}.
819 '''
820 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
821 '''New L{HeightIDWcosineLaw} interpolator.
823 @kwarg radius_wrap: Optional keyword arguments for function
824 L{pygeodesy.cosineLaw}.
826 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
827 B{C{beta}}, B{C{name}} and other exceptions.
828 '''
829 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
830 self._func = cosineLaw
832 if _FOR_DOCS:
833 __call__ = _HeightIDW.__call__
834 height = _HeightIDW.height
837class HeightIDWdistanceTo(_HeightIDW):
838 '''Height interpolator using U{Inverse Distance Weighting
839 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
840 and the distance from the points' C{LatLon.distanceTo} method,
841 conventionally in C{meter}.
842 '''
843 def __init__(self, knots, beta=2, name=NN, **distanceTo_kwds):
844 '''New L{HeightIDWdistanceTo} interpolator.
846 @kwarg distanceTo_kwds: Optional keyword arguments for the
847 B{C{knots}}' C{LatLon.distanceTo}
848 method.
850 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
851 B{C{beta}}, B{C{name}} and other exceptions.
853 @note: All B{C{points}} I{must} be instances of the same
854 ellipsoidal or spherical C{LatLon} class, I{not
855 checked however}.
856 '''
857 _HeightIDW.__init__(self, knots, beta=beta, name=name,
858 **distanceTo_kwds)
859 ks = _distanceTo(HeightError, knots=self._ks)
860 # use knots[0] class and datum to create
861 # compatible points in _height_called
862 # instead of class LatLon_ and datum None
863 self._datum = ks[0].datum
864 self._LLis = ks[0].classof
866 def _distances(self, x, y):
867 '''(INTERNAL) Yield distances to C{(x, y)}.
868 '''
869 kwds, ll = self._kwds, self._LLis(y, x)
870 for _, k in enumerate(self._ks):
871 yield k.distanceTo(ll, **kwds)
873 if _FOR_DOCS:
874 __call__ = _HeightIDW.__call__
875 height = _HeightIDW.height
878class HeightIDWequirectangular(_HeightIDW):
879 '''Height interpolator using U{Inverse Distance Weighting
880 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
881 the I{angular} distance in C{radians squared} like function
882 L{pygeodesy.equirectangular_}.
883 '''
884 def __init__(self, knots, beta=2, name=NN, **adjust_limit_wrap): # XXX beta=1
885 '''New L{HeightIDWequirectangular} interpolator.
887 @kwarg adjust_limit_wrap: Optional keyword arguments for
888 function L{pygeodesy.equirectangular_}.
890 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
891 B{C{beta}}, B{C{name}} and other exceptions.
892 '''
893 _HeightIDW.__init__(self, knots, beta=beta, name=name,
894 **adjust_limit_wrap)
896 def _distances(self, x, y):
897 '''(INTERNAL) Yield distances to C{(x, y)}.
898 '''
899 _f, kwds = equirectangular_, self._kwds
900 for _, k in enumerate(self._ks):
901 yield _f(y, x, k.lat, k.lon, **kwds).distance2
903 if _FOR_DOCS:
904 __call__ = _HeightIDW.__call__
905 height = _HeightIDW.height
908class HeightIDWeuclidean(_HeightIDW):
909 '''Height interpolator using U{Inverse Distance Weighting
910 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
911 I{angular} distance in C{radians} from function L{pygeodesy.euclidean_}.
912 '''
913 def __init__(self, knots, beta=2, name=NN, **adjust_radius_wrap):
914 '''New L{HeightIDWeuclidean} interpolator.
916 @kwarg adjust_radius_wrap: Optional keyword arguments for
917 function L{pygeodesy.euclidean_}.
919 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
920 B{C{beta}}, B{C{name}} and other exceptions.
921 '''
922 _HeightIDW.__init__(self, knots, beta=beta, name=name,
923 **adjust_radius_wrap)
924 self._func = euclidean
926 if _FOR_DOCS:
927 __call__ = _HeightIDW.__call__
928 height = _HeightIDW.height
931class HeightIDWexact(_HeightIDW):
932 '''Height interpolator using U{Inverse Distance Weighting
933 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
934 and the I{angular} distance in C{degrees} from method
935 L{GeodesicExact}C{.Inverse}.
936 '''
937 def __init__(self, knots, beta=2, name=NN, datum=None, **wrap):
938 '''New L{HeightIDWexact} interpolator.
940 @kwarg datum: Datum to override the default C{Datums.WGS84} and
941 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
942 L{Ellipsoid2} or L{a_f2Tuple}).
943 @kwarg wrap: Optional keyword argument for method C{Inverse1}
944 of class L{geodesicx.GeodesicExact}.
946 @raise TypeError: Invalid B{C{datum}}.
948 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
949 B{C{beta}}, B{C{name}} and other exceptions.
950 '''
951 _HeightIDW.__init__(self, knots, beta=beta, name=name, **wrap)
952 self._datum_setter(datum)
953 self._func = self.datum.ellipsoid.geodesicx.Inverse1
955 if _FOR_DOCS:
956 __call__ = _HeightIDW.__call__
957 height = _HeightIDW.height
960class HeightIDWflatLocal(_HeightIDW):
961 '''Height interpolator using U{Inverse Distance Weighting
962 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
963 and the I{angular} distance in C{radians squared} like function
964 L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}.
965 '''
966 def __init__(self, knots, beta=2, name=NN, **datum_hypot_scaled_wrap):
967 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
969 @kwarg datum_hypot_scaled_wrap: Optional keyword arguments for
970 function L{pygeodesy.flatLocal_}.
972 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
973 B{C{beta}}, B{C{name}} and other exceptions.
974 '''
975 _HeightIDW.__init__(self, knots, beta=beta, name=name,
976 **datum_hypot_scaled_wrap)
977 self._func = flatLocal
979 if _FOR_DOCS:
980 __call__ = _HeightIDW.__call__
981 height = _HeightIDW.height
984class HeightIDWflatPolar(_HeightIDW):
985 '''Height interpolator using U{Inverse Distance Weighting
986 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
987 and the I{angular} distance in C{radians} from function
988 L{pygeodesy.flatPolar_}.
989 '''
990 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
991 '''New L{HeightIDWflatPolar} interpolator.
993 @kwarg radius_wrap: Optional keyword arguments for function
994 L{pygeodesy.flatPolar}.
996 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
997 B{C{beta}}, B{C{name}} and other exceptions.
998 '''
999 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1000 self._func = flatPolar
1002 if _FOR_DOCS:
1003 __call__ = _HeightIDW.__call__
1004 height = _HeightIDW.height
1007class HeightIDWhaversine(_HeightIDW):
1008 '''Height interpolator using U{Inverse Distance Weighting
1009 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1010 I{angular} distance in C{radians} from function L{pygeodesy.haversine_}.
1012 @note: See note at function L{pygeodesy.vincentys_}.
1013 '''
1014 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
1015 '''New L{HeightIDWhaversine} interpolator.
1017 @kwarg radius_wrap: Optional keyword arguments for function
1018 L{pygeodesy.haversine}.
1020 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1021 B{C{beta}}, B{C{name}} and other exceptions.
1022 '''
1023 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1024 self._func = haversine
1026 if _FOR_DOCS:
1027 __call__ = _HeightIDW.__call__
1028 height = _HeightIDW.height
1031class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny
1032 if _FOR_DOCS:
1033 __doc__ = HeightIDWflatLocal.__doc__
1034 __init__ = HeightIDWflatLocal.__init__
1035 __call__ = HeightIDWflatLocal.__call__
1036 height = HeightIDWflatLocal.height
1039class HeightIDWkarney(_HeightIDW):
1040 '''Height interpolator using U{Inverse Distance Weighting
1041 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
1042 the I{angular} distance in C{degrees} from I{Karney}'s
1043 U{geographiclib<https://PyPI.org/project/geographiclib>} U{Geodesic
1044 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}
1045 Inverse method.
1046 '''
1047 def __init__(self, knots, beta=2, name=NN, datum=None, **wrap):
1048 '''New L{HeightIDWkarney} interpolator.
1050 @kwarg datum: Datum to override the default C{Datums.WGS84} and
1051 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
1052 L{Ellipsoid2} or L{a_f2Tuple}).
1053 @kwarg wrap: Optional keyword argument for method C{Inverse1}
1054 of class L{geodesicw.Geodesic}.
1056 @raise ImportError: Package U{geographiclib
1057 <https://PyPI.org/project/geographiclib>} missing.
1059 @raise TypeError: Invalid B{C{datum}}.
1061 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1062 B{C{beta}}, B{C{name}} and other exceptions.
1063 '''
1064 _HeightIDW.__init__(self, knots, beta=beta, name=name, **wrap)
1065 self._datum_setter(datum)
1066 self._func = self.datum.ellipsoid.geodesic.Inverse1
1068 if _FOR_DOCS:
1069 __call__ = _HeightIDW.__call__
1070 height = _HeightIDW.height
1073class HeightIDWthomas(_HeightIDW):
1074 '''Height interpolator using U{Inverse Distance Weighting
1075 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1076 I{angular} distance in C{radians} from function L{pygeodesy.thomas_}.
1077 '''
1078 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
1079 '''New L{HeightIDWthomas} interpolator.
1081 @kwarg datum_wrap: Optional keyword argument for function
1082 L{pygeodesy.thomas}.
1084 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1085 B{C{beta}}, B{C{name}} and other exceptions.
1086 '''
1087 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)
1088 self._func = thomas
1090 if _FOR_DOCS:
1091 __call__ = _HeightIDW.__call__
1092 height = _HeightIDW.height
1095class HeightIDWvincentys(_HeightIDW):
1096 '''Height interpolator using U{Inverse Distance Weighting
1097 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1098 I{angular} distance in C{radians} from function L{pygeodesy.vincentys_}.
1100 @note: See note at function L{pygeodesy.vincentys_}.
1101 '''
1102 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
1103 '''New L{HeightIDWvincentys} interpolator.
1105 @kwarg radius_wrap: Optional keyword arguments for function
1106 L{pygeodesy.vincentys}.
1108 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1109 B{C{beta}}, B{C{name}} and other exceptions.
1110 '''
1111 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1112 self._func = vincentys
1114 if _FOR_DOCS:
1115 __call__ = _HeightIDW.__call__
1116 height = _HeightIDW.height
1119__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightsBase)
1121# **) MIT License
1122#
1123# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1124#
1125# Permission is hereby granted, free of charge, to any person obtaining a
1126# copy of this software and associated documentation files (the "Software"),
1127# to deal in the Software without restriction, including without limitation
1128# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1129# and/or sell copies of the Software, and to permit persons to whom the
1130# Software is furnished to do so, subject to the following conditions:
1131#
1132# The above copyright notice and this permission notice shall be included
1133# in all copies or substantial portions of the Software.
1134#
1135# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1136# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1137# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1138# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1139# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1140# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1141# OTHER DEALINGS IN THE SOFTWARE.