1

2# -*- coding: utf-8 -*-

3

4 u'''Fréchet distances.

5

6 Classes L{Frechet}, L{FrechetDegrees}, L{FrechetRadians},

7 L{FrechetCosineAndoyerLambert}, L{FrechetCosineForsytheAndoyerLambert},

8 L{FrechetCosineLaw}, L{FrechetDistanceTo}< L{FrechetEquirectangular},

9 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetFlatLocal}, L{FrechetFlatPolar},

10 L{FrechetHaversine}, L{FrechetHubeny}, L{FrechetKarney}, L{FrechetThomas}

11 and L{FrechetVincentys} to compute I{discrete} U{Fréchet

12<https://WikiPedia.org/wiki/Frechet_distance>} distances between two sets

13 of C{LatLon}, C{NumPy}, C{tuples} or other types of points.

14

15 Only L{FrechetDistanceTo} -iff used with L{ellipsoidalKarney.LatLon}

16 points- and L{FrechetKarney} requires installation of I{Charles Karney}'s

17 U{geographiclib<https://PyPI.org/project/geographiclib>}.

18

19 Typical usage is as follows. First, create a C{Frechet} calculator

20 from one set of C{LatLon} points.

21

22 C{f = FrechetXyz(points1, ...)}

23

24 Get the I{discrete} Fréchet distance to another set of C{LatLon} points

25 by

26

27 C{t6 = f.discrete(points2)}

28

29 Or, use function C{frechet_} with a proper C{distance} function passed

30 as keyword arguments as follows

31

32 C{t6 = frechet_(points1, points2, ..., distance=...)}.

33

34 In both cases, the returned result C{t6} is a L{Frechet6Tuple}.

35

36 For C{(lat, lon, ...)} points in a C{NumPy} array or plain C{tuples},

37 wrap the points in a L{Numpy2LatLon} respectively L{Tuple2LatLon}

38 instance, more details in the documentation thereof.

39

40 For other points, create a L{Frechet} sub-class with the appropriate

41 C{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

54 Transcribed from the original U{Computing Discrete Fréchet Distance

55<https://www.kr.TUWien.ac.AT/staff/eiter/et-archive/cdtr9464.pdf>} by

56 Eiter, T. and Mannila, H., 1994, April 25, Technical Report CD-TR 94/64,

57 Information Systems Department/Christian Doppler Laboratory for Expert

58 Systems, Technical University Vienna, Austria.

59

60 This L{Frechet.discrete} implementation optionally generates intermediate

61 points for each point set separately. For example, using keyword argument

62 C{fraction=0.5} adds one additional point halfway between each pair of

63 points. Or using C{fraction=0.1} interpolates nine additional points

64 between each points pair.

65

66 The L{Frechet6Tuple} attributes C{fi1} and/or C{fi2} will be I{fractional}

67 indices of type C{float} if keyword argument C{fraction} is used. Otherwise,

68 C{fi1} and/or C{fi2} are simply type C{int} indices into the respective

69 points set.

70

71 For example, C{fractional} index value 2.5 means an intermediate point

72 halfway between points[2] and points[3]. Use function L{fractional}

73 to obtain the intermediate point for a I{fractional} index in the

74 corresponding set of points.

75

76 The 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>}.

79 It is a measure of similarity between curves that takes into account the

80 location and ordering of the points. Therefore, it is often a better metric

81 than the well-known C{Hausdorff} distance, see the L{hausdorff} module.

82'''

83

84 from pygeodesy.basics import isscalar, _xinstanceof

85 from pygeodesy.constants import EPS, EPS1, INF, NINF

86 from pygeodesy.datums import Datum, _WGS84

87 from pygeodesy.errors import _IsnotError, PointsError

88 from pygeodesy.fmath import hypot2

89 from pygeodesy.formy import cosineAndoyerLambert_, cosineForsytheAndoyerLambert_, \

90 cosineLaw_, euclidean_, flatPolar_, haversine_, thomas_, \

91 vincentys_, _scale_rad

92 from pygeodesy.interns import NN, _datum_, _distanceTo_, _DOT_, _n_, _points_, _units_

93 from pygeodesy.iters import points2 as _points2

94 from pygeodesy.lazily import _ALL_LAZY, _FOR_DOCS

95 from pygeodesy.named import _Named, _NamedTuple, notOverloaded, _Pass

