1
2# -*- coding: utf-8 -*-
3
4u'''Height interpolation methods.
5
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 C{known} heights.
14
15Classes L{HeightCubic} and L{HeightLinear} require package U{numpy
16<https://PyPI.org/project/numpy>}, classes L{HeightLSQBiSpline} and
17L{HeightSmoothBiSpline} require package U{scipy<https://SciPy.org>}.
18Classes L{HeightIDWkarney} and L{HeightIDWdistanceTo} -iff used with
19L{ellipsoidalKarney.LatLon} points- require I{Karney}'s U{geographiclib
20<https://PyPI.org/project/geographiclib>} to be installed.
21
22B{Typical usage} is as follows.  First, create an interpolator from a
23given set of C{LatLon} points with C{known} heights, called C{knots}.
24
25C{>>> hinterpolator = HeightXyz(knots, **options)}
26
27Then, get the interpolated height of other C{LatLon} location(s) with
28
29C{>>> h = hinterpolator(ll)}
30
31or
32
33C{>>> h0, h1, h2, ... = hinterpolator(ll0, ll1, ll2, ...)}
34
35or
36
37C{>>> hs = hinterpolator(lls)  # list, tuple, generator, ...}
38
39For separate lat- and longitudes invoke the C{.height} method
40
41C{>>> h = hinterpolator.height(lat, lon)}
42
43or
44
45C{>>> h0, h1, h2, ... = hinterpolator.height(lats, lons)  # lists, tuples, ...}
46
47
48The C{knots} do not need to be ordered for any of the height
49interpolators.
50
51Errors from C{scipy} as raised as L{SciPyError}s.  Warnings issued
52by C{scipy} can be thrown as L{SciPyWarning} exceptions, provided
53Python C{warnings} are filtered accordingly, see L{SciPyWarning}.
54
55@see: U{SciPy<https://docs.SciPy.org/doc/scipy/reference/interpolate.html>}.
56'''
57# make sure int/int division yields float quotient, see .basics
58from __future__ import division
59
60from pygeodesy.basics import isscalar, len2, map1, map2, _xnumpy, _xscipy
61from pygeodesy.datums import _ellipsoidal_datum, _WGS84
62from pygeodesy.errors import _AssertionError, LenError, PointsError, _SciPyIssue
63from pygeodesy.fmath import fidw, hypot2
64from pygeodesy.formy import cosineAndoyerLambert_, cosineForsytheAndoyerLambert_, \
65                            cosineLaw_, euclidean_, flatPolar_, haversine_, \
66                            _scale_rad, thomas_, vincentys_
67from pygeodesy.interns import EPS, NN, PI, PI2, PI_2, _cubic_, _datum_, \
68                             _distanceTo_, _knots_, _len_, _linear_, _scipy_, \
69                             _0_0
70from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _FOR_DOCS
71from pygeodesy.named import _Named, notOverloaded
72from pygeodesy.points import LatLon_
73from pygeodesy.props import Property_RO
74from pygeodesy.streprs import _boolkwds, Fmt
75from pygeodesy.units import Int_
76from pygeodesy.utily import radiansPI, radiansPI2, unrollPI
77
78__all__ = _ALL_LAZY.heights
79__version__ = '21.06.01'
80
81
82class HeightError(PointsError):
83    '''Height interpolator C{Height...} or interpolation issue.
84    '''
85    pass
86
87
88def _alist(ais):
89    # return list of floats, not numpy.float64s
90    return list(map(float, ais))
91
92
93def _allis2(llis, m=1, Error=HeightError):  # imported by .geoids
94    # dtermine return type and convert lli C{LatLon}s to list
95    if not isinstance(llis, tuple):  # llis are *args
96        raise _AssertionError('type(%s): %r' % ('*llis', llis))
97
98    n = len(llis)
99    if n == 1:  # convert single lli to 1-item list
100        llis = llis[0]
101        try:
102            n, llis = len2(llis)
103            _as = _alist  # return list of interpolated heights
104        except TypeError:  # single lli
105            n, llis = 1, [llis]
106            _as = _ascalar  # return single interpolated heights
107    else:  # of 0, 2 or more llis
108        _as = _atuple  # return tuple of interpolated heights
109
110    if n < m:
111        raise _insufficientError(m, Error=Error, llis=n)
112    return _as, llis
113
114
115def _ascalar(ais):  # imported by .geoids
116    # return single float, not numpy.float64
117    ais = list(ais)  # np.array, etc. to list
118    if len(ais) != 1:
119        raise _AssertionError('%s(%r): %s != 1' % (_len_, ais, len(ais)))
120    return float(ais[0])  # remove np.<type>
121
122
123def _atuple(ais):
124    # return tuple of floats, not numpy.float64s
125    return tuple(map(float, ais))
126
127
128def _axyllis4(atype, llis, m=1, off=True):
129    # convert lli C{LatLon}s to tuples or C{NumPy} arrays of
130    # C{SciPy} sphericals and determine the return type
131    _as, llis = _allis2(llis, m=m)
132    xis, yis, _ =  zip(*_xyhs(llis, off=off))  # PYCHOK unzip
133    return _as, atype(xis), atype(yis), llis
134
135
136def _insufficientError(need, Error=HeightError, **name_value):
137    # create an insufficient Error instance
138    t = 'insufficient, need %s' % (need,)
139    return Error(txt=t, **name_value)
140
141
142def _ordedup(ts, lo=EPS, hi=PI2-EPS):
143    # clip, order and remove duplicates
144    p, ks = 0, []
145    for k in sorted(max(lo, min(hi, t)) for t in ts):
146        if k > p:
147            ks.append(k)
148            p = k
149    return ks
150
151
152def _xyhs(lls, off=True, name='llis'):
153    # map (lat, lon, h) to (x, y, h) in radians, offset as
154    # x: 0 <= lon <= PI2, y: 0 <= lat <= PI if off is True
155    # else x: -PI <= lon <= PI, y: -PI_2 <= lat <= PI_2
156    if off:
157        xf = yf = _0_0
158    else:  # undo offset
159        xf, yf = PI, PI_2
160    try:
161        for i, ll in enumerate(lls):
162            yield (max(0.0, radiansPI2(ll.lon + 180.0)) - xf), \
163                  (max(0.0, radiansPI( ll.lat +  90.0)) - yf), ll.height
164    except AttributeError as x:
165        i = Fmt.SQUARE(name, i)
166        raise HeightError(i, ll, txt=str(x))
167
168
169def _xyhs3(atype, m, knots, off=True):
170    # convert knot C{LatLon}s to tuples or C{NumPy} arrays and C{SciPy} sphericals
171    xs, ys, hs = zip(*_xyhs(knots, off=off, name=_knots_))  # PYCHOK unzip
172    n = len(hs)
173    if n < m:
174        raise _insufficientError(m, knots=n)
175    return map1(atype, xs, ys, hs)
176
177
178class _HeightBase(_Named):  # imported by .geoids
179    '''(INTERNAL) Interpolator base class.
180    '''
181    _adjust = None     # not applicable
182    _datum  = None     # ._height datum
183    _kmin   = 2        # min number of knots
184    _LLis   = LatLon_  # ._height class
185    _np     = None     # numpy
186    _np_v   = None     # version
187    _spi    = None     # scipy.interpolate
188    _sp_v   = None     # version
189    _wrap   = None     # not applicable
190
191    def __call__(self, *args):  # PYCHOK no cover
192        '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
193        '''
194        notOverloaded(self, callername='__call__', *args)
195
196    @Property_RO
197    def adjust(self):
198        '''Get the adjust setting (C{bool} or C{None} if not applicable).
199        '''
200        return self._adjust
201
202    def _axyllis4(self, llis):
203        return _axyllis4(self._np.array, llis)
204
205    @Property_RO
206    def datum(self):
207        '''Get the datum (L{Datum} or C{None} if not applicable).
208        '''
209        return self._datum
210
211    def _ev(self, *args):  # PYCHOK no cover
212        '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
213        '''
214        notOverloaded(self, *args)
215
216    def _eval(self, llis):  # XXX single arg, not *args
217        _as, xis, yis, _ = self._axyllis4(llis)
218        try:  # SciPy .ev signature: y first, then x!
219            return _as(self._ev(yis, xis))
220        except Exception as x:
221            raise _SciPyIssue(x)
222
223    def _height(self, lats, lons, Error=HeightError):
224        LLis, d = self._LLis, self.datum
225        if isscalar(lats) and isscalar(lons):
226            llis = LLis(lats, lons, datum=d)
227        else:
228            n, lats = len2(lats)
229            m, lons = len2(lons)
230            if n != m:
231                # format a LenError, but raise an Error
232                e = LenError(self.__class__, lats=n, lons=m, txt=None)
233                raise e if Error is LenError else Error(str(e))
234            llis = [LLis(*ll, datum=d) for ll in zip(lats, lons)]
235        return self(llis)  # __call__(lli) or __call__(llis)
236
237    @Property_RO
238    def kmin(self):
239        '''Get the minimum number of knots (C{int}).
240        '''
241        return self._kmin
242
243    def _NumSciPy(self, throwarnings=False):
244        '''(INTERNAL) Import C{numpy} and C{scipy}, once.
245        '''
246        if throwarnings:  # raise SciPyWarnings, but ...
247            # ... not if scipy has been imported already
248            import sys
249            if _scipy_ not in sys.modules:
250                import warnings
251                warnings.filterwarnings('error')
252
253        np  = _HeightBase._np
254        spi = _HeightBase._spi
255        if None in (np, spi):
256
257            sp = _xscipy(self.__class__, 1, 2)
258            import scipy.interpolate as spi
259            np = _xnumpy(self.__class__, 1, 9)
260
261            _HeightBase._np   = np
262            _HeightBase._np_v = np.__version__
263            _HeightBase._spi  = spi
264            _HeightBase._sp_v = sp.__version__
265
266        return np, spi  # XXX spi not sp!
267
268    def _xyhs3(self, knots):
269        return _xyhs3(self._np.array, self._kmin, knots)
270
271    @Property_RO
272    def wrap(self):
273        '''Get the wrap setting (C{bool} or C{None} if not applicable).
274        '''
275        return self._wrap
276
277
278class HeightCubic(_HeightBase):
279    '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
280       doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
281       C{kind='cubic'}.
282    '''
283    _interp2d =  None
284    _kind     = _cubic_
285    _kmin     =  16
286
287    def __init__(self, knots, name=NN):
288        '''New L{HeightCubic} interpolator.
289
290           @arg knots: The points with known height (C{LatLon}s).
291           @kwarg name: Optional name for this height interpolator (C{str}).
292
293           @raise HeightError: Insufficient number of B{C{knots}} or
294                               invalid B{C{knot}}.
295
296           @raise ImportError: Package C{numpy} or C{scipy} not found
297                               or not installed.
298
299           @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
300
301           @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
302                                as exception.
303        '''
304        _, spi = self._NumSciPy()
305
306        xs, ys, hs = self._xyhs3(knots)
307        try:  # SciPy.interpolate.interp2d kind 'linear' or 'cubic'
308            self._interp2d = spi.interp2d(xs, ys, hs, kind=self._kind)
309        except Exception as x:
310            raise _SciPyIssue(x)
311
312        if name:
313            self.name = name
314
315    def __call__(self, *llis):
316        '''Interpolate the height for one or several locations.
317
318           @arg llis: The location or locations (C{LatLon}, ... or
319                      C{LatLon}s).
320
321           @return: A single interpolated height (C{float}) or a list
322                    or tuple of interpolated heights (C{float}s).
323
324           @raise HeightError: Insufficient number of B{C{llis}} or
325                               an invalid B{C{lli}}.
326
327           @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
328
329           @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
330                                as exception.
331        '''
332        return _HeightBase._eval(self, llis)
333
334    def _ev(self, yis, xis):  # PYCHOK expected
335        # to make SciPy .interp2d signature(x, y), single (x, y)
336        # match SciPy .ev signature(ys, xs), flipped multiples
337        return map(self._interp2d, xis, yis)
338
339    def height(self, lats, lons):
340        '''Interpolate the height for one or several lat-/longitudes.
341
342           @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
343           @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
344
345           @return: A single interpolated height (C{float}) or a list of
346                    interpolated heights (C{float}s).
347
348           @raise HeightError: Insufficient or non-matching number of
349                               B{C{lats}} and B{C{lons}}.
350
351           @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
352
353           @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
354                                as exception.
355        '''
356        return _HeightBase._height(self, lats, lons)
357
358
359class HeightLinear(HeightCubic):
360    '''Height interpolator based on C{SciPy} U{interp2d<https://docs.SciPy.org/
361       doc/scipy/reference/generated/scipy.interpolate.interp2d.html>}
362       C{kind='linear'}.
363    '''
364    _kind = _linear_
365    _kmin =  2
366
367    def __init__(self, knots, name=NN):
368        '''New L{HeightLinear} interpolator.
369
370           @arg knots: The points with known height (C{LatLon}s).
371           @kwarg name: Optional name for this height interpolator (C{str}).
372
373           @raise HeightError: Insufficient number of B{C{knots}} or
374                               an invalid B{C{knot}}.
375
376           @raise ImportError: Package C{numpy} or C{scipy} not found
377                               or not installed.
378
379           @raise SciPyError: A C{scipy.interpolate.interp2d} issue.
380
381           @raise SciPyWarning: A C{scipy.interpolate.interp2d} warning
382                                as exception.
383        '''
384        HeightCubic.__init__(self, knots, name=name)
385
386    if _FOR_DOCS:
387        __call__ = HeightCubic.__call__
388        height   = HeightCubic.height
389
390
391class _HeightIDW(_HeightBase):
392    '''(INTERNAL) Base class for U{Inverse Distance Weighting
393       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
394       height interpolators.
395    '''
396    _beta = 0     # inverse distance power
397    _hs   = ()    # known heights
398    _xs   = ()    # knot lons
399    _ys   = ()    # knot lats
400
401    def __init__(self, knots, beta=2, name=NN, **wrap_adjust):
402        '''New L{_HeightIDW} interpolator.
403        '''
404        self._xs, self._ys, self._hs = _xyhs3(tuple, self._kmin, knots, off=False)
405        self.beta = beta
406        if name:
407            self.name = name
408        if wrap_adjust:
409            _boolkwds(self, **wrap_adjust)
410
411    def __call__(self, *llis):
412        '''Interpolate the height for one or several locations.
413
414           @arg llis: The location or locations (C{LatLon}, ... or
415                      C{LatLon}s).
416
417           @return: A single interpolated height (C{float}) or a list
418                    or tuple of interpolated heights (C{float}s).
419
420           @raise HeightError: Insufficient number of B{C{llis}},
421                               an invalid B{C{lli}} or an L{fidw}
422                               issue.
423        '''
424        _as, xis, yis, _ = _axyllis4(tuple, llis, off=False)
425        return _as(map(self._hIDW, xis, yis))
426
427    def _datum_setter(self, datum, knots):
428        '''(INTERNAL) Set the datum.
429
430           @raise TypeError: Invalid B{C{datum}}.
431        '''
432        d = datum or getattr(knots[0], _datum_, datum)
433        if d not in (None, self._datum):
434            self._datum = _ellipsoidal_datum(d, name=self.name)
435
436    def _distances(self, x, y):  # PYCHOK unused (x, y) radians
437        '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
438        '''
439        notOverloaded(self, x, y)
440
441    def _distances_angular_(self, func_, x, y):
442        # return angular distances from func_
443        for xk, yk in zip(self._xs, self._ys):
444            r, _ = unrollPI(xk, x, wrap=self._wrap)
445            yield func_(yk, y, r)
446
447    def _distances_angular_datum_(self, func_, x, y):
448        # return angular distances from func_
449        for xk, yk in zip(self._xs, self._ys):
450            r, _ = unrollPI(xk, x, wrap=self._wrap)
451            yield func_(yk, y, r, datum=self._datum)
452
453    def _hIDW(self, x, y):
454        # interpolate height at (x, y) radians or degrees
455        try:
456            ds = self._distances(x, y)
457            return fidw(self._hs, ds, beta=self._beta)
458        except (TypeError, ValueError) as x:
459            raise HeightError(str(x))
460
461    @property
462    def beta(self):
463        '''Get the inverse distance power (C{int}).
464        '''
465        return self._beta
466
467    @beta.setter  # PYCHOK setter!
468    def beta(self, beta):
469        '''Set the inverse distance power.
470
471           @arg beta: New inverse distance power (C{int} 1, 2, or 3).
472
473           @raise HeightError: Invalid B{C{beta}}.
474        '''
475        self._beta = Int_(beta=beta, Error=HeightError, low=1, high=3)
476
477    def height(self, lats, lons):
478        '''Interpolate the height for one or several lat-/longitudes.
479
480           @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
481           @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
482
483           @return: A single interpolated height (C{float}) or a list of
484                    interpolated heights (C{float}s).
485
486           @raise HeightError: Insufficient or non-matching number of
487                               B{C{lats}} and B{C{lons}} or an L{fidw}
488                               issue.
489        '''
490        return _HeightBase._height(self, lats, lons)
491
492
493class HeightIDWcosineAndoyerLambert(_HeightIDW):
494    '''Height interpolator using U{Inverse Distance Weighting
495       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
496       and the I{angular} distance in C{radians} from function
497       L{cosineAndoyerLambert_}.
498
499       @see: L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWdistanceTo},
500             L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas},
501             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
502             geostatistics/Inverse-Distance-Weighting/index.html>} and
503             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
504             shepard_interp_2d/shepard_interp_2d.html>}.
505    '''
506    _datum = _WGS84
507    _wrap  =  False
508
509    def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
510        '''New L{HeightIDWcosineAndoyerLambert} interpolator.
511
512           @arg knots: The points with known height (C{LatLon}s).
513           @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
514                         and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
515                         L{Ellipsoid2} or L{a_f2Tuple}).
516           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
517           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
518           @kwarg name: Optional name for this height interpolator (C{str}).
519
520           @raise HeightError: Insufficient number of B{C{knots}} or
521                               an invalid B{C{knot}} or B{C{beta}}.
522
523           @raise TypeError: Invalid B{C{datum}}.
524        '''
525        _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap)
526        self._datum_setter(datum, knots)
527
528    def _distances(self, x, y):  # (x, y) radians
529        return self._distances_angular_datum_(cosineAndoyerLambert_, x, y)
530
531    if _FOR_DOCS:
532        __call__ = _HeightIDW.__call__
533        height   = _HeightIDW.height
534
535
536class HeightIDWcosineForsytheAndoyerLambert(_HeightIDW):
537    '''Height interpolator using U{Inverse Distance Weighting
538       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
539       and the I{angular} distance in C{radians} from function
540       L{cosineForsytheAndoyerLambert_}.
541
542       @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWdistanceTo},
543             L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas},
544             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
545             geostatistics/Inverse-Distance-Weighting/index.html>} and
546             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
547             shepard_interp_2d/shepard_interp_2d.html>}.
548    '''
549    _datum = _WGS84
550    _wrap  =  False
551
552    def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
553        '''New L{HeightIDWcosineForsytheAndoyerLambert} interpolator.
554
555           @arg knots: The points with known height (C{LatLon}s).
556           @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
557                         and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
558                         L{Ellipsoid2} or L{a_f2Tuple}).
559           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
560           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
561           @kwarg name: Optional name for this height interpolator (C{str}).
562
563           @raise HeightError: Insufficient number of B{C{knots}} or
564                               an invalid B{C{knot}} or B{C{beta}}.
565
566           @raise TypeError: Invalid B{C{datum}}.
567        '''
568        _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap)
569        self._datum_setter(datum, knots)
570
571    def _distances(self, x, y):  # (x, y) radians
572        return self._distances_angular_datum_(cosineForsytheAndoyerLambert_, x, y)
573
574    if _FOR_DOCS:
575        __call__ = _HeightIDW.__call__
576        height   = _HeightIDW.height
577
578
579class HeightIDWcosineLaw(_HeightIDW):
580    '''Height interpolator using U{Inverse Distance Weighting
581       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
582       and the I{angular} distance in C{radians} from function L{cosineLaw_}.
583
584       @see: L{HeightIDWequirectangular}, L{HeightIDWeuclidean},
585             L{HeightIDWflatPolar}, L{HeightIDWhaversine}, L{HeightIDWvincentys},
586             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
587             geostatistics/Inverse-Distance-Weighting/index.html>} and
588             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
589             shepard_interp_2d/shepard_interp_2d.html>}.
590
591       @note: See note at function L{vincentys_}.
592    '''
593    _wrap = False
594
595    def __init__(self, knots, beta=2, wrap=False, name=NN):
596        '''New L{HeightIDWcosineLaw} interpolator.
597
598           @arg knots: The points with known height (C{LatLon}s).
599           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
600           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
601           @kwarg name: Optional name for this height interpolator (C{str}).
602
603           @raise HeightError: Insufficient number of B{C{knots}} or
604                               an invalid B{C{knot}} or B{C{beta}}.
605        '''
606        _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap)
607
608    def _distances(self, x, y):  # (x, y) radians
609        return self._distances_angular_(cosineLaw_, x, y)
610
611    if _FOR_DOCS:
612        __call__ = _HeightIDW.__call__
613        height   = _HeightIDW.height
614
615
616class HeightIDWdistanceTo(_HeightIDW):
617    '''Height interpolator using U{Inverse Distance Weighting
618       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
619       and the distance from the points' C{LatLon.distanceTo} method,
620       conventionally in C{meter}.
621
622       @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert},
623             L{HeightIDWflatPolar}, L{HeightIDWkarney}, L{HeightIDWthomas},
624             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
625             geostatistics/Inverse-Distance-Weighting/index.html>} and
626             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
627             shepard_interp_2d/shepard_interp_2d.html>}.
628    '''
629    _distanceTo_kwds = {}
630    _ks              = ()
631
632    def __init__(self, knots, beta=2, name=NN, **distanceTo_kwds):
633        '''New L{HeightIDWdistanceTo} interpolator.
634
635           @arg knots: The points with known height (C{LatLon}s).
636           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
637           @kwarg name: Optional name for this height interpolator (C{str}).
638           @kwarg distanceTo_kwds: Optional keyword arguments for the
639                                   B{C{points}}' C{LatLon.distanceTo}
640                                   method.
641
642           @raise HeightError: Insufficient number of B{C{knots}} or
643                               an invalid B{C{knot}} or B{C{beta}}.
644
645           @raise ImportError: Package U{geographiclib
646                  <https://PyPI.org/project/geographiclib>} missing
647                  iff B{C{points}} are L{ellipsoidalKarney.LatLon}s.
648
649           @note: All B{C{points}} I{must} be instances of the same
650                  ellipsoidal or spherical C{LatLon} class, I{not
651                  checked however}.
652        '''
653        n, self._ks = len2(knots)
654        if n < self._kmin:
655            raise _insufficientError(self._kmin, knots=n)
656        for i, k in enumerate(self._ks):
657            if not callable(getattr(k, _distanceTo_, None)):
658                i = Fmt.SQUARE(_knots_, i)
659                raise HeightError(i, k, txt=_distanceTo_)
660
661        # use knots[0] class and datum to create
662        # compatible points in _HeightBase._height
663        # instead of class LatLon_ and datum None
664        self._datum = self._ks[0].datum
665        self._LLis  = self._ks[0].classof
666
667        self.beta = beta
668        if name:
669            self.name = name
670        if distanceTo_kwds:
671            self._distanceTo_kwds = distanceTo_kwds
672
673    def __call__(self, *llis):
674        '''Interpolate the height for one or several locations.
675
676           @arg llis: The location or locations (C{LatLon}, ... or
677                      C{LatLon}s).
678
679           @return: A single interpolated height (C{float}) or a list
680                    or tuple of interpolated heights (C{float}s).
681
682           @raise HeightError: Insufficient number of B{C{llis}},
683                               an invalid B{C{lli}} or an L{fidw}
684                               issue.
685        '''
686        _as, llis = _allis2(llis)
687        return _as(map(self._hIDW, llis))
688
689    if _FOR_DOCS:
690        height = _HeightIDW.height
691
692    def _hIDW(self, lli):  # PYCHOK expected
693        # interpolate height at point lli
694        try:
695            kwds = self._distanceTo_kwds
696            ds = (k.distanceTo(lli, **kwds) for k in self._ks)
697            return fidw(self._hs, ds, beta=self._beta)
698        except (TypeError, ValueError) as x:
699            raise HeightError(str(x))
700
701    @property  # NOT _RO, no caching
702    def _hs(self):  # see HeightIDWkarney
703        for k in self._ks:
704            yield k.height
705
706
707class HeightIDWequirectangular(_HeightIDW):
708    '''Height interpolator using U{Inverse Distance Weighting
709       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
710       the I{angular} distance in C{radians squared} like function
711       L{equirectangular_}.
712
713       @see: L{HeightIDWeuclidean}, L{HeightIDWflatPolar},
714             L{HeightIDWhaversine}, L{HeightIDWvincentys},
715             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
716             geostatistics/Inverse-Distance-Weighting/index.html>} and
717             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
718             shepard_interp_2d/shepard_interp_2d.html>}.
719    '''
720    _adjust = True
721    _wrap   = False
722
723    def __init__(self, knots, adjust=True, wrap=False, name=NN):
724        '''New L{HeightIDWequirectangular} interpolator.
725
726           @arg knots: The points with known height (C{LatLon}s).
727           @kwarg adjust: Adjust the wrapped, unrolled longitudinal
728                          delta by the cosine of the mean latitude (C{bool}).
729           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
730           @kwarg name: Optional name for this height interpolator (C{str}).
731
732           @raise HeightError: Insufficient number of B{C{knots}} or
733                               an invalid B{C{knot}}.
734        '''
735        _HeightIDW.__init__(self, knots, beta=1, name=name, wrap=wrap,
736                                                          adjust=adjust)
737
738    def _distances(self, x, y):  # (x, y) radians**2
739        for xk, yk in zip(self._xs, self._ys):
740            d, _ = unrollPI(xk, x, wrap=self._wrap)
741            if self._adjust:
742                d *= _scale_rad(yk, y)
743            yield hypot2(d, yk - y)  # like equirectangular_ distance2
744
745    if _FOR_DOCS:
746        __call__ = _HeightIDW.__call__
747        height   = _HeightIDW.height
748
749
750class HeightIDWeuclidean(_HeightIDW):
751    '''Height interpolator using U{Inverse Distance Weighting
752       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
753       and the I{angular} distance in C{radians} from function L{euclidean_}.
754
755       @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular},
756             L{HeightIDWhaversine}, L{HeightIDWvincentys},
757             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
758             geostatistics/Inverse-Distance-Weighting/index.html>} and
759             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
760             shepard_interp_2d/shepard_interp_2d.html>}.
761    '''
762    _adjust = True
763
764    def __init__(self, knots, adjust=True, beta=2, name=NN):
765        '''New L{HeightIDWeuclidean} interpolator.
766
767           @arg knots: The points with known height (C{LatLon}s).
768           @kwarg adjust: Adjust the longitudinal delta by the cosine
769                          of the mean latitude for B{C{adjust}}=C{True}.
770           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
771           @kwarg name: Optional name for this height interpolator (C{str}).
772
773           @raise HeightError: Insufficient number of B{C{knots}} or
774                               an invalid B{C{knot}} or B{C{beta}}.
775        '''
776        _HeightIDW.__init__(self, knots, beta=beta, name=name, adjust=adjust)
777
778    def _distances(self, x, y):  # (x, y) radians
779        for xk, yk in zip(self._xs, self._ys):
780            yield euclidean_(yk, y, xk - x, adjust=self._adjust)
781
782    if _FOR_DOCS:
783        __call__ = _HeightIDW.__call__
784        height   = _HeightIDW.height
785
786
787class HeightIDWflatLocal(_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 squared} like function
791       L{flatLocal_}/L{hubeny_}.
792
793       @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert},
794             L{HeightIDWdistanceTo}, L{HeightIDWhubeny}, L{HeightIDWthomas},
795             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
796             geostatistics/Inverse-Distance-Weighting/index.html>} and
797             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
798             shepard_interp_2d/shepard_interp_2d.html>}.
799    '''
800    _datum = _WGS84
801    _wrap  =  False
802
803    def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
804        '''New L{HeightIDWflatLocal}/L{HeightIDWhubeny} interpolator.
805
806           @arg knots: The points with known height (C{LatLon}s).
807           @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
808                         and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
809                         L{Ellipsoid2} or L{a_f2Tuple}).
810           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
811           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
812           @kwarg name: Optional name for this height interpolator (C{str}).
813
814           @raise HeightError: Insufficient number of B{C{knots}} or
815                               an invalid B{C{knot}} or B{C{beta}}.
816
817           @raise TypeError: Invalid B{C{datum}}.
818        '''
819        _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap)
820        self._datum_setter(datum, knots)
821
822    def _distances(self, x, y):  # (x, y) radians
823        _r2_ = self._datum.ellipsoid._hubeny2_
824        return self._distances_angular_(_r2_, x, y)  # radians**2
825
826    if _FOR_DOCS:
827        __call__ = _HeightIDW.__call__
828        height   = _HeightIDW.height
829
830
831class HeightIDWflatPolar(_HeightIDW):
832    '''Height interpolator using U{Inverse Distance Weighting
833       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
834       and the I{angular} distance in C{radians} from function L{flatPolar_}.
835
836       @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular},
837             L{HeightIDWeuclidean}, L{HeightIDWhaversine}, L{HeightIDWvincentys},
838             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
839             geostatistics/Inverse-Distance-Weighting/index.html>} and
840             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
841             shepard_interp_2d/shepard_interp_2d.html>}.
842    '''
843    _wrap = False
844
845    def __init__(self, knots, beta=2, wrap=False, name=NN):
846        '''New L{HeightIDWflatPolar} interpolator.
847
848           @arg knots: The points with known height (C{LatLon}s).
849           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
850           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
851           @kwarg name: Optional name for this height interpolator (C{str}).
852
853           @raise HeightError: Insufficient number of B{C{knots}} or
854                               an invalid B{C{knot}} or B{C{beta}}.
855        '''
856        _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap)
857
858    def _distances(self, x, y):  # (x, y) radians
859        return self._distances_angular_(flatPolar_, x, y)
860
861    if _FOR_DOCS:
862        __call__ = _HeightIDW.__call__
863        height   = _HeightIDW.height
864
865
866class HeightIDWhaversine(_HeightIDW):
867    '''Height interpolator using U{Inverse Distance Weighting
868       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
869       and the I{angular} distance in C{radians} from function L{haversine_}.
870
871       @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular}, L{HeightIDWeuclidean},
872             L{HeightIDWflatPolar}, L{HeightIDWvincentys},
873             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
874             geostatistics/Inverse-Distance-Weighting/index.html>} and
875             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
876             shepard_interp_2d/shepard_interp_2d.html>}.
877
878       @note: See note at function L{vincentys_}.
879    '''
880    _wrap = False
881
882    def __init__(self, knots, beta=2, wrap=False, name=NN):
883        '''New L{HeightIDWhaversine} interpolator.
884
885           @arg knots: The points with known height (C{LatLon}s).
886           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
887           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
888           @kwarg name: Optional name for this height interpolator (C{str}).
889
890           @raise HeightError: Insufficient number of B{C{knots}} or
891                               an  B{C{knot}} or B{C{beta}}.
892        '''
893        _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap)
894
895    def _distances(self, x, y):  # (x, y) radians
896        return self._distances_angular_(haversine_, x, y)
897
898    if _FOR_DOCS:
899        __call__ = _HeightIDW.__call__
900        height   = _HeightIDW.height
901
902
903class HeightIDWhubeny(HeightIDWflatLocal):  # for Karl Hubeny
904    if _FOR_DOCS:
905        __doc__  = HeightIDWflatLocal.__doc__
906        __init__ = HeightIDWflatLocal.__init__
907        __call__ = HeightIDWflatLocal.__call__
908        height   = HeightIDWflatLocal.height
909
910
911class HeightIDWkarney(_HeightIDW):
912    '''Height interpolator using U{Inverse Distance Weighting
913       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW) and
914       the I{angular} distance in C{degrees} from I{Karney}'s
915       U{geographiclib<https://PyPI.org/project/geographiclib>} U{Geodesic
916       <https://GeographicLib.SourceForge.io/html/python/code.html>}
917       Inverse method.
918
919       @see: L{HeightIDWcosineAndoyerLambert},
920             L{HeightIDWcosineForsytheAndoyerLambert}, L{HeightIDWdistanceTo},
921             L{HeightIDWflatLocal}, L{HeightIDWhubeny}, L{HeightIDWthomas},
922             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
923             geostatistics/Inverse-Distance-Weighting/index.html>} and
924             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
925             shepard_interp_2d/shepard_interp_2d.html>}.
926    '''
927    _datum    = _WGS84
928    _Inverse1 =  None
929    _ks       = ()
930    _wrap     =  False
931
932    def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
933        '''New L{HeightIDWkarney} interpolator.
934
935           @arg knots: The points with known height (C{LatLon}s).
936           @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
937                         and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
938                         L{Ellipsoid2} or L{a_f2Tuple}).
939           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
940           @kwarg wrap: Wrap and L{unroll180} longitudes (C{bool}).
941           @kwarg name: Optional name for this height interpolator (C{str}).
942
943           @raise HeightError: Insufficient number of B{C{knots}} or
944                               an invalid B{C{knot}}, B{C{datum}} or
945                               B{C{beta}}.
946
947           @raise ImportError: Package U{geographiclib
948                  <https://PyPI.org/project/geographiclib>} missing.
949
950           @raise TypeError: Invalid B{C{datum}}.
951        '''
952        n, self._ks = len2(knots)
953        if n < self._kmin:
954            raise _insufficientError(self._kmin, knots=n)
955        self._datum_setter(datum, self._ks)
956        self._Inverse1 = self.datum.ellipsoid.geodesic.Inverse1
957
958        self.beta = beta
959        if wrap:
960            self._wrap = True
961        if name:
962            self.name = name
963
964    def _distances(self, x, y):  # (x, y) degrees
965        for k in self._ks:
966            # non-negative I{angular} distance in C{degrees}
967            yield self._Inverse1(y, x, k.lat, k.lon, wrap=self._wrap)
968
969    @property  # NOT _RO, no caching
970    def _hs(self):  # see HeightIDWdistanceTo
971        for k in self._ks:
972            yield k.height
973
974    def __call__(self, *llis):
975        '''Interpolate the height for one or several locations.
976
977           @arg llis: The location or locations (C{LatLon}, ... or
978                      C{LatLon}s).
979
980           @return: A single interpolated height (C{float}) or a list
981                    or tuple of interpolated heights (C{float}s).
982
983           @raise HeightError: Insufficient number of B{C{llis}},
984                               an invalid B{C{lli}} or an L{fidw}
985                               issue.
986        '''
987        def _xy2(lls):
988            try:  # like _xyhs above, but keeping degrees
989                for i, ll in enumerate(lls):
990                    yield ll.lon, ll.lat
991            except AttributeError as x:
992                i = Fmt.SQUARE(llis=i)
993                raise HeightError(i, ll, txt=str(x))
994
995        _as, llis = _allis2(llis)
996        return _as(map(self._hIDW, *zip(*_xy2(llis))))
997
998    if _FOR_DOCS:
999        height = _HeightIDW.height
1000
1001
1002class HeightIDWthomas(_HeightIDW):
1003    '''Height interpolator using U{Inverse Distance Weighting
1004       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
1005       and the I{angular} distance in C{radians} from function L{thomas_}.
1006
1007       @see: L{HeightIDWcosineAndoyerLambert}, L{HeightIDWcosineForsytheAndoyerLambert},
1008             L{HeightIDWdistanceTo}, L{HeightIDWflatLocal}, L{HeightIDWhubeny},
1009             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
1010             geostatistics/Inverse-Distance-Weighting/index.html>} and
1011             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
1012             shepard_interp_2d/shepard_interp_2d.html>}.
1013    '''
1014    _datum = _WGS84
1015    _wrap  =  False
1016
1017    def __init__(self, knots, datum=None, beta=2, wrap=False, name=NN):
1018        '''New L{HeightIDWthomas} interpolator.
1019
1020           @arg knots: The points with known height (C{LatLon}s).
1021           @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
1022                         and first B{C{knots}}' datum (L{Datum}, L{Ellipsoid},
1023                         L{Ellipsoid2} or L{a_f2Tuple}).
1024           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
1025           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
1026           @kwarg name: Optional name for this height interpolator (C{str}).
1027
1028           @raise HeightError: Insufficient number of B{C{knots}} or
1029                               an invalid B{C{knot}} or B{C{beta}}.
1030
1031           @raise TypeError: Invalid B{C{datum}}.
1032        '''
1033        _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap)
1034        self._datum_setter(datum, knots)
1035
1036    def _distances(self, x, y):  # (x, y) radians
1037        return self._distances_angular_datum_(thomas_, x, y)
1038
1039    if _FOR_DOCS:
1040        __call__ = _HeightIDW.__call__
1041        height   = _HeightIDW.height
1042
1043
1044class HeightIDWvincentys(_HeightIDW):
1045    '''Height interpolator using U{Inverse Distance Weighting
1046       <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW)
1047       and the I{angular} distance in C{radians} from function L{vincentys_}.
1048
1049       @see: L{HeightIDWcosineLaw}, L{HeightIDWequirectangular},
1050             L{HeightIDWeuclidean}, L{HeightIDWflatPolar}, L{HeightIDWhaversine},
1051             U{IDW<https://www.Geo.FU-Berlin.DE/en/v/soga/Geodata-analysis/
1052             geostatistics/Inverse-Distance-Weighting/index.html>} and
1053             U{SHEPARD_INTERP_2D<https://People.SC.FSU.edu/~jburkardt/c_src/
1054             shepard_interp_2d/shepard_interp_2d.html>}.
1055
1056       @note: See note at function L{vincentys_}.
1057    '''
1058    _wrap = False
1059
1060    def __init__(self, knots, beta=2, wrap=False, name=NN):
1061        '''New L{HeightIDWvincentys} interpolator.
1062
1063           @arg knots: The points with known height (C{LatLon}s).
1064           @kwarg beta: Inverse distance power (C{int} 1, 2, or 3).
1065           @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
1066           @kwarg name: Optional name for this height interpolator (C{str}).
1067
1068           @raise HeightError: Insufficient number of B{C{knots}} or
1069                               an invalid B{C{knot}} or B{C{beta}}.
1070        '''
1071        _HeightIDW.__init__(self, knots, beta=beta, name=name, wrap=wrap)
1072
1073    def _distances(self, x, y):  # (x, y) radians
1074        return self._distances_angular_(vincentys_, x, y)
1075
1076    if _FOR_DOCS:
1077        __call__ = _HeightIDW.__call__
1078        height   = _HeightIDW.height
1079
1080
1081class HeightLSQBiSpline(_HeightBase):
1082    '''Height interpolator using C{SciPy} U{LSQSphereBivariateSpline
1083       <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
1084       interpolate.LSQSphereBivariateSpline.html>}.
1085    '''
1086    _kmin = 16  # k = 3, always
1087
1088    def __init__(self, knots, weight=None, name=NN):
1089        '''New L{HeightLSQBiSpline} interpolator.
1090
1091           @arg knots: The points with known height (C{LatLon}s).
1092           @kwarg weight: Optional weight or weights for each B{C{knot}}
1093                          (C{scalar} or C{scalar}s).
1094           @kwarg name: Optional name for this height interpolator (C{str}).
1095
1096           @raise HeightError: Insufficient number of B{C{knots}} or
1097                               an invalid B{C{knot}} or B{C{weight}}.
1098
1099           @raise LenError: Unequal number of B{C{knots}} and B{C{weight}}s.
1100
1101           @raise ImportError: Package C{numpy} or C{scipy} not found
1102                               or not installed.
1103
1104           @raise SciPyError: A C{LSQSphereBivariateSpline} issue.
1105
1106           @raise SciPyWarning: A C{LSQSphereBivariateSpline} warning
1107                                as exception.
1108        '''
1109        np, spi = self._NumSciPy()
1110
1111        xs, ys, hs = self._xyhs3(knots)
1112        n = len(hs)
1113
1114        w = weight
1115        if isscalar(w):
1116            w = float(w)
1117            if w <= 0:
1118                raise HeightError(weight=w)
1119            w = [w] * n
1120        elif w is not None:
1121            m, w = len2(w)
1122            if m != n:
1123                raise LenError(HeightLSQBiSpline, weight=m, knots=n)
1124            w = map2(float, w)
1125            m = min(w)
1126            if m <= 0:
1127                w = Fmt.SQUARE(weight=w.find(m))
1128                raise HeightError(w, m)
1129        try:
1130            T = 1.0e-4  # like SciPy example
1131            ps = np.array(_ordedup(xs, T, PI2 - T))
1132            ts = np.array(_ordedup(ys, T, PI  - T))
1133            self._ev = spi.LSQSphereBivariateSpline(ys, xs, hs,
1134                                                    ts, ps, eps=EPS, w=w).ev
1135        except Exception as x:
1136            raise _SciPyIssue(x)
1137
1138        if name:
1139            self.name = name
1140
1141    def __call__(self, *llis):
1142        '''Interpolate the height for one or several locations.
1143
1144           @arg llis: The location or locations (C{LatLon}, ... or
1145                      C{LatLon}s).
1146
1147           @return: A single interpolated height (C{float}) or a list
1148                    or tuple of interpolated heights (C{float}s).
1149
1150           @raise HeightError: Insufficient number of B{C{llis}} or
1151                               an invalid B{C{lli}}.
1152
1153           @raise SciPyError: A C{LSQSphereBivariateSpline} issue.
1154
1155           @raise SciPyWarning: A C{LSQSphereBivariateSpline} warning
1156                                as exception.
1157        '''
1158        return _HeightBase._eval(self, llis)
1159
1160    def height(self, lats, lons):
1161        '''Interpolate the height for one or several lat-/longitudes.
1162
1163           @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
1164           @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
1165
1166           @return: A single interpolated height (C{float}) or a list of
1167                    interpolated heights (C{float}s).
1168
1169           @raise HeightError: Insufficient or non-matching number of
1170                               B{C{lats}} and B{C{lons}}.
1171
1172           @raise SciPyError: A C{LSQSphereBivariateSpline} issue.
1173
1174           @raise SciPyWarning: A C{LSQSphereBivariateSpline} warning
1175                                as exception.
1176        '''
1177        return _HeightBase._height(self, lats, lons)
1178
1179
1180class HeightSmoothBiSpline(_HeightBase):
1181    '''Height interpolator using C{SciPy} U{SmoothSphereBivariateSpline
1182       <https://docs.SciPy.org/doc/scipy/reference/generated/scipy.
1183       interpolate.SmoothSphereBivariateSpline.html>}.
1184    '''
1185    _kmin = 16  # k = 3, always
1186
1187    def __init__(self, knots, s=4, name=NN):
1188        '''New L{HeightSmoothBiSpline} interpolator.
1189
1190           @arg knots: The points with known height (C{LatLon}s).
1191           @kwarg s: The spline smoothing factor (C{4}).
1192           @kwarg name: Optional name for this height interpolator (C{str}).
1193
1194           @raise HeightError: Insufficient number of B{C{knots}} or
1195                               an invalid B{C{knot}} or B{C{s}}.
1196
1197           @raise ImportError: Package C{numpy} or C{scipy} not found
1198                               or not installed.
1199
1200           @raise SciPyError: A C{SmoothSphereBivariateSpline} issue.
1201
1202           @raise SciPyWarning: A C{SmoothSphereBivariateSpline} warning
1203                                as exception.
1204        '''
1205        _, spi = self._NumSciPy()
1206
1207        s = Int_(s, name='smoothing', Error=HeightError, low=4)
1208
1209        xs, ys, hs = self._xyhs3(knots)
1210        try:
1211            self._ev = spi.SmoothSphereBivariateSpline(ys, xs, hs,
1212                                                       eps=EPS, s=s).ev
1213        except Exception as x:
1214            raise _SciPyIssue(x)
1215
1216        if name:
1217            self.name = name
1218
1219    def __call__(self, *llis):
1220        '''Interpolate the height for one or several locations.
1221
1222           @arg llis: The location or locations (C{LatLon}, ... or
1223                      C{LatLon}s).
1224
1225           @return: A single interpolated height (C{float}) or a list
1226                    or tuple of interpolated heights (C{float}s).
1227
1228           @raise HeightError: Insufficient number of B{C{llis}} or
1229                               an invalid B{C{lli}}.
1230
1231           @raise SciPyError: A C{SmoothSphereBivariateSpline} issue.
1232
1233           @raise SciPyWarning: A C{SmoothSphereBivariateSpline} warning
1234                                as exception.
1235        '''
1236        return _HeightBase._eval(self, llis)
1237
1238    def height(self, lats, lons):
1239        '''Interpolate the height for one or several lat-/longitudes.
1240
1241           @arg lats: Latitude or latitudes (C{degrees} or C{degrees}s).
1242           @arg lons: Longitude or longitudes (C{degrees} or C{degrees}s).
1243
1244           @return: A single interpolated height (C{float}) or a list of
1245                    interpolated heights (C{float}s).
1246
1247           @raise HeightError: Insufficient or non-matching number of
1248                               B{C{lats}} and B{C{lons}}.
1249
1250           @raise SciPyError: A C{SmoothSphereBivariateSpline} issue.
1251
1252           @raise SciPyWarning: A C{SmoothSphereBivariateSpline} warning
1253                                as exception.
1254        '''
1255        return _HeightBase._height(self, lats, lons)
1256
1257
1258__all__ += _ALL_DOCS(_HeightBase)
1259
1260# **) MIT License
1261#
1262# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved.
1263#
1264# Permission is hereby granted, free of charge, to any person obtaining a
1265# copy of this software and associated documentation files (the "Software"),
1266# to deal in the Software without restriction, including without limitation
1267# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1268# and/or sell copies of the Software, and to permit persons to whom the
1269# Software is furnished to do so, subject to the following conditions:
1270#
1271# The above copyright notice and this permission notice shall be included
1272# in all copies or substantial portions of the Software.
1273#
1274# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1275# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1276# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
1277# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1278# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1279# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1280# OTHER DEALINGS IN THE SOFTWARE.
1281