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