96 from pygeodesy.namedTuples import PhiLam2Tuple

97 from pygeodesy.points import _fractional

98 from pygeodesy.props import Property_RO, property_doc_

99 from pygeodesy.streprs import _boolkwds, Fmt

100 from pygeodesy.units import FIx, Float, Number_, _xUnit, _xUnits

101 from pygeodesy.unitsBase import _Str_degrees, _Str_meter, _Str_NN, \

102 _Str_radians, _Str_radians2

103 from pygeodesy.utily import unrollPI

104

105 from collections import defaultdict as _defaultdict

106 from math import radians

107

108__all__ = _ALL_LAZY.frechet

109__version__ = '22.09.24'

110

111

112 def _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

122

123

124 class FrechetError(PointsError):

125 '''Fréchet issue.

126 '''

127 pass

128

129

130 class 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

141

142 def __init__(self, points, fraction=None, name=NN, units=NN, **wrap_adjust):

143 '''New C{Frechet...} calculator/interpolator.

144

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.

157

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.

161

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)

172

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

178

179 @Property_RO

180 def datum(self):

181 '''Get the datum (L{Datum} or C{None} if not applicable).

182 '''

183 return self._datum

184

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

192

193 def discrete(self, points, fraction=None):

194 '''Compute the C{forward, discrete Fréchet} distance.

195

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.

202

203 @return: A L{Frechet6Tuple}C{(fd, fi1, fi2, r, n, units)}.

204

205 @raise FrechetError: Insufficient number of B{C{points}} or an

206 invalid B{C{point}} or B{C{fraction}}.

207

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)

212

213 f2 = _fraction(fraction, n2)

214 p2 = self.points_fraction if f2 < EPS1 else self.points_ # PYCHOK expected

215

216 f1 = self.fraction

217 p1 = self.points_fraction if f1 < EPS1 else self.points_ # PYCHOK expected

218

219 def dF(fi1, fi2):

220 return self.distance(p1(self._ps1, fi1), p2(ps2, fi2))

221

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)

227

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)

232

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

238

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.

244

245 @raise FrechetError: Invalid B{C{fraction}}.

246 '''

247 self._f1 = _fraction(fraction, self._n1)

248

249 def point(self, point):

250 '''Convert a point for the C{.distance} method.

251

252 @arg point: The point to convert ((C{LatLon}, L{Numpy2LatLon},

253 L{Tuple2LatLon} or C{other}).

254

255 @return: The converted B{C{point}}.

256 '''

257 return point # pass thru

258

259 def points_(self, points, i):

260 '''Get and convert a point for the C{.distance} method.

261

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}).

265

266 @return: The converted B{C{points[i]}}.

267 '''

268 return self.point(points[i])

269

270 def _points2(self, points):

271 '''(INTERNAL) Check a set of points.

272 '''

273 return _points2(points, closed=False, Error=FrechetError)

274

275 def points_fraction(self, points, fi):

276 '''Get and convert a I{fractional} point for the C{.distance} method.

277

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}).

281

282 @return: The interpolated, converted, intermediate B{C{points[fi]}}.

283 '''

284 return self.point(_fractional(points, fi, None)) # wrap=None?

285

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

291

292 @units.setter # PYCHOK setter!

293 def units(self, units):

294 '''Set the distance units (C{Unit} or C{str}).

295

296 @raise TypeError: Invalid B{C{units}}.

297 '''

298 self._units = _xUnits(units, Base=Float)

299

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

305

306

307 class 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

312

313 if _FOR_DOCS:

314 __init__ = Frechet.__init__

315 discrete = Frechet.discrete

316

317

318 class 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

324

325 if _FOR_DOCS:

326 __init__ = Frechet.__init__

327 discrete = Frechet.discrete

328

329 def point(self, point):

330 '''Convert C{(lat, lon)} point in degrees to C{(a, b)}

331 in radians.

332

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))

339

340

341 class FrechetCosineAndoyerLambert(FrechetRadians):

342 '''Compute the C{Frechet} distance based on the I{angular} distance

343 in C{radians} from function L{pygeodesy.cosineAndoyerLambert_}.

344

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

351

352 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):

353 '''New L{FrechetCosineAndoyerLambert} calculator/interpolator.

354

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}).

365

