Coverage for pygeodesy/geodesici.py: 91%

622 statements  

« prev     ^ index     » next       coverage.py v7.6.0, created at 2024-07-23 11:21 -0400

1 

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

3 

4u'''Classes L{Intersectool} and L{Intersector} to find the intersections of two geodesic lines or line segments. 

5 

6Class L{Intersector} is a pure Python version of I{Karney}'s C++ class U{Intersect 

7<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Intersect.html>}. 

8 

9Class L{Intersectool} is a wrapper to invoke I{Karney}'s U{IntersectTool 

10<https://GeographicLib.SourceForge.io/C++/doc/IntersectTool.1.html>} utility, but intended I{for testing purposes only}. 

11 

12Set env variable C{PYGEODESY_INTERSECTTOOL} to the (fully qualified) path of the C{IntersectTool} executable. For usage 

13and some examples run C{"env PYGEODESY_INTERSECTTOOL=<IntersectTool-path> python3 -m pygeodesy.geodesici --help"}. 

14 

15Both L{Intersectool} and L{Intersector} provide methods C{All}, C{Closest}, C{Next} and C{Segment} and produce 

16L{XDict} instances with 4 or more items. Adjacent methods C{All5}, C{Closest5}, C{Next5} and C{Segment} return 

17or yield L{Intersectool5Tuple} or L{Intersector5Tuple}s with the lat-, longitude and azimuth of each intersection 

18as an extended, geodesic C{Position}-like L{GDict} instance. 

19 

20For more details, see the C++ U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/index.html>} 

21documentation, I{Charles F.F. Karney}'s paper U{Geodesics intersections<https://arxiv.org/abs/2308.00495>} 

22and I{S. Baselga Moreno & J.C. Martinez-Llario}'s U{Intersection and point-to-line solutions for geodesics 

23on the ellipsoid<https://riunet.UPV.ES/bitstream/handle/10251/122902/Revised_Manuscript.pdf>}. 

24''' 

25# make sure int/int division yields float quotient 

26from __future__ import division as _; del _ # PYCHOK semicolon 

27 

28from pygeodesy.basics import _copy, _enumereverse, map1, \ 

29 _xinstanceof, _xor 

30from pygeodesy.constants import EPS, INF, INT0, PI, PI2, PI_4, \ 

31 _0_0, _0_5, _1_0, _1_5, _2_0, _3_0, \ 

32 _90_0, isfinite 

33from pygeodesy.ellipsoids import _EWGS84, Fmt, unstr 

34from pygeodesy.errors import GeodesicError, IntersectionError, _an, \ 

35 _xgeodesics, _xkwds_get, _xkwds_kwds, \ 

36 _xkwds_pop2 

37# from pygeodesy.errors import exception_chaining # _MODS 

38from pygeodesy.fmath import euclid, fdot 

39from pygeodesy.fsums import Fsum, fsum1_, _ceil 

40from pygeodesy.interns import NN, _A_, _B_, _c_, _COMMASPACE_, \ 

41 _HASH_, _M_, _not_, _SPACE_, _too_ 

42from pygeodesy.karney import Caps, _diff182, GDict, _sincos2de 

43from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, \ 

44 _getenv, _PYGEODESY_INTERSECTTOOL_ 

45from pygeodesy.named import ADict, _NamedBase, _NamedTuple, _Pass 

46# from pygeodesy.namedTuples import _LL4Tuple # _MODS 

47from pygeodesy.props import deprecated_method, Property, \ 

48 Property_RO, property_RO, property_ROver 

49from pygeodesy.solveBase import _SolveCapsBase, pairs 

50# from pygeodesy.streprs import pairs # from .solveBase 

51# from pygeodesy.streprs import Fmt, unstr # from .ellipsoids 

52from pygeodesy.units import Degrees, Float, Int, _isDegrees, \ 

53 Lat, Lon, Meter, Meter_ 

54from pygeodesy.utily import sincos2, atan2, fabs, radians 

55 

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

57 

58__all__ = _ALL_LAZY.geodesici 

59__version__ = '24.07.22' 

60 

61_0t = 0, # int 

62_1_1t = -1, +1 

63_1_0_1t = -1, 0, +1 

64_aAB_ = 'aAB' 

65_c__ = '-c' # PYCHOK used! 

66_cWGS84 = _EWGS84.a * PI2 # outer circumference 

67_EPS3 = EPS * _3_0 

68_EPSr5 = pow(EPS, 0.2) # PYCHOK used! 7.4e-4 or ~3" 

69_i__ = '-i' # PYCHOK used! 

70_latA_ = 'latA' 

71_lonA_ = 'lonA' 

72_n__ = '-n' # PYCHOK used! 

73_o__ = '-o' # PYCHOK used! 

74_R__ = '-R' 

75_sAB_ = 'sAB' 

76_sX0_ = 'sX0' 

77_TRIPS = 128 

78 

79 

80class Azi(Degrees): 

81 '''(INTERNAL) Azimuth C{Unit}. 

82 ''' 

83 pass 

84 

85 

86class XDict(ADict): 

87 '''4+Item result from L{Intersectool} and L{Intersector} methods 

88 C{All}, C{Closest}, C{Next} and C{Segment} with the intersection 

89 offsets C{sA}, C{sB} and C{sX0} in C{meter} and the coincidence 

90 indicator C{c}, an C{int}, +1 for parallel, -1 for anti-parallel 

91 or 0 otherwise. 

92 

93 Offsets C{sA} and C{sB} are distances measured I{along} geodesic 

94 line C{glA} respectively C{glB}, but C{sX0} is the I{L1-distance} 

95 between the intersection and the I{origin} C{X0}. 

96 

97 If present, distance C{sAB} and angular distance C{aAB} represent 

98 the difference between the intersection point on geodesic lines 

99 C{glA} and C{glB} in C{meter} respectively C{degrees}, typically 

100 below C{5e-9 meter} or C{5 nm} and C{5e-14 degrees} or C{1 n"}. 

101 

102 For segments, indicators C{kA} and C{kB} are C{0} if the segments 

103 intersect or C{-1} or C{+1} if the intersection is I{before} the 

104 start, respectively I{after} the end of the segment, similar to 

105 L{Intersection3Tuple<Intersection3Tuple>}. Segment indicator 

106 C{k} is I{Karney}'s C{segmode}, equal C{kA * 3 + kB}. 

107 ''' 

108 _Delta = EPS # default margin, see C{Intersector._Delto} 

109 

110 def __add__(self, other): 

111 X = _copy(self) 

112 X += other 

113 return X 

114 

115 def __eq__(self, other): 

116 return not self.__ne__(other) 

117 

118 def __iadd__(self, other): 

119 if isinstance(other, tuple): # and len(other) == 2: 

120 a, b = other 

121 else: 

122 # _xinstanceof(XDict, other=other) 

123 a = other.sA 

124 b = other.sB 

125 if other.c: 

126 self.c = other.c 

127 self.sA += a # PYCHOK sA 

128 self.sB += b # PYCHOK sB 

129 return self 

130 

131 def __le__(self, other): 

132 # _xinstanceof(XDict, other=other) 

133 return self == other or self < other 

134 

135 def __lt__(self, other): 

136 # _xinstanceof(XDict, other=other) 

137 return (self.sA < other.sA or (self.sA == other.sA and # PYCHOK sA 

138 self.sB < other.sB) and self != other) # PYCHOK sB 

139 

140 def __ne__(self, other): 

141 # _xinstanceof(XDict, other=other) 

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

143 

144 def _corners(self, sA, sB, T2): 

145 # yield all corners further than C{T2} 

146 a, b = self.sA, self.sB # PYCHOK sA, sB 

147 for x in (0, sA): 

148 for y in (0, sB): 

149 if _L1(x - a, y - b) >= T2: 

150 yield XDict_(x, y) 

151 

152 def _fixCoincident(self, X, c0=0): 

153 # return the mid-point if C{X} is anti-/parallel 

154 c = c0 or X.c 

155 if c: 

156 s = (self.sA - X.sA + # PYCHOK sA 

157 (self.sB - X.sB) * c) * _0_5 # PYCHOK sB 

158 X = X + (s, s * c) # NOT += 

159 return X 

160 

161 def _fixSegment(self, sA, sB): # PYCHOK no cover 

162 # modify this anti-/parallel C{XDict} 

163 a, b, c = self.sA, self.sB, self.c # PYCHOK sA, sB, c 

164 

165 def _g(): # intersection in smallest gap 

166 if c > 0: # distance to [A, B] is |(a - b) - (A - B)| 

167 t = a - b # consider corners [0, sB] and [sA, 0] 

168 t = fabs(t + sB) < fabs(t - sA) 

169 s = a + b 

170 else: # distance to [A, B] is |(a + b) - (A + B)| 

171 t = a + b # consider corner [0, 0] and [sA, sB] 

172 t = fabs(t) < fabs(t - (sA + sB)) 

173 s = sB + (a - b) 

174 return (sB if t else sA) - s 

175 

176 ta = -a 

177 tb = sA - a 

178 tc = -c * b 

179 td = -c * (b - sB) 

180 

181 ga = 0 <= (b + c * ta) <= sB 

182 gb = 0 <= (b + c * tb) <= sB 

183 gc = 0 <= (a + tc) <= sA 

184 gd = 0 <= (a + td) <= sA 

185 

186 # test opposite rectangle sides first 

187 s = ((ta + tb) if ga and gb else ( 

188 (tc + td) if gc and gd else ( 

189 (ta + tc) if ga and gc else ( 

190 (ta + td) if ga and gd else ( 

191 (tb + tc) if gb and gc else ( 

192 (tb + td) if gb and gd else _g())))))) * _0_5 

193 self += s, s * c 

194 

195 @property_RO 

196 def _is00(self): 

197 return not (self.sA or self.sB) # PYCHOK sA, sB 

198 

199 def L1(self, other=None): 

200 '''Return the C{L1} distance. 

201 ''' 

202 a, b = self.sA, self.sB # PYCHOK sA, sB 

203 if other is not None: 

204 # _xinstanceof(XDict, other=other) 

205 a -= other.sA 

206 b -= other.sB 

207 return _L1(a, b) 

208 

209 def _nD1(self, D1): 

210 # yield the C{Closest} starts 

211 D_ = 0, D1, -D1 

212 for a, b in zip((0, 1, -1, 0, 0), 

213 (0, 0, 0, 1, -1)): 

214 yield self + (D_[a], D_[b]) 

215 

216 def _nD2(self, D2): 

217 # yield the C{Next} starts 

218 D22 = D2 * _2_0 

219 D_ = 0, D2, D22, -D22, -D2 

220 for a, b in zip((-1, -1, 1, 1, -2, 0, 2, 0), 

221 (-1, 1, -1, 1, 0, 2, 0, -2)): 

222 yield self + (D_[a], D_[b]) 

223 

224 def _nmD3(self, n, m, D3): # d3 / 2 

225 # yield the C{All} starts 

226 yield self 

227 for i in range(n, m, 2): 

228 for j in range(n, m, 2): 

229 if i or j: # skip self 

230 yield self + ((i + j) * D3, 

231 (i - j) * D3) 

232 

233 def _outSide(self, sA, sB): 

234 # is this C{Xdist} outside one or both segments? 

235 a, b = self.sA, self.sB # PYCHOK sA, sB 

