Coverage for pygeodesy/geodesici.py: 84%

447 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-06-27 20:21 -0400

1 

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

3 

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. 

7 

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. 

10 

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. 

14 

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 

22 

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 

41 

42# from math import atan2, ceil as _ceil, fabs # .fsums, .utily 

43 

44__all__ = _ALL_LAZY.geodesici 

45__version__ = '24.06.27' 

46 

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 

53 

54 

55def _L1(a, b): 

56 return fabs(a) + fabs(b) 

57 

58 

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} 

66 

67 def __init__(self, sA=0, sB=0, c=0, sX0=INT0): 

68 '''New L{XDist}. 

69 

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) 

77 

78 def __add__(self, other): 

79 X = _copy(self) 

80 X += other 

81 return X 

82 

83 def __eq__(self, other): 

84 return not self.__ne__(other) 

85 

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 

98 

99 def __le__(self, other): 

100 # _xinstanceof(XDist, other=other) 

101 return self == other or self < other 

102 

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 

107 

108 def __ne__(self, other): 

109 # _xinstanceof(XDist, other=other) 

110 return self is not other and self.L1(other) > self._Delta 

111 

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 

120 

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) 

130 

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

137 

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

145 

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) 

153 

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) 

159 

160_X000 = XDist() # PYCHOK origin 

161_XINF = XDist(INF) 

162 

163 

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

169 

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 

178 

179 def __init__(self, geodesic, **name): 

180 '''New L{Intersector}. 

181 

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

186 

187 @raise GeodesicError: The eccentricity of the B{C{geodesic}}'s ellipsoid is too 

188 large or no initial convergence. 

189 

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 

198 

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 

221 

222 @Property_RO 

223 def a(self): 

224 '''Get the I{equatorial} radius, semi-axis (C{meter}). 

225 ''' 

226 return self.ellipsoid.a 

227 

228 equatoradius = a # = Requatorial 

229 

230 def All(self, glA, glB, X0=_X000, **sMax): 

231 '''Find all intersection of two geodesic lines up to a limit. 

232 

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

238 

239 @return: Yield an L{XDist} for each intersection found. 

240 

241 @raise GeodesicError: Geodesic line B{C{glA}} or B{C{glB}} 

242 invalid, incompatible or ill-configured. 

243 

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) 

250 

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 

257 

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) 

279 

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 

292 

293 X_.addend(Q, X0.L1(Q), i + 1) 

294 for X in X_[a:]: # addended Xs 

295 X._skip(S_, T2d3D) 

296 

297 return X_.sortrim(X0, sMax) # generator! 

298 

299 def All5(self, glA, glB, X0=_X000, aMax=0, **sMax): 

300 '''Find all intersection of two geodesic lines up to a limit. 

301 

302 @kwarg aMax: Upper limit for the angular distance (C{degrees}) 

303 or C{None} or C{0} for unlimited. 

304 

305 @return: Yield an L{Intersector5Tuple}C{(A, B, sAB, aAB, c)} 

306 for each intersection found. 

307 

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 

319 

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 

330 

331 raise IntersectionError(Fmt.no_convergence(S.L1(), self._Tol)) 

332 

333 @Property_RO 

334 def _C_2(self): # normalizer, semi-circumference, C++ _d 

335 return self.R * PI # ~20K Km WGS84 

336 

337 def Closest(self, glA, glB, X0=_X000): 

338 '''Find the closest intersection of two geodesic lines. 

339 

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

343 

344 @return: The intersection (L{XDist}) or C{None} if none found. 

345 

346 @raise GeodesicError: Geodesic line B{C{glA}} or B{C{glB}} 

347 invalid, incompatible or ill-configured. 

348 

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) 

364 

365 return None if Q is X0 else Q.set_(sX0=d, iteration=q) 

366 

367 def Closest5(self, glA, glB, X0=_X000): 

368 '''Find the closest intersection of two geodesic lines. 

369 

370 @return: An L{Intersector5Tuple}C{(A, B, sAB, aAB, c)} 

371 or C{None} if none found. 

372 

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) 

377 

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 

406 

407 _gl3 = None 

408 

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

413 

414 @_conjDist3s.setter # PYCHOK setter! 

415 def _conjDist3(self, gl): 

416 # _XLines(gl, gl) 

417 self._gl3 = gl 

418 

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) 

423 

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 

429 

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 

435 

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 

440 

441 @Property_RO 

442 def ellipsoid(self): 

443 '''Get the C{geodesic}'s ellipsoid (C{Ellipsoid}). 

444 ''' 

445 return self.geodesic.datum.ellipsoid 

446 

447 @Property_RO 

448 def _EPS3R(self): 

449 return _EPS3 * self.R 