366 @raise FrechetError: Insufficient number of B{C{points}} or

367 invalid B{C{fraction}}.

368

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)

374

375 if _FOR_DOCS:

376 discrete = Frechet.discrete

377

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)

383

384

385 class FrechetCosineForsytheAndoyerLambert(FrechetRadians):

386 '''Compute the C{Frechet} distance based on the I{angular} distance

387 in C{radians} from function L{pygeodesy.cosineForsytheAndoyerLambert_}.

388

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

395

396 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):

397 '''New L{FrechetCosineForsytheAndoyerLambert} calculator/interpolator.

398

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}).

409

410 @raise FrechetError: Insufficient number of B{C{points}} or

411 invalid B{C{fraction}}.

412

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)

418

419 if _FOR_DOCS:

420 discrete = Frechet.discrete

421

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)

427

428

429 class FrechetCosineLaw(FrechetRadians):

430 '''Compute the C{Frechet} distance based on the I{angular} distance

431 in C{radians} from function L{pygeodesy.cosineLaw_}.

432

433 @see: L{FrechetEquirectangular}, L{FrechetEuclidean},

434 L{FrechetExact}, L{FrechetFlatPolar}, L{FrechetHaversine}

435 and L{FrechetVincentys}.

436

437 @note: See note at function L{pygeodesy.vincentys_}.

438 '''

439 _wrap = False

440

441 def __init__(self, points, wrap=False, fraction=None, name=NN):

442 '''New L{FrechetCosineLaw} calculator/interpolator.

443

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}).

452

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)

458

459 if _FOR_DOCS:

460 discrete = Frechet.discrete

461

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)

467

468

469 class FrechetDistanceTo(Frechet):

470 '''Compute the C{Frechet} distance based on the distance from the

471 points' C{LatLon.distanceTo} method, conventionally in C{meter}.

472

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

479

480 def __init__(self, points, fraction=None, name=NN, **distanceTo_kwds):

481 '''New L{FrechetDistanceTo} calculator/interpolator.

482

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.

492

493 @raise FrechetError: Insufficient number of B{C{points}} or an

494 invalid B{C{point}} or B{C{fraction}}.

495

496 @raise ImportError: Package U{geographiclib

497 <https://PyPI.org/project/geographiclib>} missing

498 iff B{C{points}} are L{ellipsoidalKarney.LatLon}s.

499

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

506

507 if _FOR_DOCS:

508 discrete = Frechet.discrete

509

510 def distance(self, p1, p2):

511 '''Return the distance in C{meter}.

512 '''

513 return p1.distanceTo(p2, **self._distanceTo_kwds)

514

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

524

525

526 class FrechetEquirectangular(FrechetRadians):

527 '''Compute the C{Frechet} distance based on the I{equirectangular}

528 distance in C{radians squared} like function L{pygeodesy.equirectangular_}.

529

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

536

537 def __init__(self, points, adjust=True, wrap=False, fraction=None, name=NN):

538 '''New L{FrechetEquirectangular} calculator/interpolator.

539

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}).

550

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)

556

557 if _FOR_DOCS:

558 discrete = Frechet.discrete

559

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

567

568

569 class FrechetEuclidean(FrechetRadians):

570 '''Compute the C{Frechet} distance based on the I{Euclidean}

571 distance in C{radians} from function L{pygeodesy.euclidean_}.

572

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

579

580 def __init__(self, points, adjust=True, fraction=None, name=NN):

581 '''New L{FrechetEuclidean} calculator/interpolator.

582

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}).

592

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)

598

599 if _FOR_DOCS:

600 discrete = Frechet.discrete

601

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)

606

607

608 class FrechetExact(FrechetDegrees):

609 '''Compute the C{Frechet} distance based on the I{angular}

610 distance in C{degrees} from method L{GeodesicExact}C{.Inverse}.

611

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

619

620 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):

621 '''New L{FrechetExact} calculator/interpolator.

622

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}).

633

634 @raise FrechetError: Insufficient number of B{C{points}} or

635 invalid B{C{fraction}}.

636

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

643

644 if _FOR_DOCS:

645 discrete = Frechet.discrete

646

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)

651

652

653 class 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_}.

656

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

664

665 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):

666 '''New L{FrechetFlatLocal}/L{FrechetHubeny} calculator/interpolator.

