Coverage for pygeodesy/heights.py: 95%
312 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-11-12 13:23 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2023-11-12 13:23 -0500
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_, \
80 _STAR_, _version2
81from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS
82from pygeodesy.named import _Named, notOverloaded
83from pygeodesy.points import _distanceTo, LatLon_, Fmt, _Wrap
84from pygeodesy.props import Property_RO, property_RO
85# from pygeodesy.streprs import Fmt # from .points
86# from pygeodesy.units import Float_, Int_ # from .fmath
87# from pygeodesy.utily import _Wrap # from .points
89# from math import radians # from .formy
91__all__ = _ALL_LAZY.heights
92__version__ = '23.10.03'
94_error_ = 'error'
95_llis_ = 'llis'
96_smoothing_ = 'smoothing'
99class HeightError(PointsError):
100 '''Height interpolator C{Height...} or interpolation issue.
101 '''
102 pass
105def _alist(ais):
106 # return list of floats, not numpy.float64s
107 return list(map(float, ais))
110def _ascalar(ais): # in .geoids
111 # return single float, not numpy.float64
112 ais = list(ais) # np.array, etc. to list
113 if len(ais) != 1:
114 n = Fmt.PAREN(len=repr(ais))
115 t = _SPACE_(len(ais), _NOTEQUAL_, 1)
116 raise _AssertionError(n, txt=t)
117 return float(ais[0]) # remove np.<type>
120def _atuple(ais):
121 # return tuple of floats, not numpy.float64s
122 return tuple(map(float, ais))
125def _as_llis2(llis, m=1, Error=HeightError): # in .geoids
126 # determine return type and convert lli C{LatLon}s to list
127 if not isinstance(llis, tuple): # llis are *args
128 n = Fmt.PAREN(type_=_STAR_(NN, _llis_))
129 raise _AssertionError(n, txt=repr(llis))
131 n = len(llis)
132 if n == 1: # convert single lli to 1-item list
133 llis = llis[0]
134 try:
135 n, llis = len2(llis)
136 _as = _alist # return list of interpolated heights
137 except TypeError: # single lli
138 n, llis = 1, [llis]
139 _as = _ascalar # return single interpolated heights
140 else: # of 0, 2 or more llis
141 _as = _atuple # return tuple of interpolated heights
143 if n < m:
144 raise _insufficientError(m, Error=Error, llis=n)
145 return _as, llis
148def _height_called(inst, lats, lons, Error=HeightError, **wrap): # in .geoids
149 LLis, d = inst._LLis, inst.datum
150 if isscalar(lats) and isscalar(lons):
151 llis = LLis(lats, lons, datum=d)
152 else:
153 n, lats = len2(lats)
154 m, lons = len2(lons)
155 if n != m:
156 # format a LenError, but raise an Error
157 e = LenError(inst.__class__, lats=n, lons=m, txt=None)
158 raise e if Error is LenError else Error(str(e))
159 llis = [LLis(*t, datum=d) for t in zip(lats, lons)]
160 return inst(llis, **wrap) # __call__(lli) or __call__(llis)
163def _insufficientError(need, Error=HeightError, **name_value): # PYCHOK no cover
164 # create an insufficient Error instance
165 t = _COMMASPACE_(_insufficient_, str(need) + _PLUS_)
166 return Error(txt=t, **name_value)
169def _ordedup(ts, lo=EPS, hi=PI2-EPS):
170 # clip, order and remove duplicates
171 # p, ks = 0, []
172 # for k in sorted(max(lo, min(hi, t)) for t in ts):
173 # if k > p:
174 # ks.append(k)
175 # p = k
176 # return ks
177 return sorted(set(max(lo, min(hi, t)) for t in ts)) # list
180def _xyhs(lls, wrap=False, name=_llis_):
181 # map (lat, lon, h) to (x, y, h) in radians, offset
182 # x as 0 <= lon <= PI2 and y as 0 <= lat <= PI
183 _0, _90, _180 = _0_0, _90_0, _180_0
184 _m, _r, _w = max, radians, _Wrap._latlonop(wrap)
185 try:
186 for i, ll in enumerate(lls):
187 y, x = _w(ll.lat, ll.lon)
188 yield _m(_0, _r(x + _180)), \
189 _m(_0, _r(y + _90)), ll.height
190 except AttributeError as x:
191 i = Fmt.SQUARE(name, i)
192 raise HeightError(i, ll, cause=x)
195class _HeightBase(_Named): # in .geoids
196 '''(INTERNAL) Interpolator base class.
197 '''
198 _datum = _WGS84 # default
199 _kmin = 2 # min number of knots
200 _LLis = LatLon_ # ._height class
201 _np_sp = None # (numpy, scipy)
202 _wrap = None # wrap knots and llis
204 @property_RO
205 def datum(self):
206 '''Get the C{datum} setting or the default (L{Datum}).
207 '''
208 return self._datum
210 @property_RO
211 def kmin(self):
212 '''Get the minimum number of knots (C{int}).
213 '''
214 return self._kmin
216 @property_RO
217 def wrap(self):
218 '''Get the C{wrap} setting (C{bool}) or C{None}.
219 '''
220 return self._wrap
223class _HeightsBase(_HeightBase): # in .geoids
224 '''(INTERNAL) Interpolator base class.
225 '''
226 _np_sp = None # (numpy, scipy)
228 def __call__(self, *llis, **wrap): # PYCHOK no cover
229 '''Interpolate the height for one or several locations. I{Must be overloaded}.
231 @arg llis: One or more locations (C{LatLon}s), all positional.
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}.'''
259 notOverloaded(self, *args)
261 def _eval(self, llis, **wrap): # XXX single arg, not *args
262 _as, xis, yis, _ = self._as_xyllis4(llis, **wrap)
263 try: # SciPy .ev signature: y first, then x!
264 return _as(self._ev(yis, xis))
265 except Exception as x:
266 raise _SciPyIssue(x)
268 def height(self, lats, lons, **wrap): # PYCHOK no cover
269 '''Interpolate the height for one or several lat-/longitudes. I{Must be overloaded}.
271 @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
272 @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
273 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{lats}}
274 and B{C{lons}} locationts (C{bool}), overriding
275 the B{C{knots}} setting.
277 @return: A single interpolated height (C{float}) or a list of
278 interpolated heights (C{float}s).
280 @raise HeightError: Insufficient or non-matching number of
281 B{C{lats}} and B{C{lons}}.
283 @raise SciPyError: A C{scipy} issue.
285 @raise SciPyWarning: A C{scipy} warning as exception.
286 '''
287 notOverloaded(self, lats, lons, **wrap)
289 def _np_sp2(self, throwarnings=False):
290 '''(INTERNAL) Import C{numpy} and C{scipy}, once.
291 '''
292 t = _HeightsBase._np_sp
293 if not t:
294 # raise SciPyWarnings, but not if
295 # scipy has been imported already
296 if throwarnings: # PYCHOK no cover
297 import sys
298 if _scipy_ not in sys.modules:
299 import warnings
300 warnings.filterwarnings(_error_)
302 sp = _xscipy(self.__class__, 1, 2)
303 np = _xnumpy(self.__class__, 1, 9)
305 _HeightsBase._np_sp = t = np, sp
306 return t
308 @Property_RO
309 def numpy(self):
310 '''Get the C{numpy} module or C{None}.
311 '''
312 np, _ = self._np_sp2()
313 return np
315 @Property_RO
316 def scipy(self):
317 '''Get the C{scipy} module or C{None}.
318 '''
319 _, sp = self._np_sp2()
320 return sp
322 @Property_RO
323 def scipy_interpolate(self):
324 '''Get the C{scipy.interpolate} module or C{None}.
325 '''
326 _ = self.scipy
327 import scipy.interpolate as spi
328 return spi
330 def _scipy_version(self, **n):
331 '''Get the C{scipy} version as 2- or 3-tuple C{(major, minor, micro)}.
332 '''
333 return _version2(self.scipy.version.version, **n)
335 def _xyhs3(self, knots, **wrap):
336 # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals
337 xs, ys, hs = zip(*_xyhs(knots, name=_knots_, **wrap)) # PYCHOK yield
338 n = len(hs)
339 if n < self.kmin:
340 raise _insufficientError(self.kmin, knots=n)
341 return map1(self.numpy.array, xs, ys, hs)
344class HeightCubic(_HeightsBase):
345 '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
346 doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
347 C{kind='cubic'}.
348 '''
349 _interp2d = None
350 _kind = _cubic_
351 _kmin = 16
353 def __init__(self, knots, name=NN, **wrap):
354 '''New L{HeightCubic} interpolator.
356 @arg knots: The points with known height (C{LatLon}s).
357 @kwarg name: Optional name for this height interpolator (C{str}).
358 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{knots}}
359 and B{C{llis}} locations (C{bool}).
361 @raise HeightError: Insufficient number of B{C{knots}} or
362 invalid B{C{knot}}.
364 @raise ImportError: Package C{numpy} or C{scipy} not found
365 or not installed.
367 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
369 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
370 as exception.
371 '''
372 spi = self.scipy_interpolate
374 xs, ys, hs = self._xyhs3(knots, **wrap)
375 try: # SciPy.interpolate.interp2d kind 'linear' or 'cubic'
376 self._interp2d = spi.interp2d(xs, ys, hs, kind=self._kind)
377 except Exception as x:
378 raise _SciPyIssue(x)
380 if name:
381 self.name = name
383 def __call__(self, *llis, **wrap):
384 '''Interpolate the height for one or several locations.
386 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
388 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
389 as exception.
391 @see: L{_HeightsBase.__call__} for details about B{C{llis}},
392 B{C{wrap}}, return values and other exceptions.
393 '''
394 return _HeightsBase._eval(self, llis, **wrap)
396 def _ev(self, yis, xis): # PYCHOK expected
397 # to make SciPy .interp2d signature(x, y), single (x, y)
398 # match SciPy .ev signature(ys, xs), flipped multiples
399 return map(self._interp2d, xis, yis)
401 def height(self, lats, lons, **wrap):
402 '''Interpolate the height for one or several lat-/longitudes.
404 @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
406 @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
407 as exception.
409 @see: L{_HeightsBase.height} for details about B{C{lats}},
410 B{C{lons}}, B{C{wrap}}, return values and other exceptions.
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{HeightCubic.__init__} for details about B{C{knots}},
427 B{C{name}}, B{C{wrap}} and other exceptions.
428 '''
429 HeightCubic.__init__(self, knots, name=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=NN, **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: Optional name for this height interpolator (C{str}).
450 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{knots}}
451 and B{C{llis}} locations (C{bool}).
453 @raise HeightError: Insufficient number of B{C{knots}} or
454 an invalid B{C{knot}} or B{C{weight}}.
456 @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
458 @raise ImportError: Package C{numpy} or C{scipy} not found
459 or not installed.
461 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
463 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
464 warning as exception.
465 '''
466 np = self.numpy
467 spi = self.scipy_interpolate
469 xs, ys, hs = self._xyhs3(knots, **wrap)
470 n = len(hs)
472 w = weight
473 if isscalar(w):
474 w = float(w)
475 if w <= 0:
476 raise HeightError(weight=w)
477 w = [w] * n
478 elif w is not None:
479 m, w = len2(w)
480 if m != n:
481 raise LenError(HeightLSQBiSpline, weight=m, knots=n)
482 w = map2(float, w)
483 m = min(w)
484 if m <= 0: # PYCHOK no cover
485 w = Fmt.SQUARE(weight=w.find(m))
486 raise HeightError(w, m)
487 try:
488 T = 1.0e-4 # like SciPy example
489 ps = np.array(_ordedup(xs, T, PI2 - T))
490 ts = np.array(_ordedup(ys, T, PI - T))
491 self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs,
492 ts, ps, eps=EPS, w=w).ev
493 except Exception as x:
494 raise _SciPyIssue(x)
496 if name:
497 self.name = name
499 def __call__(self, *llis, **wrap):
500 '''Interpolate the height for one or several locations.
502 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
504 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
505 warning as exception.
507 @see: L{_HeightsBase.__call__} for details about B{C{llis}},
508 B{C{wrap}}, return values and other exceptions.
509 '''
510 return _HeightsBase._eval(self, llis, **wrap)
512 def height(self, lats, lons, **wrap):
513 '''Interpolate the height for one or several lat-/longitudes.
515 @raise SciPyError: A C{scipy} C{LSQSphereBivariateSpline} issue.
517 @raise SciPyWarning: A C{scipy} C{LSQSphereBivariateSpline}
518 warning as exception.
520 @see: L{_HeightsBase.height} for details about B{C{lats}},
521 B{C{lons}}, B{C{wrap}}, return values and other exceptions.
522 '''
523 return _height_called(self, lats, lons, **wrap)
526class HeightSmoothBiSpline(_HeightsBase):
527 '''Height interpolator using C{SciPy} U{SmoothSphereBivariateSpline
528 <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
529 interpolate.SmoothSphereBivariateSpline.html>}.
530 '''
531 _kmin = 16 # k = 3, always
533 def __init__(self, knots, s=4, name=NN, **wrap):
534 '''New L{HeightSmoothBiSpline} interpolator.
536 @arg knots: The points with known height (C{LatLon}s).
537 @kwarg s: The spline smoothing factor (C{scalar}), default C{4}.
538 @kwarg name: Optional name for this height interpolator (C{str}).
539 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{knots}}
540 and any called B{C{llis}} and height B{C{lats}}
541 and B{C{lons}} locations (C{bool}).
543 @raise HeightError: Insufficient number of B{C{knots}} or
544 an invalid B{C{knot}} or B{C{s}}.
546 @raise ImportError: Package C{numpy} or C{scipy} not found
547 or not installed.
549 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
551 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
552 warning as exception.
553 '''
554 spi = self.scipy_interpolate
556 s = Float_(s, name=_smoothing_, Error=HeightError, low=4)
558 xs, ys, hs = self._xyhs3(knots, **wrap)
559 try:
560 self._ev = spi.SmoothSphereBivariateSpline(ys, xs, hs,
561 eps=EPS, s=s).ev
562 except Exception as x:
563 raise _SciPyIssue(x)
565 if name:
566 self.name = name
568 def __call__(self, *llis, **wrap):
569 '''Interpolate the height for one or several locations.
571 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
573 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
574 warning as exception.
576 @see: L{_HeightsBase.__call__} for details about B{C{llis}},
577 B{C{wrap}}, return values and other exceptions.
578 '''
579 return _HeightsBase._eval(self, llis, **wrap)
581 def height(self, lats, lons, **wrap):
582 '''Interpolate the height for one or several lat-/longitudes.
584 @raise SciPyError: A C{scipy} C{SmoothSphereBivariateSpline} issue.
586 @raise SciPyWarning: A C{scipy} C{SmoothSphereBivariateSpline}
587 warning as exception.
589 @see: L{_HeightsBase.height} for details about B{C{lats}},
590 B{C{lons}}, B{C{wrap}}, return values and other exceptions.
591 '''
592 return _height_called(self, lats, lons, **wrap)
595class _HeightIDW(_HeightBase):
596 '''(INTERNAL) Base class for U{Inverse Distance Weighting
597 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) height
598 interpolators.
600 @see: U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
601 geostatistics/Inverse-Distance-Weighting/index.html>},
602 U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
603 shepard_interp_2d/shepard_interp_2d.html>} and other C{_HeightIDW*}
604 classes.
605 '''
606 _beta = 0 # fidw inverse power
607 _func = None # formy function
608 _ks = () # knots list or tuple
609 _kwds = {} # func_ options
611 def __init__(self, knots, beta=2, name=NN, **kwds):
612 '''New C{_HeightIDW*} interpolator.
614 @arg knots: The points with known height (C{LatLon}s).
615 @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
616 @kwarg name: Optional name for this height interpolator (C{str}).
617 @kwarg kwds: Optional keyword argument for distance function,
618 retrievable with property C{kwds}.
620 @raise HeightError: Insufficient number of B{C{knots}} or
621 an invalid B{C{knot}} or B{C{beta}}.
622 '''
623 if name:
624 self.name = name
625 n, self._ks = len2(knots)
626 if n < self.kmin:
627 raise _insufficientError(self.kmin, knots=n)
628 self.beta = beta
629 self._kwds = kwds or {}
631 def __call__(self, *llis, **wrap):
632 '''Interpolate the height for one or several locations.
634 @arg llis: One or more locations (C{LatLon}s), all positional.
635 @kwarg wrap: If C{True}, wrap or I{normalize} all B{C{llis}}
636 locations (C{bool}).
638 @return: A single interpolated height (C{float}) or a list
639 or tuple of interpolated heights (C{float}s).
641 @raise HeightError: Insufficient number of B{C{llis}}, an
642 invalid B{C{lli}} or L{pygeodesy.fidw}
643 issue.
644 '''
645 def _xy2(lls, wrap=False):
646 _w = _Wrap._latlonop(wrap)
647 try: # like _xyhs above, but degrees
648 for i, ll in enumerate(lls):
649 yield _w(ll.lon, ll.lat)
650 except AttributeError as x:
651 i = Fmt.SQUARE(llis=i)
652 raise HeightError(i, ll, cause=x)
654 _as, llis = _as_llis2(llis)
655 return _as(map(self._hIDW, *zip(*_xy2(llis, **wrap))))
657 @property_RO
658 def adjust(self):
659 '''Get the C{adjust} setting (C{bool}) or C{None}.
660 '''
661 return _xkwds_get(self._kwds, adjust=None)
663 @property
664 def beta(self):
665 '''Get the inverse distance power (C{int}).
666 '''
667 return self._beta
669 @beta.setter # PYCHOK setter!
670 def beta(self, beta):
671 '''Set the inverse distance power (C{int} 1, 2, or 3).
673 @raise HeightError: Invalid B{C{beta}}.
674 '''
675 self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
677 @property_RO
678 def datum(self):
679 '''Get the C{datum} setting or the default (L{Datum}).
680 '''
681 return _xkwds_get(self._kwds, datum=self._datum)
683 def _datum_setter(self, datum):
684 '''(INTERNAL) Set the default C{datum}.
685 '''
686 d = datum or _xattr(self._ks[0], datum=None)
687 if d and d is not self._datum:
688 self._datum = _ellipsoidal_datum(d, name=self.name)
690 def _distances(self, x, y):
691 '''(INTERNAL) Yield distances to C{(x, y)}.
692 '''
693 _f, kwds = self._func, self._kwds
694 if not callable(_f): # PYCHOK no cover
695 notOverloaded(self, distance_function=_f)
696 for _, k in enumerate(self._ks):
697 yield _f(y, x, k.lat, k.lon, **kwds)
699 def _hIDW(self, x, y):
700 '''(INTERNAL) Return the IDW-interpolate height at
701 location (x, y), both C{degrees} or C{radians}.
702 '''
703 try:
704 ds = self._distances(x, y)
705 hs = (k.height for k in self._ks)
706 return fidw(hs, ds, beta=self._beta)
707 except (TypeError, ValueError) as x:
708 raise HeightError(str(x), cause=x)
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
722 B{C{lats}} and B{C{lons}} or L{pygeodesy.fidw}
723 issue.
724 '''
725 return _height_called(self, lats, lons, **wrap)
727 @property_RO
728 def hypot(self):
729 '''Get the C{hypot} setting (C{callable}) or C{None}.
730 '''
731 return _xkwds_get(self._kwds, hypot=None)
733 @property_RO
734 def knots(self):
735 '''Get the B{C{knots}} (C{list} or C{tuple}).
736 '''
737 return self._ks
739 @property_RO
740 def kwds(self):
741 '''Get the optional keyword arguments (C{dict}).
742 '''
743 return self._kwds
745 @property_RO
746 def limit(self):
747 '''Get the C{limit} setting (C{degrees}) or C{None}.
748 '''
749 return _xkwds_get(self._kwds, limit=None)
751 @property_RO
752 def radius(self):
753 '''Get the C{radius} setting (C{bool}) or C{None}.
754 '''
755 return _xkwds_get(self._kwds, radius=None)
757 @property_RO
758 def scaled(self):
759 '''Get the C{scaled} setting (C{bool}) or C{None}.
760 '''
761 return _xkwds_get(self._kwds, scaled=None)
763 @property_RO
764 def wrap(self):
765 '''Get the C{wrap} setting or the default (C{bool}) or C{None}.
766 '''
767 return _xkwds_get(self._kwds, wrap=self._wrap)
770class HeightIDWcosineAndoyerLambert(_HeightIDW):
771 '''Height interpolator using U{Inverse Distance Weighting
772 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
773 and the I{angular} distance in C{radians} from function
774 L{pygeodesy.cosineAndoyerLambert_}.
775 '''
776 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
777 '''New L{HeightIDWcosineAndoyerLambert} interpolator.
779 @kwarg datum_wrap: Optional keyword arguments for function
780 L{pygeodesy.cosineAndoyerLambert}.
782 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
783 B{C{beta}}, B{C{name}} and other exceptions.
784 '''
785 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)
786 self._func = cosineAndoyerLambert
788 if _FOR_DOCS:
789 __call__ = _HeightIDW.__call__
790 height = _HeightIDW.height
793class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):
794 '''Height interpolator using U{Inverse Distance Weighting
795 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
796 and the I{angular} distance in C{radians} from function
797 L{pygeodesy.cosineForsytheAndoyerLambert_}.
798 '''
799 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
800 '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
802 @kwarg datum_wrap: Optional keyword arguments for function
803 L{pygeodesy.cosineAndoyerLambert}.
805 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
806 B{C{beta}}, B{C{name}} and other exceptions.
807 '''
808 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)
809 self._func = cosineForsytheAndoyerLambert
811 if _FOR_DOCS:
812 __call__ = _HeightIDW.__call__
813 height = _HeightIDW.height
816class HeightIDWcosineLaw(_HeightIDW):
817 '''Height interpolator using U{Inverse Distance Weighting
818 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
819 I{angular} distance in C{radians} from function L{pygeodesy.cosineLaw_}.
821 @note: See note at function L{pygeodesy.vincentys_}.
822 '''
823 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
824 '''New L{HeightIDWcosineLaw} interpolator.
826 @kwarg radius_wrap: Optional keyword arguments for function
827 L{pygeodesy.cosineLaw}.
829 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
830 B{C{beta}}, B{C{name}} and other exceptions.
831 '''
832 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
833 self._func = cosineLaw
835 if _FOR_DOCS:
836 __call__ = _HeightIDW.__call__
837 height = _HeightIDW.height
840class HeightIDWdistanceTo(_HeightIDW):
841 '''Height interpolator using U{Inverse Distance Weighting
842 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
843 and the distance from the points' C{LatLon.distanceTo} method,
844 conventionally in C{meter}.
845 '''
846 def __init__(self, knots, beta=2, name=NN, **distanceTo_kwds):
847 '''New L{HeightIDWdistanceTo} interpolator.
849 @kwarg distanceTo_kwds: Optional keyword arguments for the
850 B{C{knots}}' C{LatLon.distanceTo}
851 method.
853 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
854 B{C{beta}}, B{C{name}} and other exceptions.
856 @note: All B{C{points}} I{must} be instances of the same
857 ellipsoidal or spherical C{LatLon} class, I{not
858 checked however}.
859 '''
860 _HeightIDW.__init__(self, knots, beta=beta, name=name,
861 **distanceTo_kwds)
862 ks = _distanceTo(HeightError, knots=self._ks)
863 # use knots[0] class and datum to create
864 # compatible points in _height_called
865 # instead of class LatLon_ and datum None
866 self._datum = ks[0].datum
867 self._LLis = ks[0].classof
869 def _distances(self, x, y):
870 '''(INTERNAL) Yield distances to C{(x, y)}.
871 '''
872 kwds, ll = self._kwds, self._LLis(y, x)
873 for _, k in enumerate(self._ks):
874 yield k.distanceTo(ll, **kwds)
876 if _FOR_DOCS:
877 __call__ = _HeightIDW.__call__
878 height = _HeightIDW.height
881class HeightIDWequirectangular(_HeightIDW):
882 '''Height interpolator using U{Inverse Distance Weighting
883 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
884 the I{angular} distance in C{radians squared} like function
885 L{pygeodesy.equirectangular_}.
886 '''
887 def __init__(self, knots, beta=2, name=NN, **adjust_limit_wrap): # XXX beta=1
888 '''New L{HeightIDWequirectangular} interpolator.
890 @kwarg adjust_limit_wrap: Optional keyword arguments for
891 function L{pygeodesy.equirectangular_}.
893 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
894 B{C{beta}}, B{C{name}} and other exceptions.
895 '''
896 _HeightIDW.__init__(self, knots, beta=beta, name=name,
897 **adjust_limit_wrap)
899 def _distances(self, x, y):
900 '''(INTERNAL) Yield distances to C{(x, y)}.
901 '''
902 _f, kwds = equirectangular_, self._kwds
903 for _, k in enumerate(self._ks):
904 yield _f(y, x, k.lat, k.lon, **kwds).distance2
906 if _FOR_DOCS:
907 __call__ = _HeightIDW.__call__
908 height = _HeightIDW.height
911class HeightIDWeuclidean(_HeightIDW):
912 '''Height interpolator using U{Inverse Distance Weighting
913 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
914 I{angular} distance in C{radians} from function L{pygeodesy.euclidean_}.
915 '''
916 def __init__(self, knots, beta=2, name=NN, **adjust_radius_wrap):
917 '''New L{HeightIDWeuclidean} interpolator.
919 @kwarg adjust_radius_wrap: Optional keyword arguments for
920 function L{pygeodesy.euclidean_}.
922 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
923 B{C{beta}}, B{C{name}} and other exceptions.
924 '''
925 _HeightIDW.__init__(self, knots, beta=beta, name=name,
926 **adjust_radius_wrap)
927 self._func = euclidean
929 if _FOR_DOCS:
930 __call__ = _HeightIDW.__call__
931 height = _HeightIDW.height
934class HeightIDWexact(_HeightIDW):
935 '''Height interpolator using U{Inverse Distance Weighting
936 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
937 and the I{angular} distance in C{degrees} from method
938 L{GeodesicExact}C{.Inverse}.
939 '''
940 def __init__(self, knots, beta=2, name=NN, datum=None, **wrap):
941 '''New L{HeightIDWexact} interpolator.
943 @kwarg datum: Datum to override the default C{Datums.WGS84} and
944 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
945 L{Ellipsoid2} or L{a_f2Tuple}).
946 @kwarg wrap: Optional keyword argument for method C{Inverse1}
947 of class L{geodesicx.GeodesicExact}.
949 @raise TypeError: Invalid B{C{datum}}.
951 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
952 B{C{beta}}, B{C{name}} and other exceptions.
953 '''
954 _HeightIDW.__init__(self, knots, beta=beta, name=name, **wrap)
955 self._datum_setter(datum)
956 self._func = self.datum.ellipsoid.geodesicx.Inverse1
958 if _FOR_DOCS:
959 __call__ = _HeightIDW.__call__
960 height = _HeightIDW.height
963class HeightIDWflatLocal(_HeightIDW):
964 '''Height interpolator using U{Inverse Distance Weighting
965 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
966 and the I{angular} distance in C{radians squared} like function
967 L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}.
968 '''
969 def __init__(self, knots, beta=2, name=NN, **datum_hypot_scaled_wrap):
970 '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
972 @kwarg datum_hypot_scaled_wrap: Optional keyword arguments for
973 function L{pygeodesy.flatLocal_}.
975 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
976 B{C{beta}}, B{C{name}} and other exceptions.
977 '''
978 _HeightIDW.__init__(self, knots, beta=beta, name=name,
979 **datum_hypot_scaled_wrap)
980 self._func = flatLocal
982 if _FOR_DOCS:
983 __call__ = _HeightIDW.__call__
984 height = _HeightIDW.height
987class HeightIDWflatPolar(_HeightIDW):
988 '''Height interpolator using U{Inverse Distance Weighting
989 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
990 and the I{angular} distance in C{radians} from function
991 L{pygeodesy.flatPolar_}.
992 '''
993 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
994 '''New L{HeightIDWflatPolar} interpolator.
996 @kwarg radius_wrap: Optional keyword arguments for function
997 L{pygeodesy.flatPolar}.
999 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1000 B{C{beta}}, B{C{name}} and other exceptions.
1001 '''
1002 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1003 self._func = flatPolar
1005 if _FOR_DOCS:
1006 __call__ = _HeightIDW.__call__
1007 height = _HeightIDW.height
1010class HeightIDWhaversine(_HeightIDW):
1011 '''Height interpolator using U{Inverse Distance Weighting
1012 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1013 I{angular} distance in C{radians} from function L{pygeodesy.haversine_}.
1015 @note: See note at function L{pygeodesy.vincentys_}.
1016 '''
1017 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
1018 '''New L{HeightIDWhaversine} interpolator.
1020 @kwarg radius_wrap: Optional keyword arguments for function
1021 L{pygeodesy.haversine}.
1023 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1024 B{C{beta}}, B{C{name}} and other exceptions.
1025 '''
1026 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1027 self._func = haversine
1029 if _FOR_DOCS:
1030 __call__ = _HeightIDW.__call__
1031 height = _HeightIDW.height
1034class HeightIDWhubeny(HeightIDWflatLocal): # for Karl Hubeny
1035 if _FOR_DOCS:
1036 __doc__ = HeightIDWflatLocal.__doc__
1037 __init__ = HeightIDWflatLocal.__init__
1038 __call__ = HeightIDWflatLocal.__call__
1039 height = HeightIDWflatLocal.height
1042class HeightIDWkarney(_HeightIDW):
1043 '''Height interpolator using U{Inverse Distance Weighting
1044 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
1045 the I{angular} distance in C{degrees} from I{Karney}'s
1046 U{geographiclib<https://PyPI.org/project/geographiclib>} U{Geodesic
1047 <https://GeographicLib.SourceForge.io/Python/doc/code.html>}
1048 Inverse method.
1049 '''
1050 def __init__(self, knots, beta=2, name=NN, datum=None, **wrap):
1051 '''New L{HeightIDWkarney} interpolator.
1053 @kwarg datum: Datum to override the default C{Datums.WGS84} and
1054 first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
1055 L{Ellipsoid2} or L{a_f2Tuple}).
1056 @kwarg wrap: Optional keyword argument for method C{Inverse1}
1057 of class L{geodesicw.Geodesic}.
1059 @raise ImportError: Package U{geographiclib
1060 <https://PyPI.org/project/geographiclib>} missing.
1062 @raise TypeError: Invalid B{C{datum}}.
1064 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1065 B{C{beta}}, B{C{name}} and other exceptions.
1066 '''
1067 _HeightIDW.__init__(self, knots, beta=beta, name=name, **wrap)
1068 self._datum_setter(datum)
1069 self._func = self.datum.ellipsoid.geodesic.Inverse1
1071 if _FOR_DOCS:
1072 __call__ = _HeightIDW.__call__
1073 height = _HeightIDW.height
1076class HeightIDWthomas(_HeightIDW):
1077 '''Height interpolator using U{Inverse Distance Weighting
1078 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1079 I{angular} distance in C{radians} from function L{pygeodesy.thomas_}.
1080 '''
1081 def __init__(self, knots, beta=2, name=NN, **datum_wrap):
1082 '''New L{HeightIDWthomas} interpolator.
1084 @kwarg datum_wrap: Optional keyword argument for function
1085 L{pygeodesy.thomas}.
1087 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1088 B{C{beta}}, B{C{name}} and other exceptions.
1089 '''
1090 _HeightIDW.__init__(self, knots, beta=beta, name=name, **datum_wrap)
1091 self._func = thomas
1093 if _FOR_DOCS:
1094 __call__ = _HeightIDW.__call__
1095 height = _HeightIDW.height
1098class HeightIDWvincentys(_HeightIDW):
1099 '''Height interpolator using U{Inverse Distance Weighting
1100 <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and the
1101 I{angular} distance in C{radians} from function L{pygeodesy.vincentys_}.
1103 @note: See note at function L{pygeodesy.vincentys_}.
1104 '''
1105 def __init__(self, knots, beta=2, name=NN, **radius_wrap):
1106 '''New L{HeightIDWvincentys} interpolator.
1108 @kwarg radius_wrap: Optional keyword arguments for function
1109 L{pygeodesy.vincentys}.
1111 @see: L{_HeightIDW.__init__} for details about B{C{knots}},
1112 B{C{beta}}, B{C{name}} and other exceptions.
1113 '''
1114 _HeightIDW.__init__(self, knots, beta=beta, name=name, **radius_wrap)
1115 self._func = vincentys
1117 if _FOR_DOCS:
1118 __call__ = _HeightIDW.__call__
1119 height = _HeightIDW.height
1122__all__ += _ALL_DOCS(_HeightBase, _HeightIDW, _HeightsBase)
1124# **) MIT License
1125#
1126# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1127#
1128# Permission is hereby granted, free of charge, to any person obtaining a
1129# copy of this software and associated documentation files (the "Software"),
1130# to deal in the Software without restriction, including without limitation
1131# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1132# and/or sell copies of the Software, and to permit persons to whom the
1133# Software is furnished to do so, subject to the following conditions:
1134#
1135# The above copyright notice and this permission notice shall be included
1136# in all copies or substantial portions of the Software.
1137#
1138# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1139# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1140# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1141# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1142# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1143# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1144# OTHER DEALINGS IN THE SOFTWARE.