236 kA = -1 if a < 0 else (+1 if a > sA else INT0) 

237 kB = -1 if b < 0 else (+1 if b > sB else INT0) 

238 self.set_(kA=kA, kB=kB, k=(kA * 3 + kB) or INT0) 

239 return bool(kA or kB) 

240 

241 def _skip(self, S_, T1_Delta): 

242 # remove starts from list C{S_} near this C{XDict} 

243 for j, S in _enumereverse(S_): 

244 if S.L1(self) < T1_Delta: 

245 S_.pop(j) 

246 

247 

248def XDict_(sA=_0_0, sB=_0_0, c=INT0, sX0=_0_0): 

249 '''(INTERNAL) New L{XDict} from positionals. 

250 ''' 

251 return XDict(sA=sA, sB=sB, c=c, sX0=sX0) 

252 

253_X000 = XDict_() # PYCHOK origin 

254_XINF = XDict_(INF) 

255 

256 

257class _IntersectBase(_NamedBase): 

258 '''(INTERNAL) Base class for L{Intersectool} and L{Intersector}. 

259 ''' 

260 # _g = None 

261 

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

263 _xinstanceof(*_EWGS84._Geodesics, geodesic=geodesic) 

264 self._g = geodesic 

265 if name: 

266 self.name = name 

267 

268 @Property_RO 

269 def a(self): 

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

271 ''' 

272 return self.ellipsoid.a 

273 

274 equatoradius = a # = Requatorial 

275 

276 def All(self, glA, glB, **kwds): # PYCHOK no cover 

277 '''(INTERNAL) I{Must be overloaded}.''' 

278 self._notOverloaded(glA, glB, **kwds) 

279 

280 @Property_RO 

281 def _cHalf(self): # normalizer, semi-circumference 

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

283 

284 @Property_RO 

285 def _cMax(self): # outer circumference 

286 return max(self.a, self.ellipsoid.b, self.R) * PI2 

287 

288 @property_RO 

289 def datum(self): 

290 '''Get the geodesic's datum (C{Datum}). 

291 ''' 

292 return self.geodesic.datum 

293 

294 @Property_RO 

295 def ellipsoid(self): 

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

297 ''' 

298 return self.geodesic.datum.ellipsoid 

299 

300 @Property_RO 

301 def f(self): 

302 '''Get the I{flattening} (C{scalar}), C{0} for spherical, negative for prolate. 

303 ''' 

304 return self.ellipsoid.f 

305 

306 flattening = f 

307 

308 @property_RO 

309 def geodesic(self): 

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

311 ''' 

312 return self._g 

313 

314 def _illz2G(self, G, il): 

315 '''(INTERNAL) Set C{InverseLine} 1-/2-attrs into C{G}, a C{GDict}. 

316 ''' 

317 try: 

318 G.set_(lat1=il.lat1, lon1=il.lon1, azi1=il.azi1, a12=il.a13, # .Arc() 

319 lat2=il.lat2, lon2=il.lon2, azi2=il.azi2, s12=il.s13) # .Distance() 

320 except AttributeError: 

321 r = il.Position(il.s13, outmask=Caps._STD_LINE) # isfinite(il.s13) 

322 G.set_(**r) 

323# for n, v in r.items(): 

324# if not hasattr(il, n): 

325# setattr(il, n, v) 

326 return G 

327 

328 def intersect7(self, start1, end1, start2, end2, X0=_X000, aMaX0=0, sMaX0=_cWGS84, 

329 **LatLon_and_kwds): 

330 '''Yield the intersection points of two lines, each defined by two (ellipsoidal) 

331 points or by an (ellipsoidal) start point and an azimuth from North. 

332 

333 @arg start1: Start point of the first line (C{LatLon}). 

334 @arg end1: End point of the first line (C{LatLon}) or the azimuth at the 

335 B{C{start1}} point (compass C{degrees360}). 

336 @arg start2: Start point of the second line (C{LatLon}). 

337 @arg end2: End point of the second line (C{LatLon}) or the azimuth at the 

338 B{C{start2}} point (compass C{degrees360}). 

339 @kwarg X0: Optional I{origin} for I{L1-distances} (L{XDict}) or C{None} for 

340 the L{Middle<Intersector.Middle>}, otherwise C{XDiff_(0, 0)}. 

341 @kwarg aMaX0: Upper limit for the I{angular L1-distance} 

342 (C{degrees}) or C{None} or C{0} for unlimited. 

343 @kwarg sMaX0_C: Optional, upper limit C{B{sMaX0}=2*PI*R} for the 

344 I{L1-distance} to B{C{X0}} (C{meter}). 

345 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=None} to return intersection 

346 points and optional, additional B{C{LatLon}} keyword arguments. 

347 

348 @note: The C{lat} and C{lon} attr of B{C{start1}}, B{C{end1}}, B{C{start2}} and 

349 B{C{end2}} are used I{verbatim}, ignoring C{datum} or C{ellipsoid}. 

350 

351 @return: Yield an L{Intersect7Tuple}C{(A, B, sAB, aAB, c, kA, kB)} for every 

352 intersection found, with C{A} and C{B} each a B{C{LatLon}} or if 

353 C{B{LatLon} is None} or not specified, a L{LatLon4Tuple}C{(lat, lon, 

354 height, datum)} with C{height 0} and this C{datum}. 

355 

356 @raise GeodesicError: Invalid B{C{start1}}, B{C{end1}}, B{C{start2}} or 

357 B{C{end2}} or B{C{end1}} and B{C{end2}} differ in type. 

358 

359 @raise IntersectionError: No convergence. 

360 ''' 

361 

362 def _args(s, e): 

363 t = (e,) if _isDegrees(e) else (e.lat, e.lon) 

364 return (s.lat, s.lon) + t 

365 

366 try: 

367 glA = self.Line(*_args(start1, end1)) 

368 glB = self.Line(*_args(start2, end2)) 

369 except Exception as x: 

370 raise GeodesicError(start1=start1, end1=end1, start2=start2, end2=end2, cause=x) 

371 

372 LL, kwds = _xkwds_pop2(LatLon_and_kwds, LatLon=None) 

373 d, kwds = _xkwds_pop2(kwds, datum=self.datum) 

374 h, kwds = _xkwds_pop2(kwds, height=0) 

375 

376 _LL4T = _MODS.namedTuples._LL4Tuple 

377 for X in self.All(glA, glB, X0=X0, aMaX0=aMaX0, sMaX0=sMaX0, _C=True): 

378 A = B = _LL4T(X.latA, X.lonA, h, d, LL, kwds, iteration=X.iteration) 

379 if X.sAB or X.latA != X.latB or X.lonA != X.lonB: 

380 B = _LL4T(X.latB, X.lonB, h, d, LL, kwds, iteration=X.iteration) 

381 yield Intersect7Tuple(A, B, X.sAB, X.aAB, X.c, _xkwds_get(X, kA=0), 

382 _xkwds_get(X, kB=0)) 

383 

384 def _Inversa12(self, A, B=None): 

385 lls = (0, 0, A, 0) if B is None else (A.lat2, A.lon2, 

386 B.lat2, B.lon2) 

387 r = self._g.Inverse(*lls, outmask=Caps.DISTANCE) 

388 return r.s12, r.a12 # .a12 always in r 

389 

390 def k2kAkB(self, k): 

391 '''Unravel C{k} into C{kA} and C{kB}. 

392 

393 @arg k: Segment indicator C{kA * 3 + kB} (C{int}). 

394 

395 @return: An C{ADict(k=k, kA=kA, kB=kB)}. 

396 

397 @raise GeodesicError: Invalid B{C{k}}. 

398 ''' 

399 for kA in range(-1, 2): 

400 for kB in range(-1, 2): 

401 if (kA * 3 + kB) == k: 

402 return ADict(k=k, kA=kA, kB=kB) 

403 raise GeodesicError(k=k) 

404 

405# def k2kAkB(self, k): 

406# # unravel C{k} into C{kA} and C{kB}. 

407# kA, kB = divmod(k, 3) 

408# if kB > 1: 

409# kA += 1 

410# kB -= 3 

411# return kA, kB 

412 

413 def Line(self, lat1, lon1, azi1_lat2, *lon2, **name): # PYCHOK no cover 

414 '''(INTERNAL) I{Must be overloaded}.''' 

415 self._notOverloaded(lat1, lon1, azi1_lat2, *lon2, **name) 

416 

417 def _ll3z4ll(self, lat1, lon1, azi1_lat2, *lon2): 

418 t = Lat(lat1=lat1), Lon(lon1=lon1) 

419 if lon2: # get azis for All, keep lat-/lons 

420 t += Lat(lat2=azi1_lat2), Lon(lon2=lon2[0]) 

421 else: 

422 t += Azi(azi1=azi1_lat2), 

423 return t 

424 

425 @deprecated_method 

426 def Next5s(self, glA, glB, X0=_X000, aMax=1801, sMax=0, **unused): # PYCHOK no cover 

427 '''DEPRECATED on 2024.07.02, use method C{All5}.''' 

428 return self.All5(glA, glB, X0=X0, aMaX0=aMax, sMaX0=sMax) # PYCHOK attr 

429 

430 @Property_RO 

431 def R(self): 

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

433 ''' 

434 return self.ellipsoid.R2 

435 

436 def _sMaX0_C2(self, aMaX0=0, **sMaX0_C): 

437 _g = _xkwds_get 

438 s = _g(sMaX0_C, sMaX0=self._cMax) 

439 s = _g(sMaX0_C, sMax=s) # for backward ... 

440 a = _g(sMaX0_C, aMax=aMaX0) # ... compatibility 

441 if a: # degrees to meter, approx. 

442 s = min(s, self.R * radians(a)) # ellipsoid.degrees2m(a) 

443 s = _g(sMaX0_C, _R=s) 

444 if s < _EPS3: 

445 s = _EPS3 # raise GeodesicError(sMaX0=s) 

446 return s, _g(sMaX0_C, _C=False) 

447 

448 def _xNext(self, glA, glB, eps1, **eps_C): # PYCHOK no cover 

449 eps1 = _xkwds_get(eps_C, eps=eps1) # eps for backward compatibility 

450 if eps1 is not None: 

451 a = glA.lat1 - glB.lat1 

452 b = glA.lon1 - glB.lon1 

453 if euclid(a, b) > eps1: 

454 raise GeodesicError(lat_=a, lon_=b, eps1=eps1) 

455 return _xkwds_kwds(eps_C, _C=False) 

456 

457 

458class Intersectool(_IntersectBase, _SolveCapsBase): 

459 '''Wrapper to invoke I{Karney}'s utility U{IntersectTool 

460 <https://GeographicLib.SourceForge.io/C++/doc/IntersectTool.1.html>} 

461 similar to class L{Intersector<geodesici.Intersector>}. 

462 

463 @note: Use property C{IntersectTool} or env variable C{PYGEODESY_INTERSECTTOOL} 

464 to specify the (fully qualified) path to the C{IntersectTool} executable. 

465 

466 @note: This C{Intersectool} is intended I{for testing purposes only}, it invokes 

467 the C{IntersectTool} executable for I{every} method call. 

468 ''' 

469 _c_alt = _c__, # Closest latA lonA aziA latB lonB aziB 

470 _C_option = '-C', 

471 _Error = GeodesicError 

