Coverage for pygeodesy/frechet.py: 97%
243 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-05 15:46 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-05 15:46 -0400
2# -*- coding: utf-8 -*-
4u'''Fréchet distances.
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.
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>}.
19Typical usage is as follows. First, create a C{Frechet} calculator
20from one set of C{LatLon} points.
22C{f = FrechetXyz(points1, ...)}
24Get the I{discrete} Fréchet distance to another set of C{LatLon} points
25by
27C{t6 = f.discrete(points2)}
29Or, use function C{frechet_} with a proper C{distance} function passed
30as keyword arguments as follows
32C{t6 = frechet_(points1, points2, ..., distance=...)}.
34In both cases, the returned result C{t6} is a L{Frechet6Tuple}.
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.
40For other points, create a L{Frechet} sub-class with the appropriate
41C{distance} method overloading L{Frechet.distance} as in this example.
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)
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.
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.
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.
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.
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'''
84from pygeodesy.basics import isscalar, _xinstanceof
85from pygeodesy.constants import EPS, EPS1, INF, NINF
86from pygeodesy.datums import Datum, _WGS84
87from pygeodesy.errors import _IsnotError, PointsError
88from pygeodesy.fmath import hypot2
89from pygeodesy.formy import cosineAndoyerLambert_, cosineForsytheAndoyerLambert_, \
90 cosineLaw_, euclidean_, flatPolar_, haversine_, thomas_, \
91 vincentys_, _scale_rad
92from pygeodesy.interns import NN, _datum_, _distanceTo_, _DOT_, _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_, _xUnit, _xUnits
101from pygeodesy.unitsBase import _Str_degrees, _Str_meter, _Str_NN, \
102 _Str_radians, _Str_radians2
103from pygeodesy.utily import unrollPI
105from collections import defaultdict as _defaultdict
106from math import radians
108__all__ = _ALL_LAZY.frechet
109__version__ = '22.09.24'
112def _fraction(fraction, n):
113 f = 1 # int, no fractional indices
114 if fraction in (None, 1):
115 pass
116 elif not (isscalar(fraction) and EPS < fraction < EPS1
117 and (float(n) - fraction) < n):
118 raise FrechetError(fraction=fraction)
119 elif fraction < EPS1:
120 f = float(fraction)
121 return f
124class FrechetError(PointsError):
125 '''Fréchet issue.
126 '''
127 pass
130class Frechet(_Named):
131 '''Frechet base class, requires method L{Frechet.distance} to
132 be overloaded.
133 '''
134 _adjust = None # not applicable
135 _datum = None # not applicable
136 _f1 = 1
137 _n1 = 0
138 _ps1 = None
139 _units = _Str_NN # XXX Str to _Pass and for backward compatibility
140 _wrap = None # not applicable
142 def __init__(self, points, fraction=None, name=NN, units=NN, **wrap_adjust):
143 '''New C{Frechet...} calculator/interpolator.
145 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
146 L{Tuple2LatLon}[] or C{other}[]).
147 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
148 interpolate intermediate B{C{points}} or use
149 C{None}, C{0} or C{1} for no intermediate
150 B{C{points}} and no I{fractional} indices.
151 @kwarg name: Optional calculator/interpolator name (C{str}).
152 @kwarg units: Optional, the distance units (C{Unit} or C{str}).
153 @kwarg wrap_adjust: Optionally, C{wrap} and unroll longitudes, iff
154 applicable (C{bool}) and C{adjust} wrapped,
155 unrolled longitudinal delta by the cosine
156 of the mean latitude, iff applicable.
158 @raise FrechetError: Insufficient number of B{C{points}} or
159 invalid B{C{fraction}} or B{C{wrap}} or
160 B{C{ajust}} not applicable.
162 '''
163 self._n1, self._ps1 = self._points2(points)
164 if fraction:
165 self.fraction = fraction
166 if name:
167 self.name = name
168 if units: # and not self.units:
169 self.units = units
170 if wrap_adjust:
171 _boolkwds(self, **wrap_adjust)
173 @Property_RO
174 def adjust(self):
175 '''Get the adjust setting (C{bool} or C{None} if not applicable).
176 '''
177 return self._adjust
179 @Property_RO
180 def datum(self):
181 '''Get the datum (L{Datum} or C{None} if not applicable).
182 '''
183 return self._datum
185 def _datum_setter(self, datum):
186 '''(INTERNAL) Set the datum.
187 '''
188 d = datum or getattr(self._ps1[0], _datum_, datum)
189 if d and d != self.datum: # PYCHOK no cover
190 _xinstanceof(Datum, datum=d)
191 self._datum = d
193 def discrete(self, points, fraction=None):
194 '''Compute the C{forward, discrete Fréchet} distance.
196 @arg points: Second set of points (C{LatLon}[], L{Numpy2LatLon}[],
197 L{Tuple2LatLon}[] or C{other}[]).
198 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
199 interpolate intermediate B{C{points}} or use
200 C{None}, C{0} or C{1} for no intermediate
201 B{C{points}} and no I{fractional} indices.
203 @return: A L{Frechet6Tuple}C{(fd, fi1, fi2, r, n, units)}.
205 @raise FrechetError: Insufficient number of B{C{points}} or an
206 invalid B{C{point}} or B{C{fraction}}.
208 @raise RecursionError: Recursion depth exceeded, see U{sys.getrecursionlimit
209 <https://docs.Python.org/3/library/sys.html#sys.getrecursionlimit>}.
210 '''
211 n2, ps2 = self._points2(points)
213 f2 = _fraction(fraction, n2)
214 p2 = self.points_fraction if f2 < EPS1 else self.points_ # PYCHOK expected
216 f1 = self.fraction
217 p1 = self.points_fraction if f1 < EPS1 else self.points_ # PYCHOK expected
219 def dF(fi1, fi2):
220 return self.distance(p1(self._ps1, fi1), p2(ps2, fi2))
222 try:
223 return _frechet_(self._n1, f1, n2, f2, dF, self.units)
224 except TypeError as x:
225 t = _DOT_(self.classname, self.discrete.__name__)
226 raise FrechetError(t, cause=x)
228 def distance(self, point1, point2): # PYCHOK no cover
229 '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
230 '''
231 notOverloaded(self, point1, point2)
233 @property_doc_(''' the index fraction (C{float}).''')
234 def fraction(self):
235 '''Get the index fraction (C{float} or C{1}).
236 '''
237 return self._f1
239 @fraction.setter # PYCHOK setter!
240 def fraction(self, fraction):
241 '''Set the index fraction (C{float} in C{EPS}..L{EPS1}) to interpolate
242 intermediate B{C{points}} or use C{None}, C{0} or C{1} for no
243 intermediate B{C{points}} and no I{fractional} indices.
245 @raise FrechetError: Invalid B{C{fraction}}.
246 '''
247 self._f1 = _fraction(fraction, self._n1)
249 def point(self, point):
250 '''Convert a point for the C{.distance} method.
252 @arg point: The point to convert ((C{LatLon}, L{Numpy2LatLon},
253 L{Tuple2LatLon} or C{other}).
255 @return: The converted B{C{point}}.
256 '''
257 return point # pass thru
259 def points_(self, points, i):
260 '''Get and convert a point for the C{.distance} method.
262 @arg points: The orignal B{C{points}} to convert ((C{LatLon}[],
263 L{Numpy2LatLon}[], L{Tuple2LatLon}[] or C{other}[]).
264 @arg i: The B{C{points}} index (C{int}).
266 @return: The converted B{C{points[i]}}.
267 '''
268 return self.point(points[i])
270 def _points2(self, points):
271 '''(INTERNAL) Check a set of points.
272 '''
273 return _points2(points, closed=False, Error=FrechetError)
275 def points_fraction(self, points, fi):
276 '''Get and convert a I{fractional} point for the C{.distance} method.
278 @arg points: The orignal B{C{points}} to convert ((C{LatLon}[],
279 L{Numpy2LatLon}[], L{Tuple2LatLon}[] or C{other}[]).
280 @arg fi: The I{fractional} index in B{C{points}} (C{float} or C{int}).
282 @return: The interpolated, converted, intermediate B{C{points[fi]}}.
283 '''
284 return self.point(_fractional(points, fi, None)) # wrap=None?
286 @property_doc_(''' the distance units (C{Unit} or C{str}).''')
287 def units(self):
288 '''Get the distance units (C{Unit} or C{str}).
289 '''
290 return self._units
292 @units.setter # PYCHOK setter!
293 def units(self, units):
294 '''Set the distance units (C{Unit} or C{str}).
296 @raise TypeError: Invalid B{C{units}}.
297 '''
298 self._units = _xUnits(units, Base=Float)
300 @Property_RO
301 def wrap(self):
302 '''Get the wrap setting (C{bool} or C{None} if not applicable).
303 '''
304 return self._wrap
307class FrechetDegrees(Frechet):
308 '''L{Frechet} base class for distances in C{degrees} from
309 C{LatLon} points in C{degrees}.
310 '''
311 _units = _Str_degrees
313 if _FOR_DOCS:
314 __init__ = Frechet.__init__
315 discrete = Frechet.discrete
318class FrechetRadians(Frechet):
319 '''L{Frechet} base class for distances in C{radians} or C{radians
320 squared} from C{LatLon} points converted from C{degrees} to
321 C{radians}.
322 '''
323 _units = _Str_radians
325 if _FOR_DOCS:
326 __init__ = Frechet.__init__
327 discrete = Frechet.discrete
329 def point(self, point):
330 '''Convert C{(lat, lon)} point in degrees to C{(a, b)}
331 in radians.
333 @return: An L{PhiLam2Tuple}C{(phi, lam)}.
334 '''
335 try:
336 return point.philam
337 except AttributeError:
338 return PhiLam2Tuple(radians(point.lat), radians(point.lon))
341class FrechetCosineAndoyerLambert(FrechetRadians):
342 '''Compute the C{Frechet} distance based on the I{angular} distance
343 in C{radians} from function L{pygeodesy.cosineAndoyerLambert_}.
345 @see: L{FrechetCosineForsytheAndoyerLambert}, L{FrechetDistanceTo},
346 L{FrechetExact}, L{FrechetFlatLocal}, L{FrechetHubeny},
347 L{FrechetThomas} and L{FrechetKarney}.
348 '''
349 _datum = _WGS84
350 _wrap = False
352 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):
353 '''New L{FrechetCosineAndoyerLambert} calculator/interpolator.
355 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
356 L{Tuple2LatLon}[] or C{other}[]).
357 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
358 and first B{C{points}}' datum (L{Datum}).
359 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool})
360 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
361 interpolate intermediate B{C{points}} or use C{0},
362 C{1} or C{None} to avoid intermediate B{C{points}}
363 and L{pygeodesy.fractional} indices.
364 @kwarg name: Optional calculator/interpolator name (C{str}).
366 @raise FrechetError: Insufficient number of B{C{points}} or
367 invalid B{C{fraction}}.
369 @raise TypeError: Invalid B{C{datum}}.
370 '''
371 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
372 wrap=wrap)
373 self._datum_setter(datum)
375 if _FOR_DOCS:
376 discrete = Frechet.discrete
378 def distance(self, p1, p2):
379 '''Return the L{pygeodesy.cosineAndoyerLambert_} distance in C{radians}.
380 '''
381 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
382 return cosineAndoyerLambert_(p2.phi, p1.phi, r, datum=self._datum)
385class FrechetCosineForsytheAndoyerLambert(FrechetRadians):
386 '''Compute the C{Frechet} distance based on the I{angular} distance
387 in C{radians} from function L{pygeodesy.cosineForsytheAndoyerLambert_}.
389 @see: L{FrechetCosineAndoyerLambert}, L{FrechetDistanceTo},
390 L{FrechetExact}, L{FrechetFlatLocal}, L{FrechetHubeny},
391 L{FrechetThomas} and L{FrechetKarney}.
392 '''
393 _datum = _WGS84
394 _wrap = False
396 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):
397 '''New L{FrechetCosineForsytheAndoyerLambert} calculator/interpolator.
399 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
400 L{Tuple2LatLon}[] or C{other}[]).
401 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
402 and first B{C{points}}' datum (L{Datum}).
403 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool})
404 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
405 interpolate intermediate B{C{points}} or use C{0},
406 C{1} or C{None} to avoid intermediate B{C{points}}
407 and L{pygeodesy.fractional} indices.
408 @kwarg name: Optional calculator/interpolator name (C{str}).
410 @raise FrechetError: Insufficient number of B{C{points}} or
411 invalid B{C{fraction}}.
413 @raise TypeError: Invalid B{C{datum}}.
414 '''
415 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
416 wrap=wrap)
417 self._datum_setter(datum)
419 if _FOR_DOCS:
420 discrete = Frechet.discrete
422 def distance(self, p1, p2):
423 '''Return the L{pygeodesy.cosineForsytheAndoyerLambert_} distance in C{radians}.
424 '''
425 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
426 return cosineForsytheAndoyerLambert_(p2.phi, p1.phi, r, datum=self._datum)
429class FrechetCosineLaw(FrechetRadians):
430 '''Compute the C{Frechet} distance based on the I{angular} distance
431 in C{radians} from function L{pygeodesy.cosineLaw_}.
433 @see: L{FrechetEquirectangular}, L{FrechetEuclidean},
434 L{FrechetExact}, L{FrechetFlatPolar}, L{FrechetHaversine}
435 and L{FrechetVincentys}.
437 @note: See note at function L{pygeodesy.vincentys_}.
438 '''
439 _wrap = False
441 def __init__(self, points, wrap=False, fraction=None, name=NN):
442 '''New L{FrechetCosineLaw} calculator/interpolator.
444 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
445 L{Tuple2LatLon}[] or C{other}[]).
446 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
447 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
448 interpolate intermediate B{C{points}} or use C{0},
449 C{1} or C{None} to avoid intermediate B{C{points}}
450 and L{pygeodesy.fractional} indices.
451 @kwarg name: Optional calculator/interpolator name (C{str}).
453 @raise FrechetError: Insufficient number of B{C{points}} or
454 invalid B{C{fraction}}.
455 '''
456 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
457 wrap=wrap)
459 if _FOR_DOCS:
460 discrete = Frechet.discrete
462 def distance(self, p1, p2):
463 '''Return the L{pygeodesy.cosineLaw_} distance in C{radians}.
464 '''
465 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
466 return cosineLaw_(p2.phi, p1.phi, r)
469class FrechetDistanceTo(Frechet):
470 '''Compute the C{Frechet} distance based on the distance from the
471 points' C{LatLon.distanceTo} method, conventionally in C{meter}.
473 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert},
474 L{FrechetExact}, L{FrechetFlatLocal}, L{FrechetHubeny}, L{FrechetThomas}
475 and L{FrechetKarney}.
476 '''
477 _distanceTo_kwds = {}
478 _units = _Str_meter
480 def __init__(self, points, fraction=None, name=NN, **distanceTo_kwds):
481 '''New L{FrechetDistanceTo} calculator/interpolator.
483 @arg points: First set of points (C{LatLon}[]).
484 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
485 interpolate intermediate B{C{points}} or use C{0},
486 C{1} or C{None} to avoid intermediate B{C{points}}
487 and L{pygeodesy.fractional} indices.
488 @kwarg name: Optional calculator/interpolator name (C{str}).
489 @kwarg distanceTo_kwds: Optional keyword arguments for the
490 B{C{points}}' C{LatLon.distanceTo}
491 method.
493 @raise FrechetError: Insufficient number of B{C{points}} or an
494 invalid B{C{point}} or B{C{fraction}}.
496 @raise ImportError: Package U{geographiclib
497 <https://PyPI.org/project/geographiclib>} missing
498 iff B{C{points}} are L{ellipsoidalKarney.LatLon}s.
500 @note: All B{C{points}} I{must} be instances of the same
501 ellipsoidal or spherical C{LatLon} class.
502 '''
503 FrechetRadians.__init__(self, points, fraction=fraction, name=name)
504 if distanceTo_kwds:
505 self._distanceTo_kwds = distanceTo_kwds
507 if _FOR_DOCS:
508 discrete = Frechet.discrete
510 def distance(self, p1, p2):
511 '''Return the distance in C{meter}.
512 '''
513 return p1.distanceTo(p2, **self._distanceTo_kwds)
515 def _points2(self, points):
516 '''(INTERNAL) Check a set of points.
517 '''
518 np, ps = Frechet._points2(self, points)
519 for i, p in enumerate(ps):
520 if not callable(getattr(p, _distanceTo_, None)):
521 i = Fmt.SQUARE(_points_, i)
522 raise FrechetError(i, p, txt=_distanceTo_)
523 return np, ps
526class FrechetEquirectangular(FrechetRadians):
527 '''Compute the C{Frechet} distance based on the I{equirectangular}
528 distance in C{radians squared} like function L{pygeodesy.equirectangular_}.
530 @see: L{FrechetCosineLaw}, L{FrechetEuclidean}, L{FrechetExact},
531 L{FrechetFlatPolar}, L{FrechetHaversine} and L{FrechetVincentys}.
532 '''
533 _adjust = True
534 _units = _Str_radians2
535 _wrap = False
537 def __init__(self, points, adjust=True, wrap=False, fraction=None, name=NN):
538 '''New L{FrechetEquirectangular} calculator/interpolator.
540 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
541 L{Tuple2LatLon}[] or C{other}[]).
542 @kwarg adjust: Adjust the wrapped, unrolled longitudinal
543 delta by the cosine of the mean latitude (C{bool}).
544 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
545 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
546 interpolate intermediate B{C{points}} or use C{0},
547 C{1} or C{None} to avoid intermediate B{C{points}}
548 and L{pygeodesy.fractional} indices.
549 @kwarg name: Optional calculator/interpolator name (C{str}).
551 @raise FrechetError: Insufficient number of B{C{points}} or
552 invalid B{C{adjust}} or B{C{seed}}.
553 '''
554 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
555 adjust=adjust, wrap=wrap)
557 if _FOR_DOCS:
558 discrete = Frechet.discrete
560 def distance(self, p1, p2):
561 '''Return the L{pygeodesy.equirectangular_} distance in C{radians squared}.
562 '''
563 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
564 if self._adjust:
565 r *= _scale_rad(p1.phi, p2.phi)
566 return hypot2(r, p2.phi - p1.phi) # like equirectangular_ d2
569class FrechetEuclidean(FrechetRadians):
570 '''Compute the C{Frechet} distance based on the I{Euclidean}
571 distance in C{radians} from function L{pygeodesy.euclidean_}.
573 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular},
574 L{FrechetExact}, L{FrechetFlatPolar}, L{FrechetHaversine}
575 and L{FrechetVincentys}.
576 '''
577 _adjust = True
578 _wrap = True # fixed
580 def __init__(self, points, adjust=True, fraction=None, name=NN):
581 '''New L{FrechetEuclidean} calculator/interpolator.
583 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
584 L{Tuple2LatLon}[] or C{other}[]).
585 @kwarg adjust: Adjust the wrapped, unrolled longitudinal
586 delta by the cosine of the mean latitude (C{bool}).
587 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
588 interpolate intermediate B{C{points}} or use C{0},
589 C{1} or C{None} to avoid intermediate B{C{points}}
590 and L{pygeodesy.fractional} indices.
591 @kwarg name: Optional calculator/interpolator name (C{str}).
593 @raise FrechetError: Insufficient number of B{C{points}} or
594 invalid B{C{fraction}}.
595 '''
596 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
597 adjust=adjust)
599 if _FOR_DOCS:
600 discrete = Frechet.discrete
602 def distance(self, p1, p2):
603 '''Return the L{pygeodesy.euclidean_} distance in C{radians}.
604 '''
605 return euclidean_(p2.phi, p1.phi, p2.lam - p1.lam, adjust=self._adjust)
608class FrechetExact(FrechetDegrees):
609 '''Compute the C{Frechet} distance based on the I{angular}
610 distance in C{degrees} from method L{GeodesicExact}C{.Inverse}.
612 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert},
613 L{FrechetDistanceTo}, L{FrechetFlatLocal}, L{FrechetHubeny},
614 L{FrechetKarney} and L{FrechetThomas}.
615 '''
616 _datum = _WGS84
617 _Inverse1 = None
618 _wrap = False
620 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):
621 '''New L{FrechetExact} calculator/interpolator.
623 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
624 L{Tuple2LatLon}[] or C{other}[]).
625 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
626 and first B{C{points}}' datum (L{Datum}).
627 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool})
628 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
629 interpolate intermediate B{C{points}} or use C{0},
630 C{1} or C{None} to avoid intermediate B{C{points}}
631 and L{pygeodesy.fractional} indices.
632 @kwarg name: Optional calculator/interpolator name (C{str}).
634 @raise FrechetError: Insufficient number of B{C{points}} or
635 invalid B{C{fraction}}.
637 @raise TypeError: Invalid B{C{datum}}.
638 '''
639 FrechetDegrees.__init__(self, points, fraction=fraction, name=name,
640 wrap=wrap)
641 self._datum_setter(datum)
642 self._Inverse1 = self.datum.ellipsoid.geodesicx.Inverse1 # note -x
644 if _FOR_DOCS:
645 discrete = Frechet.discrete
647 def distance(self, p1, p2):
648 '''Return the non-negative I{angular} distance in C{degrees}.
649 '''
650 return self._Inverse1(p1.lat, p1.lon, p2.lat, p2.lon, wrap=self._wrap)
653class FrechetFlatLocal(FrechetRadians):
654 '''Compute the C{Frechet} distance based on the I{angular} distance in
655 C{radians squared} like function L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_}.
657 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert},
658 L{FrechetDistanceTo}, L{FrechetExact}, L{FrechetHubeny},
659 L{FrechetKarney} and L{FrechetThomas}.
660 '''
661 _datum = _WGS84
662 _units = _Str_radians2
663 _wrap = False
665 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):
666 '''New L{FrechetFlatLocal}/L{FrechetHubeny} calculator/interpolator.
668 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
669 L{Tuple2LatLon}[] or C{other}[]).
670 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
671 and first B{C{points}}' datum (L{Datum}).
672 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool})
673 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
674 interpolate intermediate B{C{points}} or use C{0},
675 C{1} or C{None} to avoid intermediate B{C{points}}
676 and L{pygeodesy.fractional} indices.
677 @kwarg name: Optional calculator/interpolator name (C{str}).
679 @raise FrechetError: Insufficient number of B{C{points}} or
680 invalid B{C{fraction}}.
682 @raise TypeError: Invalid B{C{datum}}.
683 '''
684 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
685 wrap=wrap)
686 self._datum_setter(datum)
688 if _FOR_DOCS:
689 discrete = Frechet.discrete
691 def distance(self, p1, p2):
692 '''Return the L{pygeodesy.flatLocal_}/L{pygeodesy.hubeny_} distance
693 in C{radians squared}.
694 '''
695 d, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
696 return self._hubeny_2(p2.phi, p1.phi, d)
698 @Property_RO
699 def _hubeny_2(self):
700 '''(INTERNAL) Get and cache the C{.datum.ellipsoid._hubeny_2} method.
701 '''
702 return self.datum.ellipsoid._hubeny_2
705class FrechetFlatPolar(FrechetRadians):
706 '''Compute the C{Frechet} distance based on the I{angular} distance
707 in C{radians} from function L{flatPolar_}.
709 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular},
710 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetHaversine}
711 and L{FrechetVincentys}.
712 '''
713 _wrap = False
715 def __init__(self, points, wrap=False, fraction=None, name=NN):
716 '''New L{FrechetFlatPolar} calculator/interpolator.
718 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
719 L{Tuple2LatLon}[] or C{other}[]).
720 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
721 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
722 interpolate intermediate B{C{points}} or use C{0},
723 C{1} or C{None} to avoid intermediate B{C{points}}
724 and L{pygeodesy.fractional} indices.
725 @kwarg name: Optional calculator/interpolator name (C{str}).
727 @raise FrechetError: Insufficient number of B{C{points}} or
728 invalid B{C{fraction}}.
729 '''
730 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
731 wrap=wrap)
733 if _FOR_DOCS:
734 discrete = Frechet.discrete
736 def distance(self, p1, p2):
737 '''Return the L{flatPolar_} distance in C{radians}.
738 '''
739 d, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
740 return flatPolar_(p2.phi, p1.phi, d)
743class FrechetHaversine(FrechetRadians):
744 '''Compute the C{Frechet} distance based on the I{angular}
745 distance in C{radians} from function L{pygeodesy.haversine_}.
747 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular},
748 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetFlatPolar}
749 and L{FrechetVincentys}.
751 @note: See note at function L{pygeodesy.vincentys_}.
752 '''
753 _wrap = False
755 def __init__(self, points, wrap=False, fraction=None, name=NN):
756 '''New L{FrechetHaversine} calculator/interpolator.
758 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
759 L{Tuple2LatLon}[] or C{other}[]).
760 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
761 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
762 interpolate intermediate B{C{points}} or use C{0},
763 C{1} or C{None} to avoid intermediate B{C{points}}
764 and L{pygeodesy.fractional} indices.
765 @kwarg name: Optional calculator/interpolator name (C{str}).
767 @raise FrechetError: Insufficient number of B{C{points}} or
768 invalid B{C{fraction}}.
769 '''
770 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
771 wrap=wrap)
773 if _FOR_DOCS:
774 discrete = Frechet.discrete
776 def distance(self, p1, p2):
777 '''Return the L{pygeodesy.haversine_} distance in C{radians}.
778 '''
779 d, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
780 return haversine_(p2.phi, p1.phi, d)
783class FrechetHubeny(FrechetFlatLocal): # for Karl Hubeny
784 if _FOR_DOCS:
785 __doc__ = FrechetFlatLocal.__doc__
786 __init__ = FrechetFlatLocal.__init__
787 discrete = FrechetFlatLocal.discrete
788 distance = FrechetFlatLocal.discrete
791class FrechetKarney(FrechetExact):
792 '''Compute the C{Frechet} distance based on the I{angular}
793 distance in C{degrees} from I{Karney}'s U{geographiclib
794 <https://PyPI.org/project/geographiclib>} U{Geodesic
795 <https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}
796 Inverse method.
798 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert},
799 L{FrechetDistanceTo}, L{FrechetExact}, L{FrechetFlatLocal},
800 L{FrechetHubeny} and L{FrechetThomas}.
801 '''
802 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):
803 '''New L{FrechetKarney} calculator/interpolator.
805 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
806 L{Tuple2LatLon}[] or C{other}[]).
807 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
808 and first B{C{points}}' datum (L{Datum}).
809 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool})
810 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
811 interpolate intermediate B{C{points}} or use C{0},
812 C{1} or C{None} to avoid intermediate B{C{points}}
813 and L{pygeodesy.fractional} indices.
814 @kwarg name: Optional calculator/interpolator name (C{str}).
816 @raise FrechetError: Insufficient number of B{C{points}} or
817 invalid B{C{fraction}}.
819 @raise ImportError: Package U{geographiclib
820 <https://PyPI.org/project/geographiclib>} missing.
822 @raise TypeError: Invalid B{C{datum}}.
823 '''
824 FrechetDegrees.__init__(self, points, fraction=fraction, name=name,
825 wrap=wrap)
826 self._datum_setter(datum)
827 self._Inverse1 = self.datum.ellipsoid.geodesic.Inverse1
830class FrechetThomas(FrechetRadians):
831 '''Compute the C{Frechet} distance based on the I{angular} distance
832 in C{radians} from function L{pygeodesy.thomas_}.
834 @see: L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert},
835 L{FrechetDistanceTo}, L{FrechetExact}, L{FrechetFlatLocal},
836 L{FrechetHubeny} and L{FrechetKarney}.
837 '''
838 _datum = _WGS84
839 _wrap = False
841 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):
842 '''New L{FrechetThomas} calculator/interpolator.
844 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
845 L{Tuple2LatLon}[] or C{other}[]).
846 @kwarg datum: Optional datum overriding the default C{Datums.WGS84}
847 and first B{C{points}}' datum (L{Datum}).
848 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool})
849 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
850 interpolate intermediate B{C{points}} or use C{0},
851 C{1} or C{None} to avoid intermediate B{C{points}}
852 and L{pygeodesy.fractional} indices.
853 @kwarg name: Optional calculator/interpolator name (C{str}).
855 @raise FrechetError: Insufficient number of B{C{points}} or
856 invalid B{C{fraction}}.
858 @raise TypeError: Invalid B{C{datum}}.
859 '''
860 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
861 wrap=wrap)
862 self._datum_setter(datum)
864 if _FOR_DOCS:
865 discrete = Frechet.discrete
867 def distance(self, p1, p2):
868 '''Return the L{pygeodesy.thomas_} distance in C{radians}.
869 '''
870 r, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
871 return thomas_(p2.phi, p1.phi, r, datum=self._datum)
874class FrechetVincentys(FrechetRadians):
875 '''Compute the C{Frechet} distance based on the I{angular}
876 distance in C{radians} from function L{pygeodesy.vincentys_}.
878 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular},
879 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetFlatPolar}
880 and L{FrechetHaversine}.
882 @note: See note at function L{pygeodesy.vincentys_}.
883 '''
884 _wrap = False
886 def __init__(self, points, wrap=False, fraction=None, name=NN):
887 '''New L{FrechetVincentys} calculator/interpolator.
889 @arg points: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
890 L{Tuple2LatLon}[] or C{other}[]).
891 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
892 @kwarg fraction: Index fraction (C{float} in L{EPS}..L{EPS1}) to
893 interpolate intermediate B{C{points}} or use C{0},
894 C{1} or C{None} to avoid intermediate B{C{points}}
895 and L{pygeodesy.fractional} indices.
896 @kwarg name: Optional calculator/interpolator name (C{str}).
898 @raise FrechetError: Insufficient number of B{C{points}} or
899 invalid B{C{fraction}}.
900 '''
901 FrechetRadians.__init__(self, points, fraction=fraction, name=name,
902 wrap=wrap)
904 if _FOR_DOCS:
905 discrete = Frechet.discrete
907 def distance(self, p1, p2):
908 '''Return the L{pygeodesy.vincentys_} distance in C{radians}.
909 '''
910 d, _ = unrollPI(p1.lam, p2.lam, wrap=self._wrap)
911 return vincentys_(p2.phi, p1.phi, d)
914def _frechet_(ni, fi, nj, fj, dF, units): # MCCABE 14
915 '''(INTERNAL) Recursive core of function L{frechet_}
916 and method C{discrete} of C{Frechet...} classes.
917 '''
918 iFs = {}
920 def iF(i): # cache index, depth ints and floats
921 return iFs.setdefault(i, i)
923 cF = _defaultdict(dict)
925 def rF(i, j, r): # recursive Fréchet
926 i = iF(i)
927 j = iF(j)
928 try:
929 t = cF[i][j]
930 except KeyError:
931 r = iF(r + 1)
932 try:
933 if i > 0:
934 if j > 0:
935 t = min(rF(i - fi, j, r),
936 rF(i - fi, j - fj, r),
937 rF(i, j - fj, r))
938 elif j < 0:
939 raise IndexError
940 else: # j == 0
941 t = rF(i - fi, 0, r)
942 elif i < 0:
943 raise IndexError
945 elif j > 0: # i == 0
946 t = rF(0, j - fj, r)
947 elif j < 0: # i == 0
948 raise IndexError
949 else: # i == j == 0
950 t = (NINF, i, j, r)
952 d = dF(i, j)
953 if d > t[0]:
954 t = (d, i, j, r)
955 except IndexError:
956 t = (INF, i, j, r)
957 cF[i][j] = t
958 return t
960 t = rF(ni - 1, nj - 1, 0)
961 t += (sum(map(len, cF.values())), units)
962# del cF, iFs
963 return Frechet6Tuple(t) # *t
966def frechet_(points1, points2, distance=None, units=NN):
967 '''Compute the I{discrete} U{Fréchet<https://WikiPedia.org/wiki/Frechet_distance>}
968 distance between two paths given as sets of points.
970 @arg points1: First set of points (C{LatLon}[], L{Numpy2LatLon}[],
971 L{Tuple2LatLon}[] or C{other}[]).
972 @arg points2: Second set of points (C{LatLon}[], L{Numpy2LatLon}[],
973 L{Tuple2LatLon}[] or C{other}[]).
974 @kwarg distance: Callable returning the distance between a B{C{points1}}
975 and a B{C{points2}} point (signature C{(point1, point2)}).
976 @kwarg units: Optional, the distance units (C{Unit} or C{str}).
978 @return: A L{Frechet6Tuple}C{(fd, fi1, fi2, r, n, units)} where C{fi1}
979 and C{fi2} are type C{int} indices into B{C{points1}} respectively
980 B{C{points2}}.
982 @raise FrechetError: Insufficient number of B{C{points1}} or B{C{points2}}.
984 @raise RecursionError: Recursion depth exceeded, see U{sys.getrecursionlimit()
985 <https://docs.Python.org/3/library/sys.html#sys.getrecursionlimit>}.
987 @raise TypeError: If B{C{distance}} is not a callable.
989 @note: Function L{frechet_} does not support I{fractional} indices for
990 intermediate B{C{points1}} and B{C{points2}}.
991 '''
992 if not callable(distance):
993 raise _IsnotError(callable.__name__, distance=distance)
995 n1, ps1 = _points2(points1, closed=False, Error=FrechetError)
996 n2, ps2 = _points2(points2, closed=False, Error=FrechetError)
998 def dF(i1, i2):
999 return distance(ps1[i1], ps2[i2])
1001 return _frechet_(n1, 1, n2, 1, dF, units)
1004class Frechet6Tuple(_NamedTuple):
1005 '''6-Tuple C{(fd, fi1, fi2, r, n, units)} with the I{discrete}
1006 U{Fréchet<https://WikiPedia.org/wiki/Frechet_distance>} distance
1007 C{fd}, I{fractional} indices C{fi1} and C{fi2} as C{FIx}, the
1008 recursion depth C{r}, the number of distances computed C{n} and
1009 the L{units} class or class or name of the distance C{units}.
1011 If I{fractional} indices C{fi1} and C{fi2} are C{int}, the
1012 returned C{fd} is the distance between C{points1[fi1]} and
1013 C{points2[fi2]}. For C{float} indices, the distance is
1014 between an intermediate point along C{points1[int(fi1)]} and
1015 C{points1[int(fi1) + 1]} respectively an intermediate point
1016 along C{points2[int(fi2)]} and C{points2[int(fi2) + 1]}.
1018 Use function L{fractional} to compute the point at a
1019 I{fractional} index.
1020 '''
1021 _Names_ = ('fd', 'fi1', 'fi2', 'r', _n_, _units_)
1022 _Units_ = (_Pass, FIx, FIx, Number_, Number_, _Pass)
1024 def toUnits(self, **Error): # PYCHOK expected
1025 '''Overloaded C{_NamedTuple.toUnits} for C{fd} units.
1026 '''
1027 U = _xUnit(self.units, Float) # PYCHOK expected
1028 self._Units_ = (U,) + Frechet6Tuple._Units_[1:]
1029 return _NamedTuple.toUnits(self, **Error)
1031# def __gt__(self, other):
1032# _xinstanceof(Frechet6Tuple, other=other)
1033# return self if self.fd > other.fd else other # PYCHOK .fd=[0]
1034#
1035# def __lt__(self, other):
1036# _xinstanceof(Frechet6Tuple, other=other)
1037# return self if self.fd < other.fd else other # PYCHOK .fd=[0]
1039# **) MIT License
1040#
1041# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1042#
1043# Permission is hereby granted, free of charge, to any person obtaining a
1044# copy of this software and associated documentation files (the "Software"),
1045# to deal in the Software without restriction, including without limitation
1046# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1047# and/or sell copies of the Software, and to permit persons to whom the
1048# Software is furnished to do so, subject to the following conditions:
1049#
1050# The above copyright notice and this permission notice shall be included
1051# in all copies or substantial portions of the Software.
1052#
1053# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1054# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1055# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1056# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1057# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1058# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1059# OTHER DEALINGS IN THE SOFTWARE.