667

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}).

678

679 @raise FrechetError: Insufficient number of B{C{points}} or

680 invalid B{C{fraction}}.

681

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)

687

688 if _FOR_DOCS:

689 discrete = Frechet.discrete

690

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)

697

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

703

704

705 class FrechetFlatPolar(FrechetRadians):

706 '''Compute the C{Frechet} distance based on the I{angular} distance

707 in C{radians} from function L{flatPolar_}.

708

709 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular},

710 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetHaversine}

711 and L{FrechetVincentys}.

712 '''

713 _wrap = False

714

715 def __init__(self, points, wrap=False, fraction=None, name=NN):

716 '''New L{FrechetFlatPolar} calculator/interpolator.

717

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}).

726

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)

732

733 if _FOR_DOCS:

734 discrete = Frechet.discrete

735

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)

741

742

743 class FrechetHaversine(FrechetRadians):

744 '''Compute the C{Frechet} distance based on the I{angular}

745 distance in C{radians} from function L{pygeodesy.haversine_}.

746

747 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular},

748 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetFlatPolar}

749 and L{FrechetVincentys}.

750

751 @note: See note at function L{pygeodesy.vincentys_}.

752 '''

753 _wrap = False

754

755 def __init__(self, points, wrap=False, fraction=None, name=NN):

756 '''New L{FrechetHaversine} calculator/interpolator.

757

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}).

766

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)

772

773 if _FOR_DOCS:

774 discrete = Frechet.discrete

775

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)

781

782

783 class 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

789

790

791 class 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.

797

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.

804

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}).

815

816 @raise FrechetError: Insufficient number of B{C{points}} or

817 invalid B{C{fraction}}.

818

819 @raise ImportError: Package U{geographiclib

820 <https://PyPI.org/project/geographiclib>} missing.

821

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

828

829

830 class FrechetThomas(FrechetRadians):

831 '''Compute the C{Frechet} distance based on the I{angular} distance

832 in C{radians} from function L{pygeodesy.thomas_}.

833

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

840

841 def __init__(self, points, datum=None, wrap=False, fraction=None, name=NN):

842 '''New L{FrechetThomas} calculator/interpolator.

843

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}).

854

855 @raise FrechetError: Insufficient number of B{C{points}} or

856 invalid B{C{fraction}}.

857

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)

863

864 if _FOR_DOCS:

865 discrete = Frechet.discrete

866

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)

872

873

874 class FrechetVincentys(FrechetRadians):

875 '''Compute the C{Frechet} distance based on the I{angular}

876 distance in C{radians} from function L{pygeodesy.vincentys_}.

877

878 @see: L{FrechetCosineLaw}, L{FrechetEquirectangular},

879 L{FrechetEuclidean}, L{FrechetExact}, L{FrechetFlatPolar}

880 and L{FrechetHaversine}.

881

882 @note: See note at function L{pygeodesy.vincentys_}.

883 '''

884 _wrap = False

885

886 def __init__(self, points, wrap=False, fraction=None, name=NN):

887 '''New L{FrechetVincentys} calculator/interpolator.

888

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}).

897

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)

903

904 if _FOR_DOCS:

905 discrete = Frechet.discrete

906

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)

912

913

914 def _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 = {}

919

920 def iF(i): # cache index, depth ints and floats

921 return iFs.setdefault(i, i)

922

923 cF = _defaultdict(dict)

924

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

944

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)

951

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

959

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

964

965

966 def 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.

969

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}).

977

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}}.

981

982 @raise FrechetError: Insufficient number of B{C{points1}} or B{C{points2}}.

983

984 @raise RecursionError: Recursion depth exceeded, see U{sys.getrecursionlimit()

985 <https://docs.Python.org/3/library/sys.html#sys.getrecursionlimit>}.

986

987 @raise TypeError: If B{C{distance}} is not a callable.

988

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)

994

995 n1, ps1 = _points2(points1, closed=False, Error=FrechetError)

996 n2, ps2 = _points2(points2, closed=False, Error=FrechetError)

997

998 def dF(i1, i2):

999 return distance(ps1[i1], ps2[i2])

1000

1001 return _frechet_(n1, 1, n2, 1, dF, units)

1002

1003

1004 class 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}.

1010

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]}.

1017

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)

1023

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)

1030

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]

1038

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.