472 _i_alt = _i__, # Segment latA1 lonA1 latA2 lonA2 latB1 lonB1 latB2 lonB2 

473 _linelimit = 1200 # line printer width X 10 

474 _n_alt = _n__, # Next latA lonA aziA aziB 

475 _Names_ABs = _latA_, _lonA_, 'latB', 'lonB', _sAB_ # -C to stderr 

476 _Names_XDict = 'sA', 'sB', _c_ # plus 'k' from -i or 'sX0' from -R 

477 _o_alt = _o__, # Offset latA lonA aziA latB lonB aziB x0 y0 

478 _Xable_name = 'IntersectTool' 

479 _Xable_path = _getenv(_PYGEODESY_INTERSECTTOOL_, _PYGEODESY_INTERSECTTOOL_) 

480 

481 def __init__(self, a_geodesic=None, f=None, **name): 

482 '''New L{IntersectTool}. 

483 

484 @arg a_geodesic: Earth' equatorial axis (C{meter}) or a geodesic 

485 (L{GeodesicExact<pygeodesy.geodesicx.GeodesicExact>}, 

486 wrapped L{Geodesic<pygeodesy.geodesicw.Geodesic>} or 

487 L{GeodesicSolve<pygeodesy.geodsolve.GeodesicSolve>}). 

488 @kwarg f: Earth' flattening (C{scalar}), required if B{C{a_geodesic}} 

489 is in C{meter}, ignored otherwise. 

490 @kwarg name: Optional C{B{name}=NN} (C{str}). 

491 

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

493 large or no initial convergence. 

494 

495 @see: The B{Note} at I{Karney}'s C++ U{Intersect<https://GeographicLib.sourceforge.io/ 

496 C++/doc/classGeographicLib_1_1Intersect.html#ae41f54c9a44836f6c8f140f6994930cf>}. 

497 ''' 

498 g = self._GeodesicExact() if a_geodesic is None else (a_geodesic if f is None else 

499 self._GeodesicExact(a_geodesic, f)) 

500 _IntersectBase.__init__(self, g, **name) 

501 

502 def All(self, glA, glB, X0=_X000, eps1=_0_0, aMaX0=0, **sMaX0_C): # PYCHOK signature 

503 '''Yield all intersection of two geodesic lines up to a limit. 

504 

505 @kwarg eps1: Optional margin for the L{euclid<pygeodesy.euclid>}ean distance 

506 (C{degrees}) between the C{(lat1, lon1)} points of both lines for 

507 using the L{IntersectTool<Intersectool.IntersectTool>}'s C{"-n"} 

508 option, unless C{B{eps1}=None}. 

509 

510 @return: An L{XDict} for each intersection. 

511 ''' 

512 for X, _ in self._All2(glA, glB, X0, eps1, aMaX0=aMaX0, **sMaX0_C): 

513 yield X 

514 

515 def _All2(self, glA, glB, X0, eps1, **aMaX0_sMaX0_C): # MCCABE 13 

516 '''(INTERNAL) Helper for methods C{.All} and C{.All5}. 

517 ''' 

518 def _xz2(**gl): 

519 try: 

520 n, gl = gl.popitem() # _xkwds_item2(gl) 

521 try: 

522 return self._c_alt, (gl.azi1,) 

523 except (AttributeError, KeyError): 

524 return self._i_alt, (gl.lat2, gl.lon2) 

525 except Exception as x: 

526 raise GeodesicError(n, gl, cause=x) 

527 

528 _t, a = _xz2(glA=glA) 

529 _x, b = _xz2(glB=glB) 

530 if _x is not _t: 

531 raise GeodesicError(glA=glA, glB=glB) 

532 

533 A = glA.lat1, glA.lon1 

534 B = glB.lat1, glB.lon1 

535 if _x is self._c_alt: 

536 if X0 is _X000 or X0._is00: 

537 if eps1 is not None and \ 

538 euclid(glA.lat1 - glB.lat1, 

539 glA.lon1 - glB.lon1) <= eps1: 

540 _x, B = self._n_alt, () 

541 else: # non-zero offset 

542 _x = self._o_alt 

543 b += X0.sA, X0.sB 

544 

545 sMaX0, _C = self._sMaX0_C2(**aMaX0_sMaX0_C) 

546 for X in self._XDictInvoke(_x, _sX0_, (A + a + B + b), 

547 _C=_C, _R=sMaX0): 

548 if _C: 

549 T = self._In5T(glA, glB, X, X) 

550 if _aAB_ not in X: 

551 X.set_(sAB=T.sAB, aAB=T.aAB) 

552 else: 

553 T = None 

554 yield X.set_(c=int(X.c)), T 

555 

556 def All5(self, glA, glB, X0=_X000, **aMaX0_sMaX0): 

557 '''Yield all intersection of two geodesic lines up to a limit. 

558 

559 @return: An L{Intersectool5Tuple} for each intersection. 

560 ''' 

561 for _, T in self._All2(glA, glB, X0, _0_0, _C=True, **aMaX0_sMaX0): 

562 yield T 

563 

564 @Property_RO 

565 def _cmdBasic(self): 

566 '''(INTERNAL) Get the basic C{IntersectTool} cmd (C{tuple}). 

567 ''' 

568 return (self.IntersectTool,) + (self._e_option + 

569 self._E_option + 

570 self._p_option) 

571 

572 def Closest(self, glA, glB, X0=_X000, _C=False): 

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

574 

575 @kwarg _C: Use C{B{_C}=True} to include the C{"-C"} results (C{bool}). 

576 

577 @return: An L{XDict}. 

578 ''' 

579 args = glA.lat1, glA.lon1, glA.azi1, \ 

580 glB.lat1, glB.lon1, glB.azi1 

581 if X0 is _X000 or X0._is000: 

582 _x = self._c_alt 

583 else: 

584 _x = self._o_alt 

585 args += X0.sA, X0.sB 

586 return self._XDictInvoke(_x, NN, args, _C=_C) # _R=None) 

587 

588 def Closest5(self, glA, glB, **unused): 

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

590 

591 @return: An L{Intersectool5Tuple}. 

592 ''' 

593 X = self.Closest(glA, glB, _C=True) 

594 return self._In5T(glA, glB, X, X) 

595 

596 @property_ROver 

597 def _GeodesicExact(self): 

598 '''Get the I{class} L{GeodesicExact}, I{once}. 

599 ''' 

600 return _MODS.geodesicx.GeodesicExact # overwrite propertyROver 

601 

602 def _In5T(self, glA, glB, S, X, k2=False, **_2X): 

603 A = GDict(glA).set_(lat2=X.latA, lon2=X.lonA, s12=S.sA) 

604 B = GDict(glB).set_(lat2=X.latB, lon2=X.lonB, s12=S.sB) 

605 if k2: 

606 A.set_(k2=X.kA) 

607 B.set_(k2=X.kB) 

608 s, a = self._Inversa12(A, B) 

609 sAB = _xkwds_get(X, sAB=s) 

610 if a and s and s != sAB: 

611 a *= sAB / s # adjust a 

612 return Intersectool5Tuple(A._2X(glA, **_2X), 

613 B._2X(glB, **_2X), sAB, a, X.c) 

614 

615 @Property 

616 def IntersectTool(self): 

617 '''Get the U{IntersectTool<https://GeographicLib.SourceForge.io/C++/doc/IntersectTool.1.html>} 

618 executable (C{filename}). 

619 ''' 

620 return self._Xable_path 

621 

622 @IntersectTool.setter # PYCHOK setter! 

623 def IntersectTool(self, path): 

624 '''Set the U{IntersectTool<https://GeographicLib.SourceForge.io/C++/doc/IntersectTool.1.html>} 

625 executable (C{filename}), the (fully qualified) path to the C{IntersectTool} executable. 

626 

627 @raise GeodesicError: Invalid B{C{path}}, B{C{path}} doesn't exist or isn't the 

628 C{IntersectTool} executable. 

629 ''' 

630 self._setXable(path) 

631 

632 def Line(self, lat1, lon1, azi1_lat2, *lon2, **name): 

633 '''Return a geodesic line from this C{Intersector}'s geodesic, specified by 

634 two (goedetic) points or a (goedetic) point and an (forward) azimuth. 

635 

636 @return: A 3- or 6-item, named L{GDict}. 

637 ''' 

638 args = self._ll3z4ll(lat1, lon1, azi1_lat2, *lon2) 

639 gl = GDict((u.name, u) for u in args) 

640# if lon2: # get azis for All, use lat-/lons as given 

641# r = self._g.Inverse(outmask=Caps.AZIMUTH, *args) 

642# gl.set_(azi1=Azi(azi1=r.azi1), azi2=Azi(azi2=r.azi2)) 

643 if name: 

644 gl.name= name 

645 return gl 

646 

647 def Middle(self, glA, glB, **_C): 

648 '''Get the mid-points on two geodesic line segments. 

649 

650 @kwarg _C: Use C{B{_C}=True} to include the C{"-C"} results (C{bool}). 

651 

652 @return: An L{XDict}. 

653 ''' 

654 X, _, _, _, _ = self._middle5(glA, glB, **_C) 

655 return X 

656 

657 def _middle5(self, glA, glB, _C=False, **unused): 

658 # return intersections C{A} and C{B} and the 

659 # center C{X0} of rectangle [sA, sB] 

660 

661 def _smi4(**gl): 

662 try: 

663 n, gl = gl.popitem() 

664 il = self._g.InverseLine(gl.lat1, gl.lon1, gl.lat2, gl.lon2) 

665 except Exception as x: 

666 raise GeodesicError(n, gl, cause=x) 

667 s = il.s13 

668 m = s * _0_5 

669 return s, m, il, (il.Position(m, outmask=Caps._STD_LINE) if _C else None) 

670 

671 sA, mA, iA, A = _smi4(glA=glA) 

672 sB, mB, iB, B = _smi4(glB=glB) 

673 X = XDict_(mA, mB) # centers 

674 _ = X._outSide(sA, sB) 

675 if _C: # _Names_ABs 

676 s, a = self._Inversa12(A, B) 

677 X.set_(latA=A.lat2, lonA=A.lon2, aMM=a, # assert sA == A.s12 

678 latB=B.lat2, lonB=B.lon2, sMM=s) # assert sB == B.s12 

679 return X, A, iA, B, iB 

680 

681 def Middle5(self, glA, glB, **unused): 

682 '''Get the mid-points on two geodesic line segments and their distance. 

683 

684 @return: A L{Middle5Tuple}. 

685 ''' 

686 X, A, iA, B, iB = self._middle5(glA, glB, _C=True) 

687 A, B, s, a, c = self._In5T(A, B, X, X, _2X=_M_) 

688 return Middle5Tuple(self._illz2G(A, iA), 

689 self._illz2G(B, iB), s, a, c) 

690 

691 def Next(self, glA, glB, eps1=None, **_C): # PYCHOK no cover 

692 '''Find the next intersection of two I{intersecting} geodesic lines. 

693 

694 @kwarg _C: Use C{B{_C}=True} to include the option C{"-C"} results (C{bool}). 

695 

696 @return: An L{XDict}. 

697 ''' 

698 if eps1 or _C: 

699 _C = self._xNext(glA, glB, eps1, **_C) 

700 return self._XDictInvoke(self._n_alt, NN, 

701 (glA.lat1, glA.lon1, glA.azi1, glB.azi1), 

702 **_C) # _R=None 