450 

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 

456 

457 flattening = f 

458 

459 @Property_RO 

460 def _faPI_4(self): 

461 return (self.f + _2_0) * self.a * PI_4 

462 

463 @property_RO 

464 def geodesic(self): 

465 '''Get the C{geodesic} (C{Geodesic...}). 

466 ''' 

467 return self._g 

468 

469 @Property_RO 

470 def _GeodesicLines(self): 

471 '''(INTERNAL) Get the C{Geodesic...Line} class(es). 

472 ''' 

473 return type(self._Line()), 

474 

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 

482 

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 

488 

489 def _Inverse(self, A, B): # caps=Caps.STANDARD 

490 return self._g.Inverse(A.lat2, A.lon2, B.lat2, B.lon2) 

491 

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. 

495 

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

503 

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 

513 

514 def _Line(self, lat1=0, lon1=0, azi1=0): 

515 return self._g.Line(lat1, lon1, azi1, caps=Caps.LINE_CAPS) 

516 

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 

520 

521 def Next(self, glA, glB, **eps): 

522 '''Yield the next intersection of two I{intersecting} geodesic lines. 

523 

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

527 

528 @return: The intersection (L{XDist}) or C{None} if none found. 

529 

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. 

533 

534 @raise IntersectionError: No convergence. 

535 

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) 

545 

546 def Next5(self, glA, glB, eps=_EPS3): 

547 '''Yield the next intersection of two I{intersecting} geodesic lines. 

548 

549 @return: An L{Intersector5Tuple}C{(A, B, sAB, aAB, c)} or C{None} 

550 if none found. 

551 

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) 

556 

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. 

559 

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

568 

569 @return: Yield an L{Intersector5Tuple}C{(A, B, sAB, aAB, c)} for 

570 every intersection found. 

571 

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) 

595 

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

612 

613 for T, t in Tt_: 

614 if t < d or Q is _XINF: 

615 Q, d, q = T, t, i 

616 i += 1 

617 

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) 

625 

626 return None if Q is _XINF else Q.set_(sX0=d, iteration=q) 

627 

628 def _obliqDist4(self): 

629 zx = 45.0 

630 if self.f: 

631 _abs, _cjD5 = fabs, self._conjDist5 

632 

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 

650 

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 

656 

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 

682 

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 

689 

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 

697 

698 @Property_RO 

699 def R(self): 

700 '''Get the I{authalic} earth radius (C{meter}). 

701 ''' 

702 return self.ellipsoid.R2 

703 

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) 

711 

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) 

720 

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 

742 

743 @Property_RO 

744 def _T2D1Delta(self): 

745 return self._T2d3Delta(self._D1) 

746 

747 @Property_RO 

748 def _T2D2Delta(self): 

749 return self._T2d3Delta(self._D2) 

750 

751 def _T2d3Delta(self, d3): 

752 return self._T2 - d3 - self.Delta 

753 

754 @Property_RO 

755 def _Tol(self): # convergence tolerance 

756 return self._C_2 * pow(EPS, 0.75) # _0_75 

757 

758 def toStr(self, **prec_sep_name): # PYCHOK signature 

759 '''Return this C{Intersector} as string. 

760 

761 @see: L{Ellipsoid.toStr<pygeodesy.ellipsoids.Ellipsoid.toStr>} 

762 for further details. 

763 

764 @return: C{Intersector} (C{str}). 

765 ''' 

766 return self._instr(props=(Intersector.geodesic,), **prec_sep_name) 

767 

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) 

782 

783 

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

789 

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) 

799 

800 

801class _List(list): 

802 

803 _Delta = 0 # equality margin 

804 

805 def __init__(self, Delta): 

806 self._Delta = Delta 

807# list.__init__(self) 

808 

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 

817 

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 

825 

826 def sortrim(self, X0, sMax): 

827 # trim and sort the X items 

828 

829 def _key(Xk): 

830 _, k = Xk 

831 return k # rank of X 

832 

833 for X, _ in sorted(self.trim(X0, sMax), key=_key): 

834 yield X # de-tuple (X, k) 

835 

836 def trim(self, X0, sMax): 

837 # trim and yield 2-tuple (X, rank) 

838 a, b, _eu = X0.sA, X0.sB, euclid 

839 

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 

845 

846 

847if __name__ == '__main__': 

848 

849 from pygeodesy import GeodesicExact, printf 

850 

851 I = Intersector(GeodesicExact(), name='Test') # PYCHOK I 

852 

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) 

858 

859 for i, t in enumerate(I.Next5s(a, b)): 

860 printf('Next5s %s: %r (%s)', i, t, t.iteration) 

861 

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) 

867 

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 

879 

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

888 

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.