1 2# -*- coding: utf-8 -*- 3 4u'''Fréchet distances. 5 6Classes L{Frechet}, L{FrechetDegrees}, L{FrechetRadians}, 7L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert}, 8L{FrechetCosineLaw}, L{FrechetDistanceTo}< L{FrechetEquirectangular}, 9L{FrechetEuclidean}, L{FrechetExact}, L{FrechetFlatLocal}, L{FrechetFlatPolar}, 10L{FrechetHaversine}, L{FrechetHubeny}, L{FrechetKarney}, L{FrechetThomas} 11and L{FrechetVincentys} to compute I{discrete} U{Fréchet 12<https://WikiPedia.org/wiki/Frechet_distance>} distances between two sets 13of C{LatLon}, C{NumPy}, C{tuples} or other types of points. 14 15Only L{FrechetDistanceTo} -iff used with L{ellipsoidalKarney.LatLon} 16points- and L{FrechetKarney} requires installation of I{Charles Karney}'s 17U{geographiclib<https://PyPI.org/project/geographiclib>}. 18 19Typical usage is as follows. First, create a C{Frechet} calculator 20from one set of C{LatLon} points. 21 22C{f = FrechetXyz(points1, ...)} 23 24Get the I{discrete} Fréchet distance to another set of C{LatLon} points 25by 26 27C{t6 = f.discrete(points2)} 28 29Or, use function C{frechet_} with a proper C{distance} function passed 30as keyword arguments as follows 31 32C{t6 = frechet_(points1, points2, ..., distance=...)}. 33 34In both cases, the returned result C{t6} is a L{Frechet6Tuple}. 35 36For C{(lat, lon, ...)} points in a C{NumPy} array or plain C{tuples}, 37wrap the points in a L{Numpy2LatLon} respectively L{Tuple2LatLon} 38instance, more details in the documentation thereof. 39 40For other points, create a L{Frechet} sub-class with the appropriate 41C{distance} method overloading L{Frechet.distance} as in this example. 42 43 >>> from pygeodesy import Frechet, hypot_ 44 >>> 45 >>> class F3D(Frechet): 46 >>> """Custom Frechet example. 47 >>> """ 48 >>> def distance(self, p1, p2): 49 >>> return hypot_(p1.x - p2.x, p1.y - p2.y, p1.z - p2.z) 50 >>> 51 >>> f3D = F3D(xyz1, ..., units="...") 52 >>> t6 = f3D.discrete(xyz2) 53 54Transcribed from the original U{Computing Discrete Fréchet Distance 55<https://www.kr.TUWien.ac.AT/staff/eiter/et-archive/cdtr9464.pdf>} by 56Eiter, T. and Mannila, H., 1994, April 25, Technical Report CD-TR 94/64, 57Information Systems Department/Christian Doppler Laboratory for Expert 58Systems, Technical University Vienna, Austria. 59 60This L{Frechet.discrete} implementation optionally generates intermediate 61points for each point set separately. For example, using keyword argument 62C{fraction=0.5} adds one additional point halfway between each pair of 63points. Or using C{fraction=0.1} interpolates nine additional points 64between each points pair. 65 66The L{Frechet6Tuple} attributes C{fi1} and/or C{fi2} will be I{fractional} 67indices of type C{float} if keyword argument C{fraction} is used. Otherwise, 68C{fi1} and/or C{fi2} are simply type C{int} indices into the respective 69points set. 70 71For example, C{fractional} index value 2.5 means an intermediate point 72halfway between points[2] and points[3]. Use function L{fractional} 73to obtain the intermediate point for a I{fractional} index in the 74corresponding set of points. 75 76The C{Fréchet} distance was introduced in 1906 by U{Maurice Fréchet 77<https://WikiPedia.org/wiki/Maurice_Rene_Frechet>}, see U{reference 78[6]<https://www.kr.TUWien.ac.AT/staff/eiter/et-archive/cdtr9464.pdf>}. 79It is a measure of similarity between curves that takes into account the 80location and ordering of the points. Therefore, it is often a better metric 81than the well-known C{Hausdorff} distance, see the L{hausdorff} module. 82''' 83 84from pygeodesy.basics import isscalar, _xinstanceof 85from pygeodesy.datums import Datum, _WGS84 86from pygeodesy.errors import _IsnotError, PointsError 87from pygeodesy.fmath import hypot2 88from pygeodesy.formy import cosineAndoyerLambert_, cosineForsytheAndoyerLambert_, \ 89 cosineLaw_, euclidean_, flatPolar_, haversine_, \ 90 thomas_, vincentys_, _scale_rad 91from pygeodesy.interns import EPS, EPS1, INF, NN, _datum_, _distanceTo_, _DOT_, \ 92 _n_, _points_, _units_ 93from pygeodesy.iters import points2 as _points2 94from pygeodesy.lazily import _ALL_LAZY, _FOR_DOCS 95from pygeodesy.named import _Named, _NamedTuple, notOverloaded, _Pass 96from pygeodesy.namedTuples import PhiLam2Tuple 97from pygeodesy.points import _fractional 98from pygeodesy.props import Property_RO, property_doc_ 99from pygeodesy.streprs import _boolkwds, Fmt 100from pygeodesy.units import FIx, Float, Number_, _Str_degrees, _Str_meter, \ 101 _Str_NN, _Str_radians, _Str_radians2, _xUnit, _xUnits 102from pygeodesy.utily import unrollPI 103 104from collections import defaultdict as _defaultdict 105from math import radians 106 107__all__ = _ALL_LAZY.frechet 108__version__ = '21.09.14' 109 110 111def _fraction(fraction, n): 112 f = 1 # int, no fractional indices 113 if fraction in (None, 1): 114 pass 115 elif not (isscalar(fraction) and EPS < fraction < EPS1 116 and (float(n) - fraction) < n): 117 raise FrechetError(fraction=fraction) 118 elif fraction < EPS1: 119 f = float(fraction) 120 return f 121 122 123class FrechetError(PointsError): 124 '''Fréchet issue. 125 ''' 126 pass 127 128 129class Frechet(_Named): 130 '''Frechet base class, requires method L{Frechet.distance} to 131 be overloaded. 132 ''' 133 _adjust = None # not applicable 134 _datum = None # not applicable 135 _f1 = 1 136 _n1 = 0 137 _ps1 = None 138 _units = _Str_NN # XXX Str to _Pass and for backward compatibility 139 _wrap = None # not applicable 140 141 def __init__(self, points, fraction=None, name=NN, units=NN, **wrap_adjust): 142 '''New C{Frechet...} calculator/interpolator. 143 144 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 145 L{Tuple2LatLon}[] or C{other}[]). 146 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 147 interpolate intermediate B{C{points}} or use 148 C{None}, C{0} or C{1} for no intermediate 149 B{C{points}} and no I{fractional} indices. 150 @kwarg name: Optional calculator/interpolator name (C{str}). 151 @kwarg units: Optional, the distance units (C{Unit} or C{str}). 152 @kwarg wrap_adjust: Optionally, C{wrap} and unroll longitudes, iff 153 applicable (C{bool}) and C{adjust} wrapped, 154 unrolled longitudinal delta by the cosine 155 of the mean latitude, iff applicable. 156 157 @raise FrechetError: Insufficient number of B{C{points}} or 158 invalid B{C{fraction}} or B{C{wrap}} or 159 B{C{ajust}} not applicable. 160 161 ''' 162 self._n1, self._ps1 = self._points2(points) 163 if fraction: 164 self.fraction = fraction 165 if name: 166 self.name = name 167 if units: # and not self.units: 168 self.units = units 169 if wrap_adjust: 170 _boolkwds(self, **wrap_adjust) 171 172 @Property_RO 173 def adjust(self): 174 '''Get the adjust setting (C{bool} or C{None} if not applicable). 175 ''' 176 return self._adjust 177 178 @Property_RO 179 def datum(self): 180 '''Get the datum (L{Datum} or C{None} if not applicable). 181 ''' 182 return self._datum 183 184 def _datum_setter(self, datum): 185 '''(INTERNAL) Set the datum. 186 ''' 187 d = datum or getattr(self._ps1[0], _datum_, datum) 188 if d and d != self.datum: # PYCHOK no cover 189 _xinstanceof(Datum, datum=d) 190 self._datum = d 191 192 def discrete(self, points, fraction=None): 193 '''Compute the C{forward, discrete Fréchet} distance. 194 195 @arg points: Second set of points (C{LatLon}[], L{Numpy2LatLon}[], 196 L{Tuple2LatLon}[] or C{other}[]). 197 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 198 interpolate intermediate B{C{points}} or use 199 C{None}, C{0} or C{1} for no intermediate 200 B{C{points}} and no I{fractional} indices. 201 202 @return: A L{Frechet6Tuple}C{(fd, fi1, fi2, r, n, units)}. 203 204 @raise FrechetError: Insufficient number of B{C{points}} or an 205 invalid B{C{point}} or B{C{fraction}}. 206 207 @raise RecursionError: Recursion depth exceeded, see U{sys.getrecursionlimit() 208 <https://docs.Python.org/3/library/sys.html#sys.getrecursionlimit>}. 209 ''' 210 n2, ps2 = self._points2(points) 211 212 f2 = _fraction(fraction, n2) 213 p2 = self.points_fraction if f2 < EPS1 else self.points_ # PYCHOK expected 214 215 f1 = self.fraction 216 p1 = self.points_fraction if f1 < EPS1 else self.points_ # PYCHOK expected 217 218 def dF(fi1, fi2): 219 return self.distance(p1(self._ps1, fi1), p2(ps2, fi2)) 220 221 try: 222 return _frechet_(self._n1, f1, n2, f2, dF, self.units) 223 except TypeError as x: 224 t = _DOT_(self.classname, self.discrete.__name__) 225 raise FrechetError(t, txt=str(x)) 226 227 def distance(self, point1, point2): # PYCHOK no cover 228 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}. 229 ''' 230 notOverloaded(self, point1, point2) 231 232 @property_doc_(''' the index fraction (C{float}).''') 233 def fraction(self): 234 '''Get the index fraction (C{float} or C{1}). 235 ''' 236 return self._f1 237 238 @fraction.setter # PYCHOK setter! 239 def fraction(self, fraction): 240 '''Set the index fraction (C{float} or C{1}). 241 242 @arg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 243 interpolate intermediate B{C{points}} or use 244 C{None}, C{0} or C{1} for no intermediate 245 B{C{points}} and no I{fractional} indices. 246 247 @raise FrechetError: Invalid B{C{fraction}}. 248 ''' 249 self._f1 = _fraction(fraction, self._n1) 250 251 def point(self, point): 252 '''Convert a point for the C{.distance} method. 253 254 @arg point: The point to convert ((C{LatLon}, L{Numpy2LatLon}, 255 L{Tuple2LatLon} or C{other}). 256 257 @return: The converted B{C{point}}. 258 ''' 259 return point # pass thru 260 261 def points_(self, points, i): 262 '''Get and convert a point for the C{.distance} method. 263 264 @arg points: The orignal B{C{points}} to convert ((C{LatLon}[], 265 L{Numpy2LatLon}[], L{Tuple2LatLon}[] or C{other}[]). 266 @arg i: The B{C{points}} index (C{int}). 267 268 @return: The converted B{C{points[i]}}. 269 ''' 270 return self.point(points[i]) 271 272 def _points2(self, points): 273 '''(INTERNAL) Check a set of points. 274 ''' 275 return _points2(points, closed=False, Error=FrechetError) 276 277 def points_fraction(self, points, fi): 278 '''Get and convert a I{fractional} point for the C{.distance} method. 279 280 @arg points: The orignal B{C{points}} to convert ((C{LatLon}[], 281 L{Numpy2LatLon}[], L{Tuple2LatLon}[] or C{other}[]). 282 @arg fi: The I{fractional} index in B{C{points}} (C{float} or C{int}). 283 284 @return: The interpolated, converted, intermediate B{C{points[fi]}}. 285 ''' 286 return self.point(_fractional(points, fi)) 287 288 @property_doc_(''' the distance units (C{Unit} or C{str}).''') 289 def units(self): 290 '''Get the distance units (C{Unit} or C{str}). 291 ''' 292 return self._units 293 294 @units.setter # PYCHOK setter! 295 def units(self, units): 296 '''Set the distance units. 297 298 @arg units: New units (C{Unit} or C{str}). 299 300 @raise TypeError: Invalid B{C{units}}. 301 ''' 302 self._units = _xUnits(units, Base=Float) 303 304 @Property_RO 305 def wrap(self): 306 '''Get the wrap setting (C{bool} or C{None} if not applicable). 307 ''' 308 return self._wrap 309 310 311class FrechetDegrees(Frechet): 312 '''L{Frechet} base class for distances in C{degrees} from 313 C{LatLon} points in C{degrees}. 314 ''' 315 _units = _Str_degrees 316 317 if _FOR_DOCS: 318 __init__ = Frechet.__init__ 319 discrete = Frechet.discrete 320 321 322class FrechetRadians(Frechet): 323 '''L{Frechet} base class for distances in C{radians} or C{radians 324 squared} from C{LatLon} points converted from C{degrees} to 325 C{radians}. 326 ''' 327 _units = _Str_radians 328 329 if _FOR_DOCS: 330 __init__ = Frechet.__init__ 331 discrete = Frechet.discrete 332 333 def point(self, point): 334 '''Convert C{(lat, lon)} point in degrees to C{(a, b)} 335 in radians. 336 337 @return: An L{PhiLam2Tuple}C{(phi, lam)}. 338 ''' 339 try: 340 return point.philam 341 except AttributeError: 342 return PhiLam2Tuple(radians(point.lat), radians(point.lon)) 343 344 345class FrechetCosineAndoyerLambert(FrechetRadians): 346 '''Compute the C{Frechet} distance based on the I{angular} distance 347 in C{radians} from function L{cosineAndoyerLambert_}. 348 349 @see: L{FrechetCosineForsytheAndoyerLambert}, L{FrechetDistanceTo}, 350 L{FrechetExact}, L{FrechetFlatLocal}, L{FrechetHubeny}, 351 L{FrechetThomas} and L{FrechetKarney}. 352 ''' 353 _datum = _WGS84 354 _wrap = False 355 356 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN): 357 '''New L{FrechetCosineAndoyerLambert} calculator/interpolator. 358 359 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 360 L{Tuple2LatLon}[] or C{other}[]). 361 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 362 and first B{C{points}}' datum (L{Datum}). 363 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}) 364 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 365 interpolate intermediate B{C{points}} or use 366 C{None}, C{0} or C{1} for no intermediate 367 B{C{points}} and no L{fractional} indices. 368 @kwarg name: Optional calculator/interpolator name (C{str}). 369 370 @raise FrechetError: Insufficient number of B{C{points}} or 371 invalid B{C{fraction}}. 372 373 @raise TypeError: Invalid B{C{datum}}. 374 ''' 375 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 376 wrap=wrap) 377 self._datum_setter(datum) 378 379 if _FOR_DOCS: 380 discrete = Frechet.discrete 381 382 def distance(self, p1, p2): 383 '''Return the L{cosineAndoyerLambert_} distance in C{radians}. 384 ''' 385 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 386 return cosineAndoyerLambert_(p2.phi, p1.phi, r, datum=self._datum) 387 388 389class FrechetCosineForsytheAndoyerLambert(FrechetRadians): 390 '''Compute the C{Frechet} distance based on the I{angular} distance 391 in C{radians} from function L{cosineForsytheAndoyerLambert_}. 392 393 @see: L{FrechetCosineAndoyerLambert}, L{FrechetDistanceTo}, 394 L{FrechetExact}, L{FrechetFlatLocal}, L{FrechetHubeny}, 395 L{FrechetThomas} and L{FrechetKarney}. 396 ''' 397 _datum = _WGS84 398 _wrap = False 399 400 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN): 401 '''New L{FrechetCosineForsytheAndoyerLambert} calculator/interpolator. 402 403 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 404 L{Tuple2LatLon}[] or C{other}[]). 405 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 406 and first B{C{points}}' datum (L{Datum}). 407 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}) 408 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 409 interpolate intermediate B{C{points}} or use 410 C{None}, C{0} or C{1} for no intermediate 411 B{C{points}} and no L{fractional} indices. 412 @kwarg name: Optional calculator/interpolator name (C{str}). 413 414 @raise FrechetError: Insufficient number of B{C{points}} or 415 invalid B{C{fraction}}. 416 417 @raise TypeError: Invalid B{C{datum}}. 418 ''' 419 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 420 wrap=wrap) 421 self._datum_setter(datum) 422 423 if _FOR_DOCS: 424 discrete = Frechet.discrete 425 426 def distance(self, p1, p2): 427 '''Return the L{cosineForsytheAndoyerLambert_} distance in C{radians}. 428 ''' 429 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 430 return cosineForsytheAndoyerLambert_(p2.phi, p1.phi, r, datum=self._datum) 431 432 433class FrechetCosineLaw(FrechetRadians): 434 '''Compute the C{Frechet} distance based on the I{angular} distance 435 in C{radians} from function L{cosineLaw_}. 436 437 @see: L{FrechetEquirectangular}, L{FrechetEuclidean}, 438 L{FrechetExact}, L{FrechetFlatPolar}, L{FrechetHaversine} 439 and L{FrechetVincentys}. 440 441 @note: See note at function L{vincentys_}. 442 ''' 443 _wrap = False 444 445 def __init__(self, points, wrap=False, fraction=None, name=NN): 446 '''New L{FrechetCosineLaw} calculator/interpolator. 447 448 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 449 L{Tuple2LatLon}[] or C{other}[]). 450 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}). 451 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 452 interpolate intermediate B{C{points}} or use 453 C{None}, C{0} or C{1} for no intermediate 454 B{C{points}} and no I{fractional} indices. 455 @kwarg name: Optional calculator/interpolator name (C{str}). 456 457 @raise FrechetError: Insufficient number of B{C{points}} or 458 invalid B{C{fraction}}. 459 ''' 460 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 461 wrap=wrap) 462 463 if _FOR_DOCS: 464 discrete = Frechet.discrete 465 466 def distance(self, p1, p2): 467 '''Return the L{cosineLaw_} distance in C{radians}. 468 ''' 469 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 470 return cosineLaw_(p2.phi, p1.phi, r) 471 472 473class FrechetDistanceTo(Frechet): 474 '''Compute the C{Frechet} distance based on the distance from the 475 points' C{LatLon.distanceTo} method, conventionally in C{meter}. 476 477 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert}, 478 L{FrechetExact}, L{FrechetFlatLocal}, L{FrechetHubeny}, L{FrechetThomas} 479 and L{FrechetKarney}. 480 ''' 481 _distanceTo_kwds = {} 482 _units = _Str_meter 483 484 def __init__(self, points, fraction=None, name=NN, **distanceTo_kwds): 485 '''New L{FrechetDistanceTo} calculator/interpolator. 486 487 @arg points: First set of points (C{LatLon}[]). 488 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 489 interpolate intermediate B{C{points}} or use 490 C{None}, C{0} or C{1} for no intermediate 491 B{C{points}} and no I{fractional} indices. 492 @kwarg name: Optional calculator/interpolator name (C{str}). 493 @kwarg distanceTo_kwds: Optional keyword arguments for the 494 B{C{points}}' C{LatLon.distanceTo} 495 method. 496 497 @raise FrechetError: Insufficient number of B{C{points}} or an 498 invalid B{C{point}} or B{C{fraction}}. 499 500 @raise ImportError: Package U{geographiclib 501 <https://PyPI.org/project/geographiclib>} missing 502 iff B{C{points}} are L{ellipsoidalKarney.LatLon}s. 503 504 @note: All B{C{points}} I{must} be instances of the same 505 ellipsoidal or spherical C{LatLon} class. 506 ''' 507 FrechetRadians.__init__(self, points, fraction=fraction, name=name) 508 if distanceTo_kwds: 509 self._distanceTo_kwds = distanceTo_kwds 510 511 if _FOR_DOCS: 512 discrete = Frechet.discrete 513 514 def distance(self, p1, p2): 515 '''Return the distance in C{meter}. 516 ''' 517 return p1.distanceTo(p2, **self._distanceTo_kwds) 518 519 def _points2(self, points): 520 '''(INTERNAL) Check a set of points. 521 ''' 522 np, ps = Frechet._points2(self, points) 523 for i, p in enumerate(ps): 524 if not callable(getattr(p, _distanceTo_, None)): 525 i = Fmt.SQUARE(_points_, i) 526 raise FrechetError(i, p, txt=_distanceTo_) 527 return np, ps 528 529 530class FrechetEquirectangular(FrechetRadians): 531 '''Compute the C{Frechet} distance based on the C{equirectangular} 532 distance in C{radians squared} like function L{equirectangular_}. 533 534 @see: L{FrechetCosineLaw}, L{FrechetEuclidean}, L{FrechetExact}, 535 L{FrechetFlatPolar}, L{FrechetHaversine} and L{FrechetVincentys}. 536 ''' 537 _adjust = True 538 _units = _Str_radians2 539 _wrap = False 540 541 def __init__(self, points, adjust=True, wrap=False, fraction=None, name=NN): 542 '''New L{FrechetEquirectangular} calculator/interpolator. 543 544 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 545 L{Tuple2LatLon}[] or C{other}[]). 546 @kwarg adjust: Adjust the wrapped, unrolled longitudinal 547 delta by the cosine of the mean latitude (C{bool}). 548 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}). 549 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 550 interpolate intermediate B{C{points}} or use 551 C{None}, C{0} or C{1} for no intermediate 552 B{C{points}} and no L{fractional} indices. 553 @kwarg name: Optional calculator/interpolator name (C{str}). 554 555 @raise FrechetError: Insufficient number of B{C{points}} or 556 invalid B{C{adjust}} or B{C{seed}}. 557 ''' 558 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 559 adjust=adjust, wrap=wrap) 560 561 if _FOR_DOCS: 562 discrete = Frechet.discrete 563 564 def distance(self, p1, p2): 565 '''Return the L{equirectangular_} distance in C{radians squared}. 566 ''' 567 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 568 if self._adjust: 569 r *= _scale_rad(p1.phi, p2.phi) 570 return hypot2(r, p2.phi - p1.phi) # like equirectangular_ d2 571 572 573class FrechetEuclidean(FrechetRadians): 574 '''Compute the C{Frechet} distance based on the C{Euclidean} 575 distance in C{radians} from function L{euclidean_}. 576 577 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular}, 578 L{FrechetExact}, L{FrechetFlatPolar}, L{FrechetHaversine} 579 and L{FrechetVincentys}. 580 ''' 581 _adjust = True 582 _wrap = True # fixed 583 584 def __init__(self, points, adjust=True, fraction=None, name=NN): 585 '''New L{FrechetEuclidean} calculator/interpolator. 586 587 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 588 L{Tuple2LatLon}[] or C{other}[]). 589 @kwarg adjust: Adjust the wrapped, unrolled longitudinal 590 delta by the cosine of the mean latitude (C{bool}). 591 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 592 interpolate intermediate B{C{points}} or use 593 C{None}, C{0} or C{1} for no intermediate 594 B{C{points}} and no L{fractional} indices. 595 @kwarg name: Optional calculator/interpolator name (C{str}). 596 597 @raise FrechetError: Insufficient number of B{C{points}} or 598 invalid B{C{fraction}}. 599 ''' 600 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 601 adjust=adjust) 602 603 if _FOR_DOCS: 604 discrete = Frechet.discrete 605 606 def distance(self, p1, p2): 607 '''Return the L{euclidean_} distance in C{radians}. 608 ''' 609 return euclidean_(p2.phi, p1.phi, p2.lam - p1.lam, adjust=self._adjust) 610 611 612class FrechetExact(FrechetDegrees): 613 '''Compute the C{Frechet} distance based on the I{angular} 614 distance in C{degrees} from method L{GeodesicExact}C{.Inverse}. 615 616 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert}, 617 L{FrechetDistanceTo}, L{FrechetFlatLocal}, L{FrechetHubeny}, 618 L{FrechetKarney} and L{FrechetThomas}. 619 ''' 620 _datum = _WGS84 621 _Inverse1 = None 622 _wrap = False 623 624 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN): 625 '''New L{FrechetExact} calculator/interpolator. 626 627 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 628 L{Tuple2LatLon}[] or C{other}[]). 629 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 630 and first B{C{points}}' datum (L{Datum}). 631 @kwarg wrap: Wrap and L{unroll180} longitudes (C{bool}) 632 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 633 interpolate intermediate B{C{points}} or use 634 C{None}, C{0} or C{1} for no intermediate 635 B{C{points}} and no L{fractional} indices. 636 @kwarg name: Optional calculator/interpolator name (C{str}). 637 638 @raise FrechetError: Insufficient number of B{C{points}} or 639 invalid B{C{fraction}}. 640 641 @raise TypeError: Invalid B{C{datum}}. 642 ''' 643 FrechetDegrees.__init__(self, points, fraction=fraction, name=name, 644 wrap=wrap) 645 self._datum_setter(datum) 646 self._Inverse1 = self.datum.ellipsoid.geodesicx.Inverse1 # note -x 647 648 if _FOR_DOCS: 649 discrete = Frechet.discrete 650 651 def distance(self, p1, p2): 652 '''Return the non-negative I{angular} distance in C{degrees}. 653 ''' 654 return self._Inverse1(p1.lat, p1.lon, p2.lat, p2.lon, wrap=self._wrap) 655 656 657class FrechetFlatLocal(FrechetRadians): 658 '''Compute the C{Frechet} distance based on the I{angular} distance 659 in C{radians squared} like function L{flatLocal_}/L{hubeny_}. 660 661 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert}, 662 L{FrechetDistanceTo}, L{FrechetExact}, L{FrechetHubeny}, 663 L{FrechetKarney} and L{FrechetThomas}. 664 ''' 665 _datum = _WGS84 666 _hubeny2_ = None 667 _units = _Str_radians2 668 _wrap = False 669 670 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN): 671 '''New L{FrechetFlatLocal}/L{FrechetHubeny} calculator/interpolator. 672 673 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 674 L{Tuple2LatLon}[] or C{other}[]). 675 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 676 and first B{C{points}}' datum (L{Datum}). 677 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}) 678 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 679 interpolate intermediate B{C{points}} or use 680 C{None}, C{0} or C{1} for no intermediate 681 B{C{points}} and no L{fractional} indices. 682 @kwarg name: Optional calculator/interpolator name (C{str}). 683 684 @raise FrechetError: Insufficient number of B{C{points}} or 685 invalid B{C{fraction}}. 686 687 @raise TypeError: Invalid B{C{datum}}. 688 ''' 689 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 690 wrap=wrap) 691 self._datum_setter(datum) 692 self._hubeny2_ = self.datum.ellipsoid._hubeny2_ 693 694 if _FOR_DOCS: 695 discrete = Frechet.discrete 696 697 def distance(self, p1, p2): 698 '''Return the L{flatLocal_}/L{hubeny_} distance in C{radians squared}. 699 ''' 700 d, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 701 return self._hubeny2_(p2.phi, p1.phi, d) 702 703 704class FrechetFlatPolar(FrechetRadians): 705 '''Compute the C{Frechet} distance based on the I{angular} distance 706 in C{radians} from function L{flatPolar_}. 707 708 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular}, 709 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetHaversine} 710 and L{FrechetVincentys}. 711 ''' 712 _wrap = False 713 714 def __init__(self, points, wrap=False, fraction=None, name=NN): 715 '''New L{FrechetFlatPolar} calculator/interpolator. 716 717 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 718 L{Tuple2LatLon}[] or C{other}[]). 719 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}). 720 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 721 interpolate intermediate B{C{points}} or use 722 C{None}, C{0} or C{1} for no intermediate 723 B{C{points}} and no I{fractional} indices. 724 @kwarg name: Optional calculator/interpolator name (C{str}). 725 726 @raise FrechetError: Insufficient number of B{C{points}} or 727 invalid B{C{fraction}}. 728 ''' 729 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 730 wrap=wrap) 731 732 if _FOR_DOCS: 733 discrete = Frechet.discrete 734 735 def distance(self, p1, p2): 736 '''Return the L{flatPolar_} distance in C{radians}. 737 ''' 738 d, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 739 return flatPolar_(p2.phi, p1.phi, d) 740 741 742class FrechetHaversine(FrechetRadians): 743 '''Compute the C{Frechet} distance based on the I{angular} 744 distance in C{radians} from function L{haversine_}. 745 746 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular}, 747 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetFlatPolar} 748 and L{FrechetVincentys}. 749 750 @note: See note at function L{vincentys_}. 751 ''' 752 _wrap = False 753 754 def __init__(self, points, wrap=False, fraction=None, name=NN): 755 '''New L{FrechetHaversine} calculator/interpolator. 756 757 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 758 L{Tuple2LatLon}[] or C{other}[]). 759 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}). 760 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 761 interpolate intermediate B{C{points}} or use 762 C{None}, C{0} or C{1} for no intermediate 763 B{C{points}} and no I{fractional} indices. 764 @kwarg name: Optional calculator/interpolator name (C{str}). 765 766 @raise FrechetError: Insufficient number of B{C{points}} or 767 invalid B{C{fraction}}. 768 ''' 769 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 770 wrap=wrap) 771 772 if _FOR_DOCS: 773 discrete = Frechet.discrete 774 775 def distance(self, p1, p2): 776 '''Return the L{haversine_} distance in C{radians}. 777 ''' 778 d, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 779 return haversine_(p2.phi, p1.phi, d) 780 781 782class FrechetHubeny(FrechetFlatLocal): # for Karl Hubeny 783 if _FOR_DOCS: 784 __doc__ = FrechetFlatLocal.__doc__ 785 __init__ = FrechetFlatLocal.__init__ 786 discrete = FrechetFlatLocal.discrete 787 distance = FrechetFlatLocal.discrete 788 789 790class FrechetKarney(FrechetExact): 791 '''Compute the C{Frechet} distance based on the I{angular} 792 distance in C{degrees} from I{Karney}'s U{geographiclib 793 <https://PyPI.org/project/geographiclib>} U{Geodesic 794 <https://GeographicLib.SourceForge.io/html/python/code.html>} 795 Inverse method. 796 797 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert}, 798 L{FrechetDistanceTo}, L{FrechetExact}, L{FrechetFlatLocal}, 799 L{FrechetHubeny} and L{FrechetThomas}. 800 ''' 801 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN): 802 '''New L{FrechetKarney} calculator/interpolator. 803 804 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 805 L{Tuple2LatLon}[] or C{other}[]). 806 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 807 and first B{C{points}}' datum (L{Datum}). 808 @kwarg wrap: Wrap and L{unroll180} longitudes (C{bool}) 809 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 810 interpolate intermediate B{C{points}} or use 811 C{None}, C{0} or C{1} for no intermediate 812 B{C{points}} and no L{fractional} indices. 813 @kwarg name: Optional calculator/interpolator name (C{str}). 814 815 @raise FrechetError: Insufficient number of B{C{points}} or 816 invalid B{C{fraction}}. 817 818 @raise ImportError: Package U{geographiclib 819 <https://PyPI.org/project/geographiclib>} missing. 820 821 @raise TypeError: Invalid B{C{datum}}. 822 ''' 823 FrechetDegrees.__init__(self, points, fraction=fraction, name=name, 824 wrap=wrap) 825 self._datum_setter(datum) 826 self._Inverse1 = self.datum.ellipsoid.geodesic.Inverse1 827 828 829class FrechetThomas(FrechetRadians): 830 '''Compute the C{Frechet} distance based on the I{angular} distance 831 in C{radians} from function L{thomas_}. 832 833 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert}, 834 L{FrechetDistanceTo}, L{FrechetExact}, L{FrechetFlatLocal}, 835 L{FrechetHubeny} and L{FrechetKarney}. 836 ''' 837 _datum = _WGS84 838 _wrap = False 839 840 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN): 841 '''New L{FrechetThomas} calculator/interpolator. 842 843 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 844 L{Tuple2LatLon}[] or C{other}[]). 845 @kwarg datum: Optional datum overriding the default C{Datums.WGS84} 846 and first B{C{points}}' datum (L{Datum}). 847 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}) 848 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 849 interpolate intermediate B{C{points}} or use 850 C{None}, C{0} or C{1} for no intermediate 851 B{C{points}} and no L{fractional} indices. 852 @kwarg name: Optional calculator/interpolator name (C{str}). 853 854 @raise FrechetError: Insufficient number of B{C{points}} or 855 invalid B{C{fraction}}. 856 857 @raise TypeError: Invalid B{C{datum}}. 858 ''' 859 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 860 wrap=wrap) 861 self._datum_setter(datum) 862 863 if _FOR_DOCS: 864 discrete = Frechet.discrete 865 866 def distance(self, p1, p2): 867 '''Return the L{thomas_} distance in C{radians}. 868 ''' 869 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 870 return thomas_(p2.phi, p1.phi, r, datum=self._datum) 871 872 873class FrechetVincentys(FrechetRadians): 874 '''Compute the C{Frechet} distance based on the I{angular} 875 distance in C{radians} from function L{vincentys_}. 876 877 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular}, 878 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetFlatPolar} 879 and L{FrechetHaversine}. 880 881 @note: See note at function L{vincentys_}. 882 ''' 883 _wrap = False 884 885 def __init__(self, points, wrap=False, fraction=None, name=NN): 886 '''New L{FrechetVincentys} calculator/interpolator. 887 888 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 889 L{Tuple2LatLon}[] or C{other}[]). 890 @kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}). 891 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to 892 interpolate intermediate B{C{points}} or use 893 C{None}, C{0} or C{1} for no intermediate 894 B{C{points}} and no I{fractional} indices. 895 @kwarg name: Optional calculator/interpolator name (C{str}). 896 897 @raise FrechetError: Insufficient number of B{C{points}} or 898 invalid B{C{fraction}}. 899 ''' 900 FrechetRadians.__init__(self, points, fraction=fraction, name=name, 901 wrap=wrap) 902 903 if _FOR_DOCS: 904 discrete = Frechet.discrete 905 906 def distance(self, p1, p2): 907 '''Return the L{vincentys_} distance in C{radians}. 908 ''' 909 d, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap) 910 return vincentys_(p2.phi, p1.phi, d) 911 912 913def _frechet_(ni, fi, nj, fj, dF, units): # MCCABE 14 914 '''(INTERNAL) Recursive core of function L{frechet_} 915 and method C{discrete} of C{Frechet...} classes. 916 ''' 917 iFs = {} 918 919 def iF(i): # cache index, depth ints and floats 920 return iFs.setdefault(i, i) 921 922 cF = _defaultdict(dict) 923 924 def rF(i, j, r): # recursive Fréchet 925 i = iF(i) 926 j = iF(j) 927 try: 928 t = cF[i][j] 929 except KeyError: 930 r = iF(r + 1) 931 try: 932 if i > 0: 933 if j > 0: 934 t = min(rF(i - fi, j, r), 935 rF(i - fi, j - fj, r), 936 rF(i, j - fj, r)) 937 elif j < 0: 938 raise IndexError 939 else: # j == 0 940 t = rF(i - fi, 0, r) 941 elif i < 0: 942 raise IndexError 943 944 elif j > 0: # i == 0 945 t = rF(0, j - fj, r) 946 elif j < 0: # i == 0 947 raise IndexError 948 else: # i == j == 0 949 t = (-INF, i, j, r) 950 951 d = dF(i, j) 952 if d > t[0]: 953 t = (d, i, j, r) 954 except IndexError: 955 t = (INF, i, j, r) 956 cF[i][j] = t 957 return t 958 959 t = rF(ni - 1, nj - 1, 0) 960 t += (sum(map(len, cF.values())), units) 961# del cF, iFs 962 return Frechet6Tuple(*t) 963 964 965def frechet_(points1, points2, distance=None, units=NN): 966 '''Compute the I{discrete} U{Fréchet<https://WikiPedia.org/wiki/Frechet_distance>} 967 distance between two paths given as sets of points. 968 969 @arg points1: First set of points (C{LatLon}[], L{Numpy2LatLon}[], 970 L{Tuple2LatLon}[] or C{other}[]). 971 @arg points2: Second set of points (C{LatLon}[], L{Numpy2LatLon}[], 972 L{Tuple2LatLon}[] or C{other}[]). 973 @kwarg distance: Callable returning the distance between a B{C{points1}} 974 and a B{C{points2}} point (signature C{(point1, point2)}). 975 @kwarg units: Optional, the distance units (C{Unit} or C{str}). 976 977 @return: A L{Frechet6Tuple}C{(fd, fi1, fi2, r, n, units)} where C{fi1} 978 and C{fi2} are type C{int} indices into B{C{points1}} respectively 979 B{C{points2}}. 980 981 @raise FrechetError: Insufficient number of B{C{points1}} or B{C{points2}}. 982 983 @raise RecursionError: Recursion depth exceeded, see U{sys.getrecursionlimit() 984 <https://docs.Python.org/3/library/sys.html#sys.getrecursionlimit>}. 985 986 @raise TypeError: If B{C{distance}} is not a callable. 987 988 @note: Function L{frechet_} does not support I{fractional} indices for 989 intermediate B{C{points1}} and B{C{points2}}. 990 ''' 991 if not callable(distance): 992 raise _IsnotError(callable.__name__, distance=distance) 993 994 n1, ps1 = _points2(points1, closed=False, Error=FrechetError) 995 n2, ps2 = _points2(points2, closed=False, Error=FrechetError) 996 997 def dF(i1, i2): 998 return distance(ps1[i1], ps2[i2]) 999 1000 return _frechet_(n1, 1, n2, 1, dF, units) 1001 1002 1003class Frechet6Tuple(_NamedTuple): 1004 '''6-Tuple C{(fd, fi1, fi2, r, n, units)} with the I{discrete} 1005 U{Fréchet<https://WikiPedia.org/wiki/Frechet_distance>} distance 1006 C{fd}, I{fractional} indices C{fi1} and C{fi2} as C{FIx}, the 1007 recursion depth C{r}, the number of distances computed C{n} and 1008 the L{units} class or class or name of the distance C{units}. 1009 1010 If I{fractional} indices C{fi1} and C{fi2} are C{int}, the 1011 returned C{fd} is the distance between C{points1[fi1]} and 1012 C{points2[fi2]}. For C{float} indices, the distance is 1013 between an intermediate point along C{points1[int(fi1)]} and 1014 C{points1[int(fi1) + 1]} respectively an intermediate point 1015 along C{points2[int(fi2)]} and C{points2[int(fi2) + 1]}. 1016 1017 Use function L{fractional} to compute the point at a 1018 I{fractional} index. 1019 ''' 1020 _Names_ = ('fd', 'fi1', 'fi2', 'r', _n_, _units_) 1021 _Units_ = (_Pass, FIx, FIx, Number_, Number_, _Pass) 1022 1023 def toUnits(self, **Error): # PYCHOK expected 1024 '''Overloaded C{_NamedTuple.toUnits} for C{fd} units. 1025 ''' 1026 U = _xUnit(self.units, Float) # PYCHOK expected 1027 self._Units_ = (U,) + Frechet6Tuple._Units_[1:] 1028 return _NamedTuple.toUnits(self, **Error) 1029 1030# def __gt__(self, other): 1031# _xinstanceof(Frechet6Tuple, other=other) 1032# return self if self.fd > other.fd else other # PYCHOK .fd=[0] 1033# 1034# def __lt__(self, other): 1035# _xinstanceof(Frechet6Tuple, other=other) 1036# return self if self.fd < other.fd else other # PYCHOK .fd=[0] 1037 1038# **) MIT License 1039# 1040# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved. 1041# 1042# Permission is hereby granted, free of charge, to any person obtaining a 1043# copy of this software and associated documentation files (the "Software"), 1044# to deal in the Software without restriction, including without limitation 1045# the rights to use, copy, modify, merge, publish, distribute, sublicense, 1046# and/or sell copies of the Software, and to permit persons to whom the 1047# Software is furnished to do so, subject to the following conditions: 1048# 1049# The above copyright notice and this permission notice shall be included 1050# in all copies or substantial portions of the Software. 1051# 1052# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 1053# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 1054# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 1055# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 1056# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 1057# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 1058# OTHER DEALINGS IN THE SOFTWARE. 1059