703 

704 def Next5(self, glA, glB, **eps1): # PYCHOK no cover 

705 '''Find the next intersection of two I{intersecting} geodesic lines. 

706 

707 @return: An L{Intersectool5Tuple}. 

708 ''' 

709 X = self.Next(glA, glB, _C=True, **eps1) 

710 return self._In5T(glA, glB, X, X) 

711 

712 def _R_option(self, _R=None): 

713 '''(INTERNAL) Get the C{-R maxdist} option. 

714 ''' 

715 return () if _R is None else (_R__, str(_R)) # -R maxdist 

716 

717 def Segment(self, glA, glB, **_C_unused): 

718 '''Find the intersection between two geodesic line segments. 

719 

720 @kwarg _C: Use C{B{_C}=True} to include the option C{"-C"} results (C{bool}). 

721 

722 @return: An L{XDict}. 

723 ''' 

724 X = self._XDictInvoke(self._i_alt, 'k', 

725 (glA.lat1, glA.lon1, glA.lat2, glA.lon2, 

726 glB.lat1, glB.lon1, glB.lat2, glB.lon2), 

727 _C=_xkwds_get(_C_unused, _C=False)) # _R=None 

728 try: 

729 ks = self.k2kAkB(int(X.k)) 

730 except Exception as x: 

731 raise GeodesicError(glA=glA, glB=glB, X=str(X), cause=x) 

732 return X.set_(**ks) 

733 

734 def Segment5(self, glA, glB, **unused): 

735 '''Find the next intersection of two I{intersecting} geodesic lines. 

736 

737 @return: An L{Intersectool5Tuple}. 

738 ''' 

739 X = self.Segment(glA, glB, _C=True) 

740 return self._In5T(glA, glB, X, X, k2=True) 

741 

742 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature 

743 '''Return this C{Intersectool} as string. 

744 

745 @kwarg prec_sep: Keyword argumens C{B{prec}=6} and C{B{sep}=", "} 

746 for the C{float} C{prec}ision, number of decimal digits 

747 (0..9) and the C{sep}arator string to join. Trailing 

748 zero decimals are stripped for B{C{prec}} values of 1 

749 and above, but kept for negative B{C{prec}} values. 

750 

751 @return: Intersectool items (C{str}). 

752 ''' 

753 d = dict(geodesic=self.geodesic, invokation=self.invokation, 

754 status=self.status, 

755 IntersectTool=self.IntersectTool) 

756 return sep.join(pairs(d, prec=prec)) 

757 

758 def _XDictInvoke(self, alt, _k_sX0, args, _C=False, **_R): 

759 '''(INTERNAL) Invoke C{IntersectTool}, return results as C{XDict} or 

760 a C{generator} if keyword argument C{B{_R}=sMaX0} is specified. 

761 ''' 

762 # assert len(args) == {self._c_alt: 6, 

763 # self._i_alt: 8, 

764 # self._n_alt: 4, 

765 # self._o_alt: 8}.get(alt, len(args)) 

766 cmd = self._cmdBasic 

767 Names = self._Names_XDict # has _c_ always 

768 if _k_sX0: 

769 Names += _k_sX0, 

770 if _C: 

771 cmd += self._C_option 

772 Names += self._Names_ABs 

773 if _R: 

774 cmd += self._R_option(**_R) 

775 X, _R = self._DictInvoke2(cmd + alt, args, Names, XDict, **_R) 

776 return X if _R else X.set_(c=int(X.c)) # generator or XDict 

777 

778 

779class Intersector(_IntersectBase): 

780 '''Finder of intersections between two goedesic lines, each an instance 

781 of L{GeodesicLineExact<pygeodesy.geodesicx.GeodesicLineExact>}, 

782 wrapped L{GeodesicLine<pygeodesy.geodesicw.GeodesicLine>} or 

783 L{GeodesicLineSolve<pygeodesy.geodsolve.GeodesicLineSolve>}. 

784 

785 @see: I{Karney}'s C++ class U{Intersect<https://GeographicLib.sourceforge.io/ 

786 C++/doc/classGeographicLib_1_1Intersect.html#details>} for more details. 

787 ''' 

788 

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

790 '''New L{Intersector}. 

791 

792 @arg geodesic: The geodesic (L{GeodesicExact<pygeodesy.geodesicx.GeodesicExact>}, 

793 wrapped L{Geodesic<pygeodesy.geodesicw.Geodesic>} or 

794 L{GeodesicSolve<pygeodesy.geodsolve.GeodesicSolve>}). 

795 @kwarg name: Optional C{B{name}=NN} (C{str}). 

796 

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

798 large or no initial convergence. 

799 

800 @see: The B{Note} at I{Karney}'s C++ U{Intersect<https://GeographicLib.sourceforge.io/ 

801 C++/doc/classGeographicLib_1_1Intersect.html#ae41f54c9a44836f6c8f140f6994930cf>}. 

802 ''' 

803 _IntersectBase.__init__(self, geodesic, **name) 

804 E = self.ellipsoid 

805 t1 = E.b * PI # min distance between intersects 

806 t2 = self._polarDist2(_90_0)[0] * _2_0 # furthest, closest intersect 

807 t5 = self._Inversa12( _90_0)[0] * _2_0 # longest, shortest geodesic 

808 if self.f > 0: 

809 t3 = self._obliqDist4()[0] 

810 t4 = t1 

811 else: # PYCHOK no cover 

812 t1, t2, t3 = t2, t1, t5 

813 t4, _, _ = self._polarB3() 

814 

815 self._D1 = d1 = t2 * _0_5 # ~E.L tile spacing for Closest 

816 self._D2 = d2 = t3 / _1_5 # tile spacing for Next 

817 self._D3 = d3 = t4 - self.Delta # tile spacing for All 

818 self._T1 = t1 # min distance between intersects 

819 self._T2 = t2 = t1 * _2_0 

820# self._T5 = t5 # not used 

821 if not (d1 < d3 and d2 < d3 and d2 < t2): 

822 t = Fmt.PARENSPACED(_too_('eccentric'), E.e) 

823 raise GeodesicError(ellipsoid=E.toStr(terse=2), txt=t) 

824 

825 def All(self, glA, glB, X0=None, aMaX0=0, **sMaX0_C): # MCCABE 13 

826 '''Yield all intersection of two geodesic lines up to a limit. 

827 

828 @arg glA: A geodesic line (L{Line<Intersector.Line>}). 

829 @arg glB: An other geodesic line (L{Line<Intersector.Line>}). 

830 @kwarg X0: Optional I{origin} for I{L1-distances} (L{XDict}) or 

831 C{None} for the L{Middle<Intersector.Middle>} of both 

832 lines if both are a 4-C{args} L{Line<Intersector.Line>} 

833 or C{InverseLine}, otherwise C{XDiff_(0, 0)}. 

834 @kwarg aMaX0: Upper limit for the I{angular L1-distance} 

835 (C{degrees}) or C{None} or C{0} for unlimited. 

836 @kwarg sMaX0_C: Optional, upper limit C{B{sMaX0}=2*PI*R} for the 

837 I{L1-distance} to B{C{X0}} (C{meter}) and option 

838 C{B{_C}=False} to include the intersection lat-/ 

839 longitudes C{latA}, C{lonA}, C{latB}, C{lonB} and 

840 distances C{sAB} and C{aSB}. 

841 

842 @return: Yield an L{XDict} for each intersection found. 

843 

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

845 invalid, incompatible or ill-configured. 

846 

847 @raise IntersectionError: No convergence. 

848 ''' 

849 self._xLines(glA, glB) 

850 if X0 is None: 

851 try: # determine X0 

852 X0, _, _ = self._middle3(glA, glB, True) 

853 except GeodesicError: # no .Distance 

854 X0 = _X000 

855 sMaX0, _C = self._sMaX0_C2(aMaX0, **sMaX0_C) 

856 

857 D, _D = self.Delta, self._cHalf # C++ _d 

858 xMaX0 = sMaX0 + D 

859 m = int(_ceil(xMaX0 / self._D3)) # m x m tiles 

860 d3 = xMaX0 / m 

861 T2d3D = self._T2d3Delta(d3) 

862 

863 C_ = _List(D) # closest coincident 

864 X_ = _List(D) # intersections found 

865 c0 = 0 

866 S_ = list(X0._nmD3(1 - m, m, d3 * _0_5)) 

867 # assert len(S_) == m * m + (m - 1) % 2 

868 while S_: 

869 Q, i = self._Basic2(glA, glB, S_.pop(0)) 

870 if Q in X_: 

871 continue 

872 if Q.c: # coincident intersection # PYCHOK no cover 

873 _X0fx = X0._fixCoincident 

874 Q = _X0fx(Q) # Q = Q' 

875 if c0 and Q in C_: 

876 continue 

877 C_.addend(Q) 

878 # elimate all existing intersections 

879 # on this line (which didn't set c0) 

880 c0 = Q.c 

881 for j, X in _enumereverse(X_): 

882 if _X0fx(X, c0).L1(Q) <= D: # X' == Q 

883 X_.pop(j) 

884 

885 a, s0 = len(X_), Q.sA 

886 args = self._m12_M12_M21(glA, s0) 

887 _cjD = self._conjDist 

888 for s in (-_D, _D): 

889 s += s0 

890 sa = 0 

891 while True: 

892 i += 1 

893 sa = _cjD(glA, s + sa, *args) - s0 

894 X = Q + (sa, sa * c0) 

895 if X_.addend(X, X0.L1(X), i) > xMaX0: 

896 break 

897 

898 elif c0 and Q in C_: # Q.c == 0 

899 continue 

900 else: 

901 a = len(X_) 

902 

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

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

905 X._skip(S_, T2d3D) 

906 

907 return X_.sorter(sMaX0, self._C, glA, glB, _C=_C) # generator 

908 

909 def All5(self, glA, glB, X0=_X000, **aMaX0_sMaX0_C): 

910 '''Yield all intersection of two geodesic lines up to a limit. 

911 

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

913 for each intersection found. 

914 

915 @see: Methods L{All} for further details. 

916 ''' 

917 for X in self.All(glA, glB, X0=X0, **aMaX0_sMaX0_C): 

918 yield self._In5T(glA, glB, X, X) 

919 

920 def _Basic2(self, glA, glB, S, i=0): 

921 '''(INTERNAL) Get a basic solution. 

922 ''' 

923 X = _copy(S) 

924 for _ in range(_TRIPS): 

925 S = self._Spherical(glA, glB, X) 

926 X += S 

927 i += 1 

928 if X.c or S.L1() <= self._Tol: # or isnan 

929 return self._Delto(X), i 

930 

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

932 

933 def _C(self, X, glA, glB, _C=False, _MM=False): 

934 # add the C{_C} items to C{X}, if requested. 

935 if _C: 

936 A = self._Position(glA, X.sA) 

937 B = self._Position(glB, X.sB) 

938 s, a = self._Inversa12(A, B) 

939 X.set_(latA=A.lat2, lonA=A.lon2, 

940 latB=B.lat2, lonB=B.lon2) 

941 if _MM: # in .Middle5 

942 X.set_(sMM=s, aMM=a) 

943 else: 

944 X.set_(sAB=s, aAB=a) 

945 return X 

946 

947 def Closest(self, glA, glB, X0=_X000, **_C): 

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

949 

