Coverage for pygeodesy/geodesici.py: 84%
447 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-27 20:21 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-27 20:21 -0400
2# -*- coding: utf-8 -*-
4u'''Class L{Intersector}, a pure Python version of parts of I{Karney}'s C++ class U{Intersect
5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Intersect.html>} to intersect
6geodesic lines.
8Only C++ member functions C{All}, C{Closest} and C{All} have been transcoded into Python as methods
9L{Intersector.All}, L{Intersector.Closest} and L{Intersector.Next} producing 4-item L{XDist}s.
11Adjacent methods L{Intersector.All5}, L{Intersector.Closest5}, L{Intersector.Next5} and
12L{Intersector.Next5s} return or yield L{Intersector5Tuple}s with the lat-, longitude, azimuth of
13each intersection as a C{Position} L{GDict} on each geodesic line.
15For more details, see the C++ U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/index.html>}
16documentation, I{Charles F.F. Karney}'s paper U{Geodesics intersections<https://arxiv.org/abs/2308.00495>}
17and I{S. Baselga Moreno & J.C. Martinez-Llario}'s U{Intersection and point-to-line solutions for geodesics
18on the ellipsoid<https://riunet.UPV.ES/bitstream/handle/10251/122902/Revised_Manuscript.pdf>}.
19'''
20# make sure int/int division yields float quotient
21from __future__ import division as _; del _ # PYCHOK semicolon
23from pygeodesy.basics import _copy, _enumereverse, map1, \
24 _xinstanceof, _xor
25from pygeodesy.constants import EPS, INF, INT0, PI, PI2, PI_4, _0_0, \
26 _0_5, _1_0, _1_5, _2_0, _3_0, _90_0
27from pygeodesy.ellipsoids import _EWGS84, Fmt
28from pygeodesy.errors import GeodesicError, IntersectionError, \
29 _xgeodesics, _xkwds_get
30from pygeodesy.fmath import euclid, favg, fdot
31from pygeodesy.fsums import Fsum, fsum1_, _ceil
32from pygeodesy.interns import _A_, _B_, _c_, _SPACE_, _too_
33from pygeodesy.karney import Caps, _diff182, _sincos2de
34from pygeodesy.lazily import _ALL_LAZY # _ALL_MODS as _MODS
35from pygeodesy.named import ADict, _NamedBase, _NamedTuple
36from pygeodesy.namedTuples import Degrees, Int, Meter, _Pass
37from pygeodesy.props import Property, Property_RO, property_RO
38# from pygeodesy.streprs import Fmt # from .ellipsoids
39# from pygeodesy.units import Degrees, Int, Meter # from .namedTuples
40from pygeodesy.utily import sincos2, atan2, fabs
42# from math import atan2, ceil as _ceil, fabs # .fsums, .utily
44__all__ = _ALL_LAZY.geodesici
45__version__ = '24.06.27'
47_0t = 0, # int
48_1_1t = -1, +1
49_1_0_1t = -1, 0, +1
50_EPS3 = EPS * _3_0
51_EPSr5 = pow(EPS, 0.2) # PYCHOK used!
52_TRIPS = 128
55def _L1(a, b):
56 return fabs(a) + fabs(b)
59class XDist(ADict):
60 '''4-Item result from L{Intersector.All}, L{Intersector.Closest} and
61 L{Intersector.Next} with the intersection offsets C{sA}, C{sB} and
62 C{sX0} in C{meter} and the coincidence indicator C{c}, an C{int},
63 +1 for parallel, -1 for anti-parallel, 0 otherwise.
64 '''
65 _Delta = EPS # default margin, see C{Intersector._Delto}
67 def __init__(self, sA=0, sB=0, c=0, sX0=INT0):
68 '''New L{XDist}.
70 @kwarg sA: Offset on geodesic line A (C{meter}).
71 @kwarg sB: Offset on geodesic line B (C{meter}).
72 @kwarg c: Coincidence indicator (C{int}, +1 for parallel
73 -1 for anti-parallel, 0 otherwise.
74 @kwarg sX0: Offset to C{X0} ({Cmeter}) or L{INT0}.
75 '''
76 self.set_(sA=sA, sB=sB, c=c, sX0=sX0)
78 def __add__(self, other):
79 X = _copy(self)
80 X += other
81 return X
83 def __eq__(self, other):
84 return not self.__ne__(other)
86 def __iadd__(self, other):
87 if isinstance(other, tuple): # and len(other) == 2:
88 a, b = other
89 else:
90 # _xinstanceof(XDist, other=other)
91 a = other.sA
92 b = other.sB
93 if other.c:
94 self.c = other.c
95 self.sA += a # PYCHOK sA
96 self.sB += b # PYCHOK sB
97 return self
99 def __le__(self, other):
100 # _xinstanceof(XDist, other=other)
101 return self == other or self < other
103 def __lt__(self, other):
104 # _xinstanceof(XDist, other=other)
105 return (self.sA < other.sA or (self.sA == other.sA and # PYCHOK sA
106 self.sB < other.sB) and self != other) # PYCHOK sB
108 def __ne__(self, other):
109 # _xinstanceof(XDist, other=other)
110 return self is not other and self.L1(other) > self._Delta
112 def _fixCoincident(self, X, *c0):
113 # return the mid-point if C{X} is anti-/parallel
114 c = c0[0] if c0 else X.c
115 if c:
116 s = (self.sA - X.sA + # PYCHOK sA
117 (self.sB - X.sB) * c) * _0_5 # PYCHOK sB
118 X = X + (s, s * c) # NOT +=
119 return X
121 def L1(self, other=None):
122 '''Return the C{L1} distance.
123 '''
124 a, b = self.sA, self.sB # PYCHOK sA, sB
125 if other is not None:
126 # _xinstanceof(XDist, other=other)
127 a -= other.sA
128 b -= other.sB
129 return _L1(a, b)
131 def _nD1(self, D1):
132 # yield the C{Closest} starts
133 D_ = 0, D1, -D1
134 for a, b in zip((0, 1, -1, 0, 0),
135 (0, 0, 0, 1, -1)):
136 yield self + (D_[a], D_[b])
138 def _nD2(self, D2):
139 # yield the C{Next} starts
140 D22 = D2 * _2_0
141 D_ = 0, D2, D22, -D22, -D2
142 for a, b in zip((-1, -1, 1, 1, -2, 0, 2, 0),
143 (-1, 1, -1, 1, 0, 2, 0, -2)):
144 yield self + (D_[a], D_[b])
146 def _nmD3(self, n, m, D3):
147 # yield the C{All} starts
148 for i in range(n, m, 2):
149 for j in range(n, m, 2):
150 if i or j:
151 yield self + ((i + j) * D3,
152 (i - j) * D3)
154 def _skip(self, S_, T1_Delta):
155 # remove starts from C{S_} near this C{XDist}
156 for j, S in _enumereverse(S_):
157 if S.L1(self) < T1_Delta:
158 S_.pop(j)
160_X000 = XDist() # PYCHOK origin
161_XINF = XDist(INF)
164class Intersector(_NamedBase):
165 '''Finder of intersections between two goedesic lines, each an instance
166 of L{GeodesicLineExact<pygeodesy.geodesicx.GeodesicLineExact>},
167 wrapped L{GeodesicLine<pygeodesy.geodesicw.GeodesicLine>} or
168 L{GeodesicLineSolve<pygeodesy.geodsolve.GeodesicLineSolve>}.
170 @see: I{Karney}'s C++ class U{Intersect<https://GeographicLib.sourceforge.io/
171 C++/doc/classGeographicLib_1_1Intersect.html#details>} for more details.
172 '''
173 # _D1 = 0
174 # _D2 = 0
175 # _g = None
176 # _T1 = 0
177 # _T5 = 0
179 def __init__(self, geodesic, **name):
180 '''New L{Intersector}.
182 @arg geodesic: The geodesic (L{GeodesicExact<pygeodesy.geodesicx.GeodesicExact>},
183 wrapped L{Geodesic<pygeodesy.geodesicw.Geodesic>} or
184 L{GeodesicSolve<pygeodesy.geodsolve.GeodesicSolve>}).
185 @kwarg name: Optional C{B{name}=NN} (C{str}).
187 @raise GeodesicError: The eccentricity of the B{C{geodesic}}'s ellipsoid is too
188 large or no initial convergence.
190 @see: The B{Note} at I{Karney}'s C++ U{Intersect<https://GeographicLib.sourceforge.io/
191 C++/doc/classGeographicLib_1_1Intersect.html#ae41f54c9a44836f6c8f140f6994930cf>}.
192 '''
193 _xinstanceof(*_EWGS84._Geodesics, geodesic=geodesic)
194 self._g = geodesic
195 if name:
196 self.name = name
197 E = self.ellipsoid
199 t1 = E.b * PI # min distance between intersects
200 t2 = self._polarDist2(_90_0)[0] * _2_0 # furthest closest intersect
201 t5 = self._Inversa12(_90_0)[0] * _2_0 # longest shortest geodesic
202 if self.f > 0:
203 t3 = self._obliqDist4()[0]
204 t4 = t1
205 else: # PYCHOK no cover
206 t1, t2, t3 = t2, t1, t5
207 t4 = self._polarB3()[0]
208 d1 = t2 * _0_5
209 d2 = t3 / _1_5
210 d3 = t4 - self.Delta
211 t2 = t1 * _2_0
212 if not (d1 < d3 and d2 < d3 and d2 < t2):
213 t = Fmt.PARENSPACED(_too_('eccentric'), E.e)
214 raise GeodesicError(ellipsoid=E.toStr(terse=2), txt=t)
215 self._D1 = d1 # tile spacing for Closest
216 self._D2 = d2 # tile spacing for Next
217 self._D3 = d3 # tile spacing for All
218 self._T1 = t1 # min distance between intersects
219 self._T2 = t2
220# self._T5 = t5
222 @Property_RO
223 def a(self):
224 '''Get the I{equatorial} radius, semi-axis (C{meter}).
225 '''
226 return self.ellipsoid.a
228 equatoradius = a # = Requatorial
230 def All(self, glA, glB, X0=_X000, **sMax):
231 '''Find all intersection of two geodesic lines up to a limit.
233 @arg glA: A geodesic line (L{Line<Intersector.Line>}).
234 @arg glB: An other geodesic line (L{Line<Intersector.Line>}).
235 @kwarg X0: Optional offsets along the geodesic lines (L{XDist}).
236 @kwarg sMax: Optional, upper limit C{B{sMax}=2*PI*R} for the
237 distance (C{meter}).
239 @return: Yield an L{XDist} for each intersection found.
241 @raise GeodesicError: Geodesic line B{C{glA}} or B{C{glB}}
242 invalid, incompatible or ill-configured.
244 @raise IntersectionError: No convergence.
245 '''
246 self._xLines(glA, glB)
247 sMax = _xkwds_get(sMax, sMax=self.R * PI2)
248 if sMax < _EPS3:
249 sMax = _EPS3 # raise GeodesicError(sMax=sMax)
251 D, _D = self.Delta, self._C_2
252 xMax = sMax + D
253 m = int(_ceil(xMax / self._D3)) # m x m tiles
254 d3 = xMax / m
255 T2d3D = self._T2d3Delta(d3)
256 _X0fx = X0._fixCoincident
258 c0 = 0
259 C_ = _List(D) # closest coincident
260 X_ = _List(D) # intersections found
261 S_ = list(X0._nmD3(1 - m, m, d3 * _0_5))
262 # assert len(s_) + 1 == m * m + (m - 1) % 2
263 while S_:
264 Q, i = self._Basic2(glA, glB, S_.pop(0))
265 if Q in X_:
266 continue
267 # assert Q.c == c0 or not c0
268 a, c0 = len(X_), Q.c
269 if c0: # coincident intersection
270 Q = _X0fx(Q)
271 if Q in C_:
272 continue
273 C_.addend(Q)
274 # elimate all existing intersections
275 # on this line (which didn't set c0)
276 for j, X in _enumereverse(X_):
277 if _X0fx(X, c0).L1(Q) <= D: # X' == Q
278 X_.pop(j)
280 a, s0 = len(X_), Q.sA
281 args = self._m12_M12_M21(glA, s0)
282 _cjD = self._conjDist
283 for s in (-_D, _D):
284 s += s0
285 sa = 0
286 while True:
287 sa = _cjD(glA, s + sa, *args) - s0
288 X = Q + (sa, sa * c0)
289 i += 1
290 if X_.addend(X, X0.L1(X), i) > xMax:
291 break
293 X_.addend(Q, X0.L1(Q), i + 1)
294 for X in X_[a:]: # addended Xs
295 X._skip(S_, T2d3D)
297 return X_.sortrim(X0, sMax) # generator!
299 def All5(self, glA, glB, X0=_X000, aMax=0, **sMax):
300 '''Find all intersection of two geodesic lines up to a limit.
302 @kwarg aMax: Upper limit for the angular distance (C{degrees})
303 or C{None} or C{0} for unlimited.
305 @return: Yield an L{Intersector5Tuple}C{(A, B, sAB, aAB, c)}
306 for each intersection found.
308 @see: Methods L{All} for further details.
309 '''
310 aA = aB = _0_0
311 for X in self.All(glA, glB, X0=X0, **sMax):
312 r = self._In5T(glA, glB, X, X)
313 yield r
314 if aMax:
315 aA += r.A.a12
316 aB += r.B.a12
317 if fabs(aA) > aMax or fabs(aB) > aMax:
318 break
320 def _Basic2(self, glA, glB, S, i=0):
321 '''(INTERNAL) Get a basic solution.
322 '''
323 X = _copy(S)
324 for _ in range(_TRIPS):
325 S = self._Spherical(glA, glB, X)
326 X += S
327 i += 1
328 if X.c or S.L1() <= self._Tol: # or isnan
329 return self._Delto(X), i
331 raise IntersectionError(Fmt.no_convergence(S.L1(), self._Tol))
333 @Property_RO
334 def _C_2(self): # normalizer, semi-circumference, C++ _d
335 return self.R * PI # ~20K Km WGS84
337 def Closest(self, glA, glB, X0=_X000):
338 '''Find the closest intersection of two geodesic lines.
340 @arg glA: A geodesic line (L{Line<Intersector.Line>}).
341 @arg glB: An other geodesic line (L{Line<Intersector.Line>}).
342 @kwarg X0: Optional offsets along the geodesic lines (L{XDist}).
344 @return: The intersection (L{XDist}) or C{None} if none found.
346 @raise GeodesicError: Geodesic line B{C{glA}} or B{C{glB}}
347 invalid, incompatible or ill-configured.
349 @raise IntersectionError: No convergence.
350 '''
351 self._xLines(glA, glB)
352 Q, d, S_, i = X0, INF, list(X0._nD1(self._D1)), 0
353 while S_:
354 X, i = self._Basic2(glA, glB, S_.pop(0), i)
355 X = X0._fixCoincident(X)
356 if X.L1(Q) > self.Delta: # X != Q
357 d0 = X.L1(X0)
358 if d0 < self._T1:
359 Q, d, q = X, d0, i
360 break
361 if d0 < d or Q is X0:
362 Q, d, q = X, d0, i
363 X._skip(S_, self._T2D1Delta)
365 return None if Q is X0 else Q.set_(sX0=d, iteration=q)
367 def Closest5(self, glA, glB, X0=_X000):
368 '''Find the closest intersection of two geodesic lines.
370 @return: An L{Intersector5Tuple}C{(A, B, sAB, aAB, c)}
371 or C{None} if none found.
373 @see: Method L{Closest} for further details.
374 '''
375 X = self.Closest(glA, glB, X0=X0)
376 return X if X is None else self._In5T(glA, glB, X, X)
378 def _conjDist(self, gl, s, m12=0, M12=1, M21=1, semi=False):
379 # Find semi-/conjugate point relative to s0 which is close to s1.
380 # if semi:
381 # solve for M23 = 0 using dM23 / ds3 = - (1 - M23 * M32) / m23
382 # else:
383 # solve for m23 = 0 using dm23 / ds3 = M32
384 _S2, _abs, _1 = Fsum(s).fsum2_, fabs, _1_0
385 for _ in range(_TRIPS):
386 m13, M13, M31 = self._m12_M12_M21(gl, s)
387 # see "Algorithms for geodesics", eqs. 31, 32, 33.
388 m23 = m13 * M12
389 M32 = M31 * M12
390 if m12:
391 m23 -= m12 * M13
392 if m13:
393 M32 += (_1 - M13 * M31) * m12 / m13
394 if semi:
395 M23 = M13 * M21
396 # when m12 -> eps, (1 - M12 * M21) -> eps^2, I suppose.
397 if m12 and m13:
398 M23 += (_1 - M12 * M21) * m13 / m12
399 d = m23 * M23 / (_1 - M23 * M32)
400 else:
401 d = -m23 / M32
402 s, d = _S2(d)
403 if _abs(d) <= self._Tol:
404 break
405 return s
407 _gl3 = None
409 @Property
410 def _conjDist3s(self):
411 gl, self._gl3, _D = self._gl3, None, self._C_2
412 return tuple(self._conjDist(gl, s) for s in (-_D, 0, _D))
414 @_conjDist3s.setter # PYCHOK setter!
415 def _conjDist3(self, gl):
416 # _XLines(gl, gl)
417 self._gl3 = gl
419 def _conjDist3Tt_(self, c, X0=_X000):
420 for s in self._conjDist3s:
421 T = XDist(s, s * c, c)
422 yield self._Delto(T), T.L1(X0)
424 def _conjDist5(self, azi):
425 gl = self._Line(azi1=azi)
426 s = self._conjDist(gl, self._C_2)
427 X, _ = self._Basic2(gl, gl, XDist(s * _0_5, -s * _1_5))
428 return s, (X.L1() - s * _2_0), azi, X.sA, X.sB
430 @Property_RO
431 def Delta(self):
432 '''Get the equality and tiling margin (C{meter}).
433 '''
434 return self._C_2 * _EPSr5 # ~15 Km WGS84
436 def _Delto(self, X):
437 # copy Delta into X, overriding X's default
438 X._Delta = self.Delta # NOT X.set_(self.Delta)
439 return X
441 @Property_RO
442 def ellipsoid(self):
443 '''Get the C{geodesic}'s ellipsoid (C{Ellipsoid}).
444 '''
445 return self.geodesic.datum.ellipsoid
447 @Property_RO
448 def _EPS3R(self):
449 return _EPS3 * self.R
451 @Property_RO
452 def f(self):
453 '''Get the I{flattening} (C{scalar}), C{0} for spherical, negative for prolate.
454 '''
455 return self.ellipsoid.f
457 flattening = f
459 @Property_RO
460 def _faPI_4(self):
461 return (self.f + _2_0) * self.a * PI_4
463 @property_RO
464 def geodesic(self):
465 '''Get the C{geodesic} (C{Geodesic...}).
466 '''
467 return self._g
469 @Property_RO
470 def _GeodesicLines(self):
471 '''(INTERNAL) Get the C{Geodesic...Line} class(es).
472 '''
473 return type(self._Line()),
475 def _In5T(self, glA, glB, S, X):
476 # Return an intersection as C{Intersector5Tuple}.
477 A = self._Position(glA, S.sA, S.sX0)
478 B = self._Position(glB, S.sB, S.sX0)
479 s, a = self._Inversa12(A, B)
480 r = Intersector5Tuple(A, B, s, a, X.c, iteration=X.iteration)
481 return r
483 def _Inversa12(self, A, B=None):
484 lls = (0, 0, A, 0) if B is None else (A.lat2, A.lon2,
485 B.lat2, B.lon2)
486 r = self._g.Inverse(*lls, outmask=Caps.DISTANCE)
487 return r.s12, r.a12 # .a12 always in r
489 def _Inverse(self, A, B): # caps=Caps.STANDARD
490 return self._g.Inverse(A.lat2, A.lon2, B.lat2, B.lon2)
492 def Line(self, lat1, lon1, azi_lat2, *lon2, **name):
493 '''Return a geodesic line from this C{Intersector}'s geodesic, specified by
494 two (goedetic) points or a (goedetic) point and an (initial) azimuth.
496 @arg lat1: Latitude of the first point (C{degrees}).
497 @arg lon1: Longitude of the first point (C{degrees}).
498 @arg azi_lat2: Azimuth at the first point (compass C{degrees}) if no
499 B{C{lon2}} argument is given, otherwise the latitude of
500 the second point (C{degrees}).
501 @arg lon2: If given, the longitude of the second point (C{degrees}).
502 @kwarg name: Optional C{B{name}=NN} (C{str}).
504 @return: A line (from L{geodesic<Intersector.geodesic>}C{.Line} or
505 C{-.InverseLine}), properly configured for L{Intersector}.
506 '''
507 args = (lat1, lon1, azi_lat2) + lon2
508 gl = self._g.InverseLine(*args, caps=Caps.LINE_CAPS) if lon2 else \
509 self._g.Line( *args, caps=Caps.LINE_CAPS)
510 if name:
511 gl.name= name
512 return gl
514 def _Line(self, lat1=0, lon1=0, azi1=0):
515 return self._g.Line(lat1, lon1, azi1, caps=Caps.LINE_CAPS)
517 def _m12_M12_M21(self, gl, s):
518 P = gl.Position(s, outmask=Caps._REDUCEDLENGTH_GEODESICSCALE)
519 return P.m12, P.M12, P.M21
521 def Next(self, glA, glB, **eps):
522 '''Yield the next intersection of two I{intersecting} geodesic lines.
524 @arg glA: A geodesic line (L{Line<Intersector.Line>}).
525 @arg glB: An other geodesic line (L{Line<Intersector.Line>}).
526 @kwarg eps: Optional equality margin C{B{eps}=Delta} (C{degrees}).
528 @return: The intersection (L{XDist}) or C{None} if none found.
530 @raise GeodesicError: Geodesic line B{C{glA}} or B{C{glB}}
531 invalid, incompatible, ill-configured or
532 C{(lat1, lon1)} not B{C{eps}}-equal.
534 @raise IntersectionError: No convergence.
536 @note: Offset C{X0} is implicit, zeros.
537 '''
538 self._xLines(glA, glB)
539 e = _xkwds_get(eps, eps=self.Delta)
540 a = glA.lat1 - glB.lat1
541 b = glA.lon1 - glB.lon1
542 if fabs(a) > e or fabs(b) > e:
543 raise GeodesicError(lat1=a, lon1=b, eps=e)
544 return self._Next(glA, glB)
546 def Next5(self, glA, glB, eps=_EPS3):
547 '''Yield the next intersection of two I{intersecting} geodesic lines.
549 @return: An L{Intersector5Tuple}C{(A, B, sAB, aAB, c)} or C{None}
550 if none found.
552 @see: Method L{Next} for further details.
553 '''
554 X = self.Next(glA, glB, eps=eps)
555 return X if X is None else self._In5T(glA, glB, X, X)
557 def Next5s(self, glA, glB, X0=_X000, aMax=1801, sMax=0, avg=False, **Delta):
558 '''Yield C{Next} intersections up to a maximal (angular) distance.
560 @kwarg aMax: Upper limit for the angular distance (C{degrees}) or
561 C{None} or C{0} for unlimited.
562 @kwarg sMax: Upper limit for the distance (C{meter}) or C{None} or
563 C{0} for unlimited.
564 @kwarg avg: If C{True}, set the next intersection lat- and longitude
565 to the mid-point of the previous ones (C{bool}).
566 @kwarg Delta: Optional, margin overrding this margin (C{meter}), see
567 prpoerty L{Delta<Intersector.Delta>}.
569 @return: Yield an L{Intersector5Tuple}C{(A, B, sAB, aAB, c)} for
570 every intersection found.
572 @see: Methods L{Next5} for further details.
573 '''
574 X = self.Closest(glA, glB, X0=X0)
575 if X is not None:
576 D = _xkwds_get(Delta, Delta=self.Delta)
577 S, _L, _abs = X, self._Line, fabs
578 while True:
579 A, B, _, _, _ = r = self._In5T(glA, glB, S, X)
580 yield r
581 if (aMax and (_abs(A.a12) > aMax or _abs(B.a12) > aMax)) or \
582 (sMax and (_abs(A.s12) > sMax or _abs(B.s12) > sMax)):
583 break
584 latA, lonA = A.lat2, A.lon2
585 latB, lonB = B.lat2, B.lon2
586 if avg:
587 latA = latB = favg(latA, latB)
588 lonA = lonB = favg(lonA, lonB)
589 X = self._Next(_L(latA, lonA, A.azi2),
590 _L(latB, lonB, B.azi2))
591 if X is None or X.L1() < D:
592 break
593 S += X.sA, X.sB
594 S.set_(sX0=X.sX0 + S.sX0)
596 def _Next(self, glA, glB):
597 '''(INTERNAL) Find the next intersection.
598 '''
599 X0, self._conjDist3s = _X000, glA
600 Q, d, S_, i = _XINF, INF, list(X0._nD2(self._D2)), 0
601 while S_:
602 X, i = self._Basic2(glA, glB, S_.pop(0), i)
603 X = X0._fixCoincident(X)
604 t = X.L1(X0) # == X.L1()
605 c, z = X.c, (t <= self.Delta) # X == X0
606 if z:
607 if not c:
608 continue
609 Tt_ = self._conjDist3Tt_(c, X0)
610 else:
611 Tt_ = (X, t),
613 for T, t in Tt_:
614 if t < d or Q is _XINF:
615 Q, d, q = T, t, i
616 i += 1
618 for s in ((_1_1t if z else _1_0_1t)
619 if c else _0t):
620 T = X
621 if s and c:
622 s *= self._D2
623 T = X + (s, s * c) # NOT +=
624 T._skip(S_, self._T2D2Delta)
626 return None if Q is _XINF else Q.set_(sX0=d, iteration=q)
628 def _obliqDist4(self):
629 zx = 45.0
630 if self.f:
631 _abs, _cjD5 = fabs, self._conjDist5
633 _, ds0, z0, _, _ = _cjD5(zx + _1_0)
634 s1, ds1, z1, sAx, sBx = _cjD5(zx - _1_0)
635 sx, dsx, zx = s1, _abs(ds1), z1
636 # find ds(azi) = 0 by secant method
637 for _ in range(16):
638 if ds1 == ds0:
639 break
640 z = (z0 * ds1 - z1 * ds0) / (ds1 - ds0)
641 _, ds0, z0 = s1, ds1, z1
642 s1, ds1, z1, a, b = _cjD5(z)
643 if _abs(ds1) < dsx:
644 sx, dsx, zx, sAx, sBx = s1, _abs(ds1), z, a, b
645 if not dsx:
646 break
647 else:
648 sx, sAx, sBx = self._C_2, _0_5, -_1_5
649 return sx, zx, sAx, sBx
651 def _polarB3(self, lats=False): # PYCHOK no cover
652 latx = 64.0
653 lat = _90_0 - latx
654 if self.f:
655 _d, _pD2 = fdot, self._polarDist2
657 s0, lat0 = _pD2(latx - _1_0)
658 s1, lat1 = _pD2(latx + _1_0)
659 s2, lat2 = \
660 sx, latx = _pD2(latx)
661 prolate = self.f < 0
662 # solve for ds(lat) / dlat = 0 with a quadratic fit
663 for _ in range(_TRIPS):
664 t = (lat1 - lat0), (lat0 - lat2), (lat2 - lat1)
665 d = _d(t, s2, s1, s0) * _2_0
666 if not d: # or isnan(d)
667 break
668 lat = _d(t, (lat1 + lat0) * s2,
669 (lat0 + lat2) * s1,
670 (lat2 + lat1) * s0) / d
671 s0, lat0 = s1, lat1
672 s1, lat1 = s2, lat2
673 s2, lat2 = _pD2(lat)
674 if (s2 < sx) if prolate else (s2 > sx):
675 sx, latx = s2, lat2
676 if lats:
677 _, lat = _pD2(latx, lat2=True)
678 sx += sx
679 else:
680 sx = self._C_2
681 return sx, latx, lat
683 def _polarDist2(self, lat1, lat2=False):
684 gl = self._Line(lat1=lat1)
685 s = self._conjDist(gl, self._faPI_4, semi=True)
686 if lat2:
687 lat1 = gl.Position(s, outmask=Caps.LATITUDE).lat2
688 return s, lat1
690 def _Position(self, gl, s, *sX0):
691 P = gl.Position(s, outmask=Caps._STD_LINE)
692 if sX0:
693 X = gl.Position(*sX0, outmask=Caps._STD_LINE)
694 P.set_(lat0=X.lat2, lon0=X.lon2,
695 azi0=X.azi2, s10=X.s12, a10=X.a12)
696 return P
698 @Property_RO
699 def R(self):
700 '''Get the I{authalic} earth radius (C{meter}).
701 '''
702 return self.ellipsoid.R2
704 def _Spherical(self, glA, glB, S):
705 '''(INTERNAL) Get solution based from a spherical triangle.
706 '''
707 # threshold for coincident geodesics/intersections ~4.3 nm WGS84.
708 A = self._Position(glA, S.sA)
709 B = self._Position(glB, S.sB)
710 D = self._Inverse(A, B)
712 a, da = _diff182(A.azi2, D.azi1) # interior angle at A
713 b, db = _diff182(B.azi2, D.azi2) # exterior angle at B
714 c, dc = _diff182(a, b)
715 if fsum1_(dc, db, -da, c) < 0: # inverted triangle
716 a, da = -a, -da
717 b, db = -b, -db
718 sa, ca = _sincos2de(a, da)
719 sb, cb = _sincos2de(b, db)
721 e, z, _abs = _EPS3, D.s12, fabs
722 if _abs(z) <= self._EPS3R: # XXX z <= ...
723 sA = sB = 0 # at intersection
724 c = 1 if _abs(sa - sb) <= e and _abs(ca - cb) <= e else (
725 -1 if _abs(sa + sb) <= e and _abs(ca + cb) <= e else 0)
726 elif _abs(sa) <= e and _abs(sb) <= e: # coincident
727 sA = ca * z * _0_5 # choose mid-point
728 sB = -cb * z * _0_5
729 c = 1 if (ca * cb) > 0 else -1
730 # alt1: sA = ca * z; sB = 0
731 # alt2: sB = -cb * z; sA = 0
732 else: # general case
733 sz, cz = sincos2(z / self.R)
734 # [SKIP: Divide args by |sz| to avoid possible underflow
735 # in {sa, sb} * sz; this is probably not necessary].
736 # Definitely need to treat sz < 0 (z > PI*R) correctly in
737 # order to avoid some convergence failures in _Basic2.
738 sA = atan2(sb * sz, sb * ca * cz - cb * sa) * self.R
739 sB = atan2(sa * sz, -sa * cb * cz + ca * sb) * self.R
740 c = 0
741 return XDist(sA, sB, c) # no ._Delto
743 @Property_RO
744 def _T2D1Delta(self):
745 return self._T2d3Delta(self._D1)
747 @Property_RO
748 def _T2D2Delta(self):
749 return self._T2d3Delta(self._D2)
751 def _T2d3Delta(self, d3):
752 return self._T2 - d3 - self.Delta
754 @Property_RO
755 def _Tol(self): # convergence tolerance
756 return self._C_2 * pow(EPS, 0.75) # _0_75
758 def toStr(self, **prec_sep_name): # PYCHOK signature
759 '''Return this C{Intersector} as string.
761 @see: L{Ellipsoid.toStr<pygeodesy.ellipsoids.Ellipsoid.toStr>}
762 for further details.
764 @return: C{Intersector} (C{str}).
765 '''
766 return self._instr(props=(Intersector.geodesic,), **prec_sep_name)
768 def _xLines(self, glA, glB):
769 # check two geodesic lines vs this geodesic
770 C, gls = Caps.LINE_CAPS, dict(glA=glA, glB=glB)
771 _xinstanceof(*self._GeodesicLines, **gls)
772 for n, gl in gls.items():
773 try:
774 _xgeodesics(gl.geodesic, self.geodesic)
775 c = gl.caps & C
776 if c != C: # not gl.caps_(C)
777 c, C, x = map1(bin, c, C, _xor(c, C))
778 x = _SPACE_(_xor.__name__, repr(x))[1:]
779 raise GeodesicError(caps=c, LINE_CAPS=C, txt=x)
780 except Exception as x:
781 raise GeodesicError(n, gl, cause=x)
784class Intersector5Tuple(_NamedTuple):
785 '''5-Tuple C{(A, B, sAB, aAB, c)} with C{A} and C{B} the C{Position}
786 of the intersection on each geodesic line, the distance C{sAB}
787 between C{A} and C{B} in C{meter}, angular distance C{aAB} in
788 C{degrees} and coincidence indicator C{c} (C{int}), see L{XDist}.
790 @note: C{A} and C{B} are each a C{GeodesicLine...Position} for
791 C{outmask=Caps.STANDARD} with the intersection location in
792 C{lat2}, C{lon2}, azimuth in C{azi2}, the distance C{s12}
793 in C{meter} and angular distance C{a12} in C{degrees} and
794 extended with the C{X0} offset location in C{lat0}, C{lon0},
795 C{azi0}, C{s10} and C{a10}.
796 '''
797 _Names_ = (_A_, _B_, 'sAB', 'aAB', _c_)
798 _Units_ = (_Pass, _Pass, Meter, Degrees, Int)
801class _List(list):
803 _Delta = 0 # equality margin
805 def __init__(self, Delta):
806 self._Delta = Delta
807# list.__init__(self)
809 def __contains__(self, other):
810 # handle C{if X in this: ...}
811 a, b = other.sA, other.sB
812 D, _D1 = self._Delta, _L1
813 for X in self:
814 if _D1(X.sA - a, X.sB - b) <= D:
815 return True
816 return False
818 def addend(self, X, *d0_i):
819 # append an item, updated
820 if d0_i:
821 d0, i = d0_i
822 X.set_(sX0=d0, iteration=i)
823 self.append(X)
824 return X.sX0
826 def sortrim(self, X0, sMax):
827 # trim and sort the X items
829 def _key(Xk):
830 _, k = Xk
831 return k # rank of X
833 for X, _ in sorted(self.trim(X0, sMax), key=_key):
834 yield X # de-tuple (X, k)
836 def trim(self, X0, sMax):
837 # trim and yield 2-tuple (X, rank)
838 a, b, _eu = X0.sA, X0.sB, euclid
840 for X in self:
841 k = X.sX0
842 if k <= sMax:
843 k += _eu(X.sA - a, X.sB - b)
844 yield X, k # rank of X
847if __name__ == '__main__':
849 from pygeodesy import GeodesicExact, printf
851 I = Intersector(GeodesicExact(), name='Test') # PYCHOK I
853 # <https://GeographicLib.sourceforge.io/C++/doc/classGeographicLib_1_1Intersect.html>
854 a = I.Line( 0, 0, 45)
855 b = I.Line(45, 10, 135)
856 printf('Closest: %r', I.Closest(a, b))
857 printf('Closest5: %r', I.Closest5(a, b), nt=1)
859 for i, t in enumerate(I.Next5s(a, b)):
860 printf('Next5s %s: %r (%s)', i, t, t.iteration)
862 # <https://GeographicLib.sourceforge.io/C++/doc/IntersectTool.1.html>
863 a = I.Line(50, -4, -147.7)
864 b = I.Line( 0, 0, 90)
865 printf('Closest: %r', I.Closest(a, b), nl=1)
866 printf('Closest5: %r', I.Closest5(a, b), nt=1)
868 a = I.Line(50, -4, -147.7)
869 b = I.Line( 0, 180, 0)
870 for i, X in enumerate(I.All(a, b)):
871 printf('All %s: %r (%s)', i, X, X.iteration)
872 if i > 9:
873 break
874 printf('')
875 for i, t in enumerate(I.All5(a, b)):
876 printf('All5 %s: %r (%s)', i, t, t.iteration)
877 if i > 9:
878 break
880 # % echo 50N 4W 147.7W 0 0 90 | IntersectTool -e 6371e3 0 -c -p 0 -C
881 # 6077191 -3318019 0
882 # -0.00000 -29.83966 -0.00000 -29.83966 0
883 I = Intersector(GeodesicExact(6371e3, 0), name='Test') # PYCHOK I
884 a = I.Line(50, -4, -147.7)
885 b = I.Line( 0, 0, 90)
886 printf('Closest: %r', I.Closest(a, b), nl=1)
887 printf('Closest5: %r', I.Closest5(a, b))
889# **) MIT License
890#
891# Copyright (C) 2024-2024 -- mrJean1 at Gmail -- All Rights Reserved.
892#
893# Permission is hereby granted, free of charge, to any person obtaining a
894# copy of this software and associated documentation files (the "Software"),
895# to deal in the Software without restriction, including without limitation
896# the rights to use, copy, modify, merge, publish, distribute, sublicense,
897# and/or sell copies of the Software, and to permit persons to whom the
898# Software is furnished to do so, subject to the following conditions:
899#
900# The above copyright notice and this permission notice shall be included
901# in all copies or substantial portions of the Software.
902#
903# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
904# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
905# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
906# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
907# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
908# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
909# OTHER DEALINGS IN THE SOFTWARE.