950 @arg glA: A geodesic line (L{Line<Intersector.Line>}). 

951 @arg glB: An other geodesic line (L{Line<Intersector.Line>}). 

952 @kwarg X0: Optional I{origin} for I{L1-closeness} (L{XDict}). 

953 @kwarg _C: If C{True}, include the lat-/longitudes C{latA}, 

954 C{lonA}, C{latB}, C{lonB} oon and distances C{sAB} 

955 and C{aSB} between the intersections. 

956 

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

958 

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

960 invalid, incompatible or ill-configured. 

961 

962 @raise IntersectionError: No convergence. 

963 ''' 

964 self._xLines(glA, glB) 

965 Q, d, S_, i = X0, INF, list(X0._nD1(self._D1)), 0 

966 while S_: 

967 X, i = self._Basic2(glA, glB, S_.pop(0), i) 

968 X = X0._fixCoincident(X) 

969 if X.L1(Q) > self.Delta: # X != Q 

970 d0 = X.L1(X0) 

971 if d0 < self._T1: 

972 Q, d, q = X, d0, i 

973 break 

974 if d0 < d or Q is X0: 

975 Q, d, q = X, d0, i 

976 X._skip(S_, self._T2D1Delta) 

977 

978 return None if Q is X0 else self._C(Q, glA, glB, **_C).set_(sX0=d, iteration=q) 

979 

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

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

982 

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

984 or C{None} if none found. 

985 

986 @see: Method L{Closest} for further details. 

987 ''' 

988 X = self.Closest(glA, glB, X0=X0) 

989 return X if X is None else self._In5T(glA, glB, X, X) 

990 

991 def _conjDist(self, gl, s, m12=0, M12=1, M21=1, semi=False): 

992 # Find semi-/conjugate point relative to s0 which is close to s1. 

993 # if semi: 

994 # solve for M23 = 0 using dM23 / ds3 = - (1 - M23 * M32) / m23 

995 # else: 

996 # solve for m23 = 0 using dm23 / ds3 = M32 

997 _S2, _abs, _1 = Fsum(s).fsum2_, fabs, _1_0 

998 for _ in range(_TRIPS): 

999 m13, M13, M31 = self._m12_M12_M21(gl, s) 

1000 # see "Algorithms for geodesics", eqs. 31, 32, 33. 

1001 m23 = m13 * M12 

1002 M32 = M31 * M12 

1003 if m12: # PYCHOK no cover 

1004 m23 -= m12 * M13 

1005 if m13: 

1006 M32 += (_1 - M13 * M31) * m12 / m13 

1007 if semi: 

1008 M23 = M13 * M21 

1009 # when m12 -> eps, (1 - M12 * M21) -> eps^2, I suppose. 

1010 if m12 and m13: 

1011 M23 += (_1 - M12 * M21) * m13 / m12 

1012 d = m23 * M23 / (_1 - M23 * M32) 

1013 else: 

1014 d = -m23 / M32 

1015 s, d = _S2(d) 

1016 if _abs(d) <= self._Tol: 

1017 break 

1018 return s 

1019 

1020 _gl3 = None 

1021 

1022 @Property 

1023 def _conjDist3s(self): 

1024 gl, self._gl3, _D = self._gl3, None, self._cHalf 

1025 return tuple(self._conjDist(gl, s) for s in (-_D, 0, _D)) 

1026 

1027 @_conjDist3s.setter # PYCHOK setter! 

1028 def _conjDist3(self, gl): 

1029 # _XLines(gl, gl) 

1030 self._gl3 = gl 

1031 

1032 def _conjDist3Tt_(self, c, X0=_X000): 

1033 for s in self._conjDist3s: 

1034 T = XDict_(s, s * c, c) 

1035 yield self._Delto(T), T.L1(X0) 

1036 

1037 def _conjDist5(self, azi): 

1038 gl = self._Line(azi1=azi) 

1039 s = self._conjDist(gl, self._cHalf) 

1040 X, _ = self._Basic2(gl, gl, XDict_(s * _0_5, -s * _1_5)) 

1041 return s, (X.L1() - s * _2_0), azi, X.sA, X.sB 

1042 

1043 @Property_RO 

1044 def Delta(self): 

1045 '''Get the equality and tiling margin (C{meter}). 

1046 ''' 

1047 return self._cHalf * _EPSr5 # ~15 Km WGS84 

1048 

1049 def _Delto(self, X): 

1050 # copy Delta into X, overriding X's default 

1051 X._Delta = self.Delta # NOT X.set_(self.Delta) 

1052 return X 

1053 

1054 @Property_RO 

1055 def _EPS3R(self): 

1056 return _EPS3 * self.R 

1057 

1058 @Property_RO 

1059 def _faPI_4(self): 

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

1061 

1062 @Property_RO 

1063 def _GeodesicLines(self): 

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

1065 ''' 

1066 return type(self._Line()), 

1067 

1068 def _In5T(self, glA, glB, S, X, k2=False, **_2X): 

1069 # Return an intersection as C{Intersector5Tuple}. 

1070 A = self._Position(glA, S.sA) 

1071 B = self._Position(glB, S.sB) 

1072 if k2: 

1073 A.set_(k2=X.kA) 

1074 B.set_(k2=X.kB) 

1075 s, a = self._Inversa12(A, B) 

1076 return Intersector5Tuple(A._2X(glA, **_2X), 

1077 B._2X(glB, **_2X), s, a, X.c, iteration=X.iteration) 

1078 

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

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

1081 

1082 def Line(self, lat1, lon1, azi1_lat2, *lon2, **name): 

1083 '''Return a geodesic line from this C{Intersector}'s geodesic, specified by 

1084 two (goedetic) points or a (goedetic) point and an (initial) azimuth. 

1085 

1086 @arg lat1: Latitude of the first point (C{degrees}). 

1087 @arg lon1: Longitude of the first point (C{degrees}). 

1088 @arg azi1_lat2: Azimuth at the first point (compass C{degrees}) if no 

1089 B{C{lon2}} argument is given, otherwise the latitude of 

1090 the second point (C{degrees}). 

1091 @arg lon2: If given, the longitude of the second point (C{degrees}). 

1092 @kwarg name: Optional C{B{name}=NN} (C{str}). 

1093 

1094 @return: A line (from L{geodesic<Intersector.geodesic>}C{.Line} or 

1095 C{.InverseLine} method) with C{LINE_CAPS}. 

1096 ''' 

1097 args = self._ll3z4ll(lat1, lon1, azi1_lat2, *lon2) 

1098 gl = self._g.InverseLine(*args, caps=Caps.LINE_CAPS) if lon2 else \ 

1099 self._g.Line( *args, caps=Caps.LINE_CAPS) 

1100 if name: 

1101 gl.name= name 

1102 return gl 

1103 

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

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

1106 

1107 def Middle(self, glA, glB, raiser=True, **_C): 

1108 '''Get the mid-points on two geodesic line segments. 

1109 

1110 @arg glA: A geodesic line (L{Line<Intersector.Line>}, 4-C{args}). 

1111 @arg glB: An other geodesic line (L{Line<Intersector.Line>}, 4-C{args}). 

1112 @kwarg raiser: If C{True}, check that B{C{glA}} and B{C{glB}} are a 

1113 4-C{args} L{Line<Intersector.Line>} or C{InverseLine} 

1114 (C{bool}). 

1115 @kwarg _C: If C{True}, include the lat-/longitudes C{latA}, C{lonA}, 

1116 C{latB}, C{lonB} of the mid-points and half-lengths C{sA} 

1117 and C{sB} in C{meter} of the respective line segments. 

1118 

1119 @return: The mid-point and half-length of each segment (L{XDict}), 

1120 B{C{_C}} above. 

1121 

1122 @raise GeodesicError: Geodesic line B{C{glA}} or B{C{glB}} invalid, 

1123 incompatible, ill-configured or not a 4-C{args 

1124 Line} or other C{InverseLine}. 

1125 ''' 

1126 M, _, _ = self._middle3(glA, glB, raiser) 

1127 return self._C(M, glA, glB, **_C) if _C else M 

1128 

1129 def _middle3(self, glA, glB, raiser): # in .All, .Segment 

1130 # return segment length C{sA} and C{sB} and the 

1131 # center C{X0} of rectangle [sA, sB] 

1132 self._xLines(glA, glB, s13=raiser) # need .Arc, .Distance 

1133 sA = glA.Distance() 

1134 sB = glB.Distance() 

1135 X = XDict_(sA * _0_5, sB * _0_5) 

1136 # _ = X._outSide(sA, sB) 

1137 return self._Delto(X), sA, sB 

1138 

1139 def Middle5(self, glA, glB, raiser=True): 

1140 '''Get the mid-points of two geodesic line segments and distances. 

1141 

1142 @return: A L{Middle5Tuple}C{(A, B, sMM, aMM, c)}. 

1143 

1144 @see: Method L{Middle} for further details. 

1145 ''' 

1146 M, _, _ = self._middle3(glA, glB, raiser) 

1147 M = self._C(M, glA, glB, _C=True, _MM=True) 

1148 A, B, s, a, c = self._In5T(glA, glB, M, M, _2X=_M_) 

1149 return Middle5Tuple(self._illz2G(A, glA), 

1150 self._illz2G(B, glB), s, a, c) 

1151 

1152 def _m12_M12_M21(self, gl, s): 

1153 P = gl.Position(s, outmask=Caps._REDUCEDLENGTH_GEODESICSCALE) 

1154 return P.m12, P.M12, P.M21 

1155 

1156 def Next(self, glA, glB, eps1=None, **_C): # PYCHOK no cover 

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

1158 

1159 @arg glA: A geodesic line (L{Line<Intersector.Line>}). 

1160 @arg glB: An other geodesic line (L{Line<Intersector.Line>}). 

1161 @kwarg eps1: Optional margin for the L{euclid<pygeodesy.euclid>}ean 

1162 distance (C{degrees}) between the C{(lat1, lon1)} points 

1163 of both lines or C{None} for unchecked. 

1164 @kwarg _C: If C{True}, include the lat-/longitudes C{latA}, C{lonA}, 

1165 C{latB}, C{lonB} of and distances C{sAB} and C{aSB} 

1166 between the intersections. 

1167 

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

1169 

1170 @raise GeodesicError: Geodesic line B{C{glA}} or B{C{glB}} invalid, 

1171 incompatible, ill-configured or C{(lat1, lon1)} 

1172 not B{C{eps1}}-equal. 

1173 

1174 @raise IntersectionError: No convergence. 

1175 

1176 @note: Offset C{X0} is implicit, zeros. 

1177 ''' 

1178 self._xLines(glA, glB) 

1179 if eps1 or _C: # eps 

1180 _C = self._xNext(glA, glB, eps1, **_C) 

1181 

1182 X0, self._conjDist3s = _X000, glA # reset Property 

1183 Q, d, S_, i = _XINF, INF, list(X0._nD2(self._D2)), 0 

1184 while S_: 

1185 X, i = self._Basic2(glA, glB, S_.pop(0), i) 

1186 X = X0._fixCoincident(X) 

1187 t = X.L1(X0) # == X.L1() 

1188 c, z = X.c, (t <= self.Delta) # X == X0 

1189 if z: 

1190 if not c: 

1191 continue 

1192 Tt_ = self._conjDist3Tt_(c, X0) 

1193 else: 

1194 Tt_ = (X, t), 

1195 

1196 for T, t in Tt_: 

1197 if t < d or Q is _XINF: 

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

1199 i += 1 

1200 

1201 for s in ((_1_1t if z else _1_0_1t) 

1202 if c else _0t): 

1203 T = X 

1204 if s and c: 

1205 s *= self._D2 

1206 T = X + (s, s * c) # NOT += 

1207 T._skip(S_, self._T2D2Delta) 

1208 

1209 return None if Q is _XINF else self._C(Q, glA, glB, **_C).set_(sX0=d, iteration=q) 

1210 

1211 def Next5(self, glA, glB, **eps1): # PYCHOK no cover 

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

1213 

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

1215 if none found. 

1216 

1217 @see: Method L{Next} for further details. 

1218 ''' 

1219 X = self.Next(glA, glB, **eps1) 

1220 return X if X is None else self._In5T(glA, glB, X, X) 

1221 

1222 def _obliqDist4(self): 

1223 zx = 45.0 

1224 if self.f: 

1225 _abs, _cjD5 = fabs, self._conjDist5 

1226 

1227 _, ds0, z0, _, _ = _cjD5(zx + _1_0) 

1228 s1, ds1, z1, sAx, sBx = _cjD5(zx - _1_0) 

1229 sx, dsx, zx = s1, _abs(ds1), z1 

1230 # find ds(azi) = 0 by secant method 

1231 for _ in range(16): 

1232 if ds1 == ds0: 

1233 break 

1234 z = (z0 * ds1 - z1 * ds0) / (ds1 - ds0) 

1235 _, ds0, z0 = s1, ds1, z1 

1236 s1, ds1, z1, a, b = _cjD5(z) 

1237 if _abs(ds1) < dsx: 

1238 sx, dsx, zx, sAx, sBx = s1, _abs(ds1), z, a, b 

1239 if not dsx: 

1240 break 

1241 else: 

1242 sx, sAx, sBx = self._cHalf, _0_5, -_1_5 

1243 return sx, zx, sAx, sBx 

1244 

1245 def _polarB3(self, lats=False): # PYCHOK no cover 

1246 latx = 64.0 

1247 lat = _90_0 - latx 

1248 if self.f: 

1249 _d, _pD2 = fdot, self._polarDist2 

1250 

1251 s0, lat0 = _pD2(latx - _1_0) 

1252 s1, lat1 = _pD2(latx + _1_0) 

1253 s2, lat2 = \ 

1254 sx, latx = _pD2(latx) 

1255 prolate = self.f < 0 

1256 # solve for ds(lat) / dlat = 0 with a quadratic fit 

1257 for _ in range(_TRIPS): 

1258 t = (lat1 - lat0), (lat0 - lat2), (lat2 - lat1) 

1259 d = _d(t, s2, s1, s0) * _2_0 

1260 if not d: # or isnan(d) 

1261 break 

1262 lat = _d(t, (lat1 + lat0) * s2, 

1263 (lat0 + lat2) * s1, 

1264 (lat2 + lat1) * s0) / d 

1265 s0, lat0 = s1, lat1 

1266 s1, lat1 = s2, lat2 

1267 s2, lat2 = _pD2(lat) 

1268 if (s2 < sx) if prolate else (s2 > sx): 

1269 sx, latx = s2, lat2 

1270 if lats: 

1271 _, lat = _pD2(latx, lat2=True) 

1272 sx += sx 

1273 else: 

1274 sx = self._cHalf 

1275 return sx, latx, lat 

1276 

1277 def _polarDist2(self, lat1, lat2=False): 

1278 gl = self._Line(lat1=lat1) 

1279 s = self._conjDist(gl, self._faPI_4, semi=True) 

1280 if lat2: 

1281 lat1 = gl.Position(s, outmask=Caps.LATITUDE).lat2 

1282 return s, lat1 

1283 

1284 def _Position(self, gl, s): 

1285 return gl.Position(s, outmask=Caps._STD_LINE) 

1286 

1287 def Segment(self, glA, glB, proven=None, raiser=True, **_C): 

1288 '''Find the intersection between two geodesic line segments. 

1289 

1290 @kwarg proven: Conjecture is that whenever two geodesic line 

1291 segments intersect, the intersection is the 

1292 one closest to the mid-points of segments. 

1293 If so, use C{B{proven}=True}, otherwise find 

1294 intersections on the segments and specify 

1295 C{B{proven}=None} to return the first or 

1296 C{B{proven}=False} the closest (C{bool}). 

1297 @kwarg raiser: If C{True}, check that B{C{glA}} and B{C{glB}} 

1298 are a 4-C{args} L{Line<Intersector.Line>} or 

1299 C{InverseLine} (C{bool}). 

1300 @kwarg _C: If C{True}, include the lat-/longitudes C{latA}, 

1301 C{lonA}, C{latB}, C{lonB} of and distances C{sAB} 

1302 and C{aSB} between the intersections. 

1303 

1304 @return: The intersection of the segments (L{XDict}) with 

1305 indicators C{kA}, C{kB} and C{k} set or if no 

1306 intersection is found, C{None}. 

1307 

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

1309 invalid, incompatible, ill-configured or 

1310 not an C{InverseLine} or 4-C{args Line}. 

1311 

1312 @raise IntersectionError: No convergence. 

1313 

1314 @see: Method L{Middle<Intersector.Middle>} for further details. 

1315 ''' 

1316 X0, sA, sB = self._middle3(glA, glB, raiser) 

1317 Q = self.Closest(glA, glB, X0) # to X0 

1318 if Q is not None: 

1319 if Q.c: # anti-/parallel 

1320 Q._fixSegment(sA, sB) 

1321 # are rectangle [sA, sB] corners further from X0 than Q? 

1322 d0 = X0.L1(Q) 

1323 if Q._outSide(sA, sB) and d0 <= X0.L1() and not proven: 

1324 i = Q.iteration 

1325 for T in Q._corners(sA, sB, self._T2): 

1326 X, i = self._Basic2(glA, glB, T, i) 

1327 X = T._fixCoincident(X) 

1328 if not X._outSide(sA, sB): 

1329 d = X0.L1(X) 

1330 if d < d0 or proven is None: 

1331 Q, d0 = X, d 

1332 if proven is None: 

1333 break 

1334 Q.set_(iteration=i) 

1335 

1336 Q = self._C(Q, glA, glB, **_C).set_(sX0=d0) 

1337 return Q 

1338 

1339 def Segment5(self, glA, glB, **proven_raiser): 

1340 '''Find the intersection between two geodesic line segments. 

1341 

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

1343 or C{None} if none found. 

1344 

1345 @see: Method L{Segment} for further details. 

1346 ''' 

1347 X = self.Segment(glA, glB, **proven_raiser) 

1348 return X if X is None else self._In5T(glA, glB, X, X, k2=True) 

1349 

1350 def _Spherical(self, glA, glB, S): 

1351 '''(INTERNAL) Get solution based from a spherical triangle. 

1352 ''' 

1353 # threshold for coincident geodesics/intersections ~4.3 nm WGS84. 

1354 A = self._Position(glA, S.sA) 

1355 B = self._Position(glB, S.sB) 

1356 D = self._Inverse(A, B) 

1357 

1358 a, da = _diff182(A.azi2, D.azi1) # interior angle at A 

1359 b, db = _diff182(B.azi2, D.azi2) # exterior angle at B 

1360 c, dc = _diff182(a, b) 

1361 if fsum1_(dc, db, -da, c) < 0: # inverted triangle 

1362 a, da = -a, -da 

1363 b, db = -b, -db 

1364 sa, ca = _sincos2de(a, da) 

1365 sb, cb = _sincos2de(b, db) 

1366 

1367 e, z, _abs = _EPS3, D.s12, fabs 

1368 if _abs(z) <= self._EPS3R: # XXX z <= ... 

1369 sA = sB = 0 # at intersection 

1370 c = 1 if _abs(sa - sb) <= e and _abs(ca - cb) <= e else ( 

1371 -1 if _abs(sa + sb) <= e and _abs(ca + cb) <= e else 0) 

1372 elif _abs(sa) <= e and _abs(sb) <= e: # coincident 

1373 sA = ca * z * _0_5 # choose mid-point 

1374 sB = -cb * z * _0_5 

1375 c = 1 if (ca * cb) > 0 else -1 

1376 # alt1: sA = ca * z; sB = 0 

1377 # alt2: sB = -cb * z; sA = 0 

1378 else: # general case 

1379 sz, cz = sincos2(z / self.R) 

1380 # [SKIP: Divide args by |sz| to avoid possible underflow 

1381 # in {sa, sb} * sz; this is probably not necessary]. 

1382 # Definitely need to treat sz < 0 (z > PI*R) correctly in 

1383 # order to avoid some convergence failures in _Basic2. 

1384 sA = atan2(sb * sz, sb * ca * cz - sa * cb) * self.R 

1385 sB = atan2(sa * sz, -sa * cb * cz + sb * ca) * self.R 

1386 c = 0 

1387 return XDict_(sA, sB, c) # no ._Delto 

1388 

1389 @Property_RO 

1390 def _T2D1Delta(self): 

1391 return self._T2d3Delta(self._D1) 

1392 

1393 @Property_RO 

1394 def _T2D2Delta(self): 

1395 return self._T2d3Delta(self._D2) 

1396 

1397 def _T2d3Delta(self, d3): 

1398 return self._T2 - d3 - self.Delta 

1399 

1400 @Property_RO 

1401 def _Tol(self): # convergence tolerance 

1402 return self._cHalf * pow(EPS, 0.75) # _0_75 

1403 

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

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

1406 

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

1408 for further details. 

1409 

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

1411 ''' 

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

1413 

1414 def _xLines(self, glA, glB, s13=False): 

1415 # check two geodesic lines vs this geodesic 

1416 C, gls = Caps.LINE_CAPS, dict(glA=glA, glB=glB) 

1417 _xinstanceof(*self._GeodesicLines, **gls) 

1418 for n, gl in gls.items(): 

1419 try: 

1420 _xgeodesics(gl.geodesic, self.geodesic) 

1421 if s13 and not isfinite(gl.s13): # or not gl.caps & Caps.DISTANCE_IN 

1422 t = gl.geodesic.InverseLine.__name__ 

1423 raise TypeError(_not_(_an(t))) 

1424 c = gl.caps & C 

1425 if c != C: # not gl.caps_(C) 

1426 c, C, x = map1(bin, c, C, _xor(c, C)) 

1427 x = _SPACE_(_xor.__name__, repr(x))[1:] 

1428 raise GeodesicError(caps=c, LINE_CAPS=C, txt=x) 

1429 except Exception as x: 

1430 raise GeodesicError(n, gl, cause=x) 

1431 

1432 

1433class Intersect7Tuple(_NamedTuple): 

1434 '''7-Tuple C{(A, B, sAB, aAB, c, kA, kB)} with C{A} and C{B} each 

1435 a C{LatLon} or L{LatLon4Tuple}C{(lat, lon, height, datum)} of 

1436 the intersection on each geodesic line, the distance C{sAB} in 

1437 in C{meter} and angular distance C{aAB} in C{degrees} between 

1438 C{A} and C{B}, coincidence indicator C{c} and segment indicators 

1439 C{kA} and C{kB} all C{int}, see L{XDict} and method U{intersect7 

1440 <_IntersectBase.intersect7>}. 

1441 ''' 

1442 _Names_ = (_A_, _B_, _sAB_, _aAB_, _c_, 'kA', 'kB') 

1443 _Units_ = (_Pass, _Pass, Meter, Degrees, Int, Int, Int) 

1444 

1445 

1446class Intersectool5Tuple(_NamedTuple): 

1447 '''5-Tuple C{(A, B, sAB, aAB, c)} with C{A} and C{B} the C{Position} 

1448 of the intersection on each geodesic line, the distance C{sAB} 

1449 between C{A} and C{B} in C{meter}, the angular distance C{aAB} in 

1450 C{degrees} and coincidence indicator C{c} (C{int}), see L{XDict}. 

1451 

1452 @note: C{A} and C{B} are each a C{GDict} with C{lat1}, C{lon1} and 

1453 C{azi1} or C{lat2}, C{lon2} from the geodesic line C{glA} 

1454 respectively C{glB} and the intersection location in C{latX}, 

1455 C{lonX}, distance C{s1X} in C{meter} and angular distance 

1456 C{a1M} in C{degrees} and the segment indicator C{kX}. See 

1457 L{XDict} for more details. 

1458 ''' 

1459 _Names_ = Intersect7Tuple._Names_[:5] 

1460 _Units_ = Intersect7Tuple._Units_[:5] 

1461 

1462 

1463class Intersector5Tuple(Intersectool5Tuple): 

1464 '''5-Tuple C{(A, B, sAB, aAB, c)} with C{A} and C{B} the C{Position} 

1465 of the intersection on each geodesic line, the distance C{sAB} 

1466 between C{A} and C{B} in C{meter}, angular distance C{aAB} in 

1467 C{degrees} and coincidence indicator C{c} (C{int}), see L{XDict}. 

1468 

1469 @note: C{A} and C{B} are each a C{GeodesicLine...Position} for 

1470 C{outmask=Caps.STANDARD} with the intersection location in 

1471 C{latX}, C{lonX}, azimuth in C{aziX}, distance C{s1X} in 

1472 C{meter} and angular distance C{a1X} in C{degrees} and the 

1473 segment indicator C{kX}. See L{XDict} for more details. 

1474 ''' 

1475 _Names_ = Intersectool5Tuple._Names_ 

1476 

1477 

1478class Middle5Tuple(Intersectool5Tuple): 

1479 '''5-Tuple C{(A, B, sMM, aMM, c)} with C{A} and C{B} the I{line segments} 

1480 including the mid-point location in C{latM}, C{lonM}, distance C{s1M} 

1481 in C{meter} and angular distance C{a1M} in C{degrees}, the distance 

1482 between both mid-points C{sMM} in C{meter} and angular distance C{aMM} 

1483 in C{degrees} and coincidence indicator C{c} (C{int}). See L{XDict} 

1484 for more details. 

1485 ''' 

1486 _Names_ = (_A_, _B_, 'sMM', 'aMM', _c_) 

1487 

1488 

1489class _List(list): 

1490 

1491 _Delta = 0 # equality margin 

1492 

1493 def __init__(self, Delta): 

1494 self._Delta = Delta 

1495# list.__init__(self) 

1496 

1497 def __contains__(self, other): 

1498 # handle C{if X in this: ...} 

1499 a, b = other.sA, other.sB 

1500 D, _D1 = self._Delta, _L1 

1501 for X in self: 

1502 if _D1(X.sA - a, X.sB - b) <= D: 

1503 return True 

1504 return False 

1505 

1506 def addend(self, X, *d0_i): 

1507 # append an item, updated 

1508 if d0_i: 

1509 d0, i = d0_i 

1510 X.set_(sX0=d0, iteration=i) 

1511 self.append(X) 

1512 return X.sX0 

1513 

1514 def sorter(self, sMaX0, dot_C, glA, glB, **_C): 

1515 # trim and sort the X items 

1516 

1517 def _key(X): 

1518 return X.sX0 # rank of X 

1519 

1520 t = (X for X in self if X.sX0 <= sMaX0) 

1521 for X in sorted(t, key=_key): 

1522 yield dot_C(X, glA, glB, **_C) if _C else X 

1523 

1524 

1525def _L1(a, b): 

1526 '''(INTERNAL) Return the I{L1} distance. 

1527 ''' 

1528 return fabs(a) + fabs(b) 

1529 

1530 

1531__all__ += _ALL_DOCS(_IntersectBase) 

1532 

1533if __name__ == '__main__': # MCCABE 14 

1534 

1535 from pygeodesy import printf 

1536 __help_ = '--help' 

1537 

1538 def _main(args): 

1539 

1540 from pygeodesy import GeodesicExact 

1541 from pygeodesy.internals import _plural, _usage 

1542 from pygeodesy.interns import _COLONSPACE_, _DOT_, _EQUAL_, \ 

1543 _i_, _m_, _n_, _version_, _X_ 

1544 import re 

1545 

1546 class XY0(Float): 

1547 pass 

1548 

1549 def _opts(_h): # for _usage() 

1550 ll4 = ' latA1 lonA1' 

1551 ll4 += ll4.replace('1', '2') 

1552 ll4 += ll4.replace(_A_, _B_) 

1553 llz = _SPACE_(NN, _latA_, _lonA_, 'aziA') 

1554 llz2 = llz + llz.replace(_A_, _B_) 

1555 return dict(opts='-Verbose|V--version|v--help|h--Tool|T--Check|C-R meter-', 

1556 alts=((_c_ + llz2), 

1557 (_i_ + ll4), 

1558 (_m_ + ll4), 

1559 (_n_ + llz + ' aziB'), 

1560 ('o' + llz2 + ' x0 y0')), 

1561 help=_h if isinstance(_h, str) else NN) 

1562 

1563 def _starts(Opt, arg): 

1564 return arg == Opt[1:3] or (len(arg) > 2 and Opt.startswith(arg)) 

1565 

1566 _isopt = re.compile('^[-]+[a-z]*$', flags=re.IGNORECASE).match 

1567 

1568 I = Intersector(GeodesicExact()) # PYCHOK I 

1569 M = m = _R = None 

1570 _T = _V = _h = _C = False 

1571 

1572 while args and _isopt(args[0]): 

1573 arg = args.pop(0) 

1574 if arg == _c__: 

1575 M, m = I.Closest, 6 # latA lonA aziA latB lonB aziB 

1576 elif _starts('--Check', arg): 

1577 _C = True 

1578 elif _starts(__help_, arg): 

1579 _h = args[0] if args and _isopt(args[0]) else True 

1580 elif arg == _i__: 

1581 M, m = I.Segment, 8 # latA1 lonA1 latA2 lonA2 latB1 lonB1 latB2 lonB2 

1582 elif arg == '-m': 

1583 M, m = I.Middle, 8 # latA1 lonA1 latA2 lonA2 latB1 lonB1 latB2 lonB2 

1584 _R = None # zap -R 

1585 elif arg == _n__: 

1586 M, m = I.Next, 4 # latA lonA aziA aziB 

1587 elif arg == _o__: 

1588 M, m = I.Closest, 8 # latA lonA aziA latB lonB aziB x0 y0 

1589 elif arg == _R__ and args: 

1590 _R = args.pop(0) 

1591 elif _starts('--Tool', arg): 

1592 I = Intersectool() # PYCHOK I 

1593 if _V: 

1594 I.verbose = True 

1595 if I.IntersectTool in (_PYGEODESY_INTERSECTTOOL_, None): # not set 

1596 I.IntersectTool = '/opt/local/bin/IntersectTool' # '/opt/local/Cellar/geographiclib/2.3/bin/IntersectTool' # HomeBrew 

1597 elif _V: 

1598 _ = I.version 

1599 M, _T = None, True 

1600 elif _starts('--Verbose', arg): 

1601 _V = True 

1602 if _T: 

1603 I.verbose = True 

1604 elif _starts('--version', arg): 

1605 printf(_COLONSPACE_(*((_version_, I.version) if _T else 

1606 (__version__, repr(I))))) 

1607 else: 

1608 raise ValueError('invalid option %r' % (arg,)) 

1609 

1610 if _h or M is None: 

1611 printf(_usage(__file__, **_opts(_h)), nl=1) 

1612 else: 

1613 n = len(args) 

1614 if n < m: 

1615 n = _plural('only %s arg' % (n,), n) if n else 'no args' 

1616 raise ValueError('%s, need %s' % (n, m)) 

1617 args[:] = args[:m] 

1618 

1619 kwds = dict(_C=True) if _C else {} 

1620 if M == I.Next: # -n 

1621 # get latA lonA aziA latA lonA aziB 

1622 args[3:] = args[:2] + args[3:4] 

1623 elif M == I.Closest and m > 6: # -o 

1624 y0 = Meter(y0=args.pop()) 

1625 x0 = Meter(x0=args.pop()) 

1626 kwds.update(X0=XDict_(x0, y0)) 

1627 if _R: 

1628 m = Meter_(_R, name=_R__, low=0) 

1629 kwds.update(sMaX0=m) 

1630 M = I.All 

1631 

1632 n = len(args) // 2 

1633 glA = I.Line(*args[:n]) 

1634 glB = I.Line(*args[n:]) 

1635 

1636 m = _DOT_(I.__class__.__name__, M.__name__) 

1637 if _V: 

1638 X = _SPACE_(_X_, _EQUAL_, m) 

1639 printf(unstr(X, glA, glB, **kwds)) 

1640 

1641 X = M(glA, glB, **kwds) 

1642 if X is None or isinstance(X, (XDict, tuple)): 

1643 printf(_COLONSPACE_(m, repr(X))) 

1644 else: 

1645 for i, X in enumerate(X): 

1646 printf(_COLONSPACE_(Fmt.INDEX(m, i), repr(X))) 

1647 

1648 from sys import argv, stderr 

1649 try: 

1650 if len(argv) == 2 and argv[1] == __help_: 

1651 from pygeodesy.internals import _usage_argv 

1652 

1653 s = _SPACE_(*_usage_argv(__file__)) 

1654 for t in ('-h', '-h -n', 

1655 '-c 0 0 45 40 10 135', 

1656 '-C -c 0 0 45 40 10 135', 

1657 '-T -R 2.6e7 -c 0 0 45 40 10 135', 

1658 '-c 50 -4 -147.7 0 0 90', 

1659 '-C -c 50 -4 -147.7 0 0 90', 

1660 '# % echo 0 0 10 10 50 -4 50S 4W | IntersectTool -i -p 0 -C', 

1661 '# -631414 5988887 0 -3', 

1662 '# -4.05187 -4.00000 -4.05187 -4.00000 0', 

1663 '-m 0 0 10 10 50 -4 50S 4W', 

1664 '-C -m 0 0 10 10 50 -4 50S 4W', 

1665 '-i 0 0 10 10 50 -4 50S 4W', 

1666 '-T -i 0 0 10 10 50 -4 50S 4W', 

1667 '-C -i 0 0 10 10 50 -4 50S 4W', 

1668 '-T -C -i 0 0 10 10 50 -4 50S 4W', 

1669 '-V -T -i 0 0 10 10 50 -4 -50 -4', 

1670 '-C -R 4e7 -c 50 -4 -147.7 0 0 90', 

1671 '-T -C -R 4e7 -c 50 -4 -147.7 0 0 90', 

1672 '-R 4e7 -i 0 0 10 10 50 -4 -50 -4', 

1673 '-T -R 4e7 -i 0 0 10 10 50 -4 -50 -4'): 

1674 if t.startswith(_HASH_): 

1675 printf(t, nl=int(t[2] == '%')) 

1676 else: 

1677 printf(_SPACE_(_HASH_, s, t), nl=1) 

1678 argv[1:] = t = t.split() 

1679 _main(t) 

1680 else: 

1681 _main(argv[1:]) 

1682 

1683 except Exception as x: 

1684 x = _SPACE_(x, NN, _HASH_, *argv) 

1685 printf(x, file=stderr, nl=1) 

1686 if '-V' in x or _MODS.errors.exception_chaining(): 

1687 raise 

1688 exit(1) 

1689 

1690# % env PYGEODESY_INTERSECTTOOL=... python3 -m pygeodesy.geodesici --help 

1691 

1692# % python3 -m pygeodesy.geodesici -h 

1693# 

1694# usage: python3 -m ....pygeodesy.geodesici [--Verbose | -V] [--version | -v] [--help | -h] [--Tool | -T] [--Check | -C] [-R meter] 

1695# [-c latA lonA aziA latB lonB aziB | 

1696# -i latA1 lonA1 latA2 lonA2 latB1 lonB1 latB2 lonB2 | 

1697# -m latA1 lonA1 latA2 lonA2 latB1 lonB1 latB2 lonB2 | 

1698# -n latA lonA aziA aziB | 

1699# -o latA lonA aziA latB lonB aziB x0 y0] 

1700 

1701# % python3 -m ....pygeodesy.geodesici -h -n 

1702# 

1703# usage: python3 -m ....pygeodesy.geodesici -n latA lonA aziA aziB 

1704 

1705# % python3 -m ....pygeodesy.geodesici -c 0 0 45 40 10 135 

1706# Intersector.Closest: XDict(c=0, sA=3862290.547855, sB=2339969.547699, sX0=6202260.095554) 

1707 

1708# % python3 -m ....pygeodesy.geodesici -C -c 0 0 45 40 10 135 

1709# Intersector.Closest: XDict(aAB=0.0, c=0, latA=23.875306, latB=23.875306, lonA=26.094096, lonB=26.094096, sA=3862290.547855, sAB=0.0, sB=2339969.547699, sX0=6202260.095554) 

1710 

1711# % env PYGEODESY_INTERSECTTOOL=...python3 -m ....pygeodesy.geodesici -T -R 2.6e7 -c 0 0 45 40 10 135 

1712# Intersectool.All[0]: XDict(c=0, sA=3862290.547855, sB=2339969.547699, sX0=6202260.095554) 

1713 

1714# % python3 -m ....pygeodesy.geodesici -c 50 -4 -147.7 0 0 90 

1715# Intersector.Closest: XDict(c=0, sA=6058048.653081, sB=-3311252.995823, sX0=9369301.648903) 

1716 

1717# % python3 -m ....pygeodesy.geodesici -C -c 50 -4 -147.7 0 0 90 

1718# Intersector.Closest: XDict(aAB=0.0, c=0, latA=0.0, latB=-0.0, lonA=-29.745492, lonB=-29.745492, sA=6058048.653081, sAB=0.0, sB=-3311252.995823, sX0=9369301.648903) 

1719 

1720# % echo 0 0 10 10 50 -4 50S 4W | IntersectTool -i -p 0 -C 

1721# -631414 5988887 0 -3 

1722# -4.05187 -4.00000 -4.05187 -4.00000 0 

1723 

1724# % python3 -m ....pygeodesy.geodesici -m 0 0 10 10 50 -4 50S 4W 

1725# Intersector.Middle: XDict(c=0, sA=782554.549609, sB=5536835.161499, sX0=0.0) 

1726 

1727# % python3 -m ....pygeodesy.geodesici -C -m 0 0 10 10 50 -4 50S 4W 

1728# Intersector.Middle: XDict(aAB=10.262308, c=0, latA=5.019509, latB=0.036282, lonA=4.961883, lonB=-4.0, sA=782554.549609, sAB=1138574.546746, sB=5536835.161499, sX0=0.0) 

1729 

1730# % python3 -m ....pygeodesy.geodesici -i 0 0 10 10 50 -4 50S 4W 

1731# Intersector.Segment: XDict(c=0, k=-3, kA=-1, kB=0, sA=-631414.26877, sB=5988887.278435, sX0=1866020.935315) 

1732 

1733# % env PYGEODESY_INTERSECTTOOL=... python3 -m ....pygeodesy.geodesici -T -i 0 0 10 10 50 -4 50S 4W 

1734# Intersectool.Segment: XDict(c=0, k=-3, kA=-1, kB=0, sA=-631414.26877, sB=5988887.278435) 

1735 

1736# % python3 -m ....pygeodesy.geodesici -C -i 0 0 10 10 50 -4 50S 4W 

1737# Intersector.Segment: XDict(aAB=0.0, c=0, k=-3, kA=-1, kB=0, latA=-4.051871, latB=-4.051871, lonA=-4.0, lonB=-4.0, sA=-631414.26877, sAB=0.0, sB=5988887.278435, sX0=1866020.935315) 

1738 

1739# % env PYGEODESY_INTERSECTTOOL=... python3 -m ....pygeodesy.geodesici -T -C -i 0 0 10 10 50 -4 50S 4W 

1740# Intersectool.Segment: XDict(c=0, k=-3, kA=-1, kB=0, latA=-4.051871, latB=-4.051871, lonA=-4.0, lonB=-4.0, sA=-631414.26877, sAB=0.0, sB=5988887.278435) 

1741 

1742# % env PYGEODESY_INTERSECTTOOL=... python3 -m ....pygeodesy.geodesici -V -T -i 0 0 10 10 50 -4 -50 -4 

1743# Intersectool@1: /opt/local/bin/IntersectTool --version (invoke) 

1744# Intersectool@1: '/opt/local/bin/IntersectTool: GeographicLib version 2.3' (0) 

1745# Intersectool@1: /opt/local/bin/IntersectTool: GeographicLib version 2.3 (0) 

1746# X = Intersectool.Segment(GDict(lat1=0.0, lat2=10.0, lon1=0.0, lon2=10.0), GDict(lat1=50.0, lat2=-50.0, lon1=-4.0, lon2=-4.0)) 

1747# Intersectool@2: /opt/local/bin/IntersectTool -E -p 10 -i \ 0.0 0.0 10.0 10.0 50.0 -4.0 -50.0 -4.0 (Segment) 

1748# Intersectool@2: '-631414.2687702414 5988887.2784352796 0 -3' (0) 

1749# Intersectool@2: sA=-631414.2687702414, sB=5988887.2784352796, c=0, k=-3 (0) 

1750# Intersectool.Segment: XDict(c=0, k=-3, kA=-1, kB=0, sA=-631414.26877, sB=5988887.278435) 

1751 

1752# % python3 -m ....pygeodesy.geodesici -C -R 4e7 -c 50 -4 -147.7 0 0 90 

1753# Intersector.All[0]: XDict(aAB=0.0, c=0, latA=0.0, latB=-0.0, lonA=-29.745492, lonB=-29.745492, sA=6058048.653081, sAB=0.0, sB=-3311252.995823, sX0=9369301.648903) 

1754# Intersector.All[1]: XDict(aAB=0.0, c=0, latA=0.0, latB=0.0, lonA=150.046964, lonB=150.046964, sA=-13941907.021445, sAB=0.0, sB=16703151.659744, sX0=30645058.681189) 

1755# Intersector.All[2]: XDict(aAB=0.0, c=0, latA=-0.0, latB=-0.0, lonA=-30.16058, lonB=-30.16058, sA=-33941862.69597, sAB=0.0, sB=-3357460.370268, sX0=37299323.066238) 

1756# Intersector.All[3]: XDict(aAB=0.0, c=0, latA=-0.0, latB=0.0, lonA=150.046964, lonB=150.046964, sA=-13941907.021445, sAB=0.0, sB=-23371865.025835, sX0=37313772.047279) 

1757 

1758# % env PYGEODESY_INTERSECTTOOL=... python3 -m ....pygeodesy.geodesici -T -C -R 4e7 -c 50 -4 -147.7 0 0 90 

1759# Intersectool.All[0]: XDict(c=0, latA=-0.0, latB=-0.0, lonA=-29.745492, lonB=-29.745492, sA=6058048.653081, sAB=0.0, sB=-3311252.995823, sX0=9369301.648903) 

1760# Intersectool.All[1]: XDict(c=0, latA=0.0, latB=0.0, lonA=150.046964, lonB=150.046964, sA=-13941907.021445, sAB=0.0, sB=16703151.659744, sX0=30645058.681189) 

1761# Intersectool.All[2]: XDict(c=0, latA=-0.0, latB=-0.0, lonA=-30.16058, lonB=-30.16058, sA=-33941862.69597, sAB=0.0, sB=-3357460.370268, sX0=37299323.066238) 

1762# Intersectool.All[3]: XDict(c=0, latA=-0.0, latB=0.0, lonA=150.046964, lonB=150.046964, sA=-13941907.021445, sAB=0.0, sB=-23371865.025835, sX0=37313772.047279) 

1763 

1764# % python3 -m ....pygeodesy.geodesici -R 4e7 -i 0 0 10 10 50 -4 -50 -4 

1765# Intersector.All[0]: XDict(c=0, sA=-631414.26877, sB=5988887.278435, sX0=1866020.935315) 

1766# Intersector.All[1]: XDict(c=0, sA=19422725.117572, sB=-14062417.105648, sX0=38239422.83511) 

1767# Intersector.All[2]: XDict(c=0, sA=19422725.117572, sB=25945445.811603, sX0=39048781.218067) 

1768# Intersector.All[3]: XDict(c=0, sA=39476927.464575, sB=5894074.699478, sX0=39051612.452944) 

1769 

1770# % env PYGEODESY_INTERSECTTOOL=... python3 -m ....pygeodesy.geodesici -T -R 4e7 -i 0 0 10 10 50 -4 -50 -4 

1771# Intersectool.All[0]: XDict(c=0, sA=-631414.26877, sB=5988887.278435, sX0=1862009.05513) 

1772# Intersectool.All[1]: XDict(c=0, sA=19422725.117572, sB=-14062417.105648, sX0=38243434.715295) 

1773# Intersectool.All[2]: XDict(c=0, sA=19422725.117572, sB=25945445.811603, sX0=39044769.337882) 

1774# Intersectool.All[3]: XDict(c=0, sA=39476927.464575, sB=5894074.699478, sX0=39047600.57276) 

1775 

1776 

1777# **) MIT License 

1778# 

1779# Copyright (C) 2024-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

1780# 

1781# Permission is hereby granted, free of charge, to any person obtaining a 

1782# copy of this software and associated documentation files (the "Software"), 

1783# to deal in the Software without restriction, including without limitation 

1784# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

1785# and/or sell copies of the Software, and to permit persons to whom the 

1786# Software is furnished to do so, subject to the following conditions: 

1787# 

1788# The above copyright notice and this permission notice shall be included 

1789# in all copies or substantial portions of the Software. 

1790# 

1791# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

1792# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

1793# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

1794# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

1795# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

1796# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

1797# OTHER DEALINGS IN THE SOFTWARE.