1

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

3

4 u'''Classes L{Intersectool} and L{Intersector}, a pure Python version of I{Karney}'s C++ class

5 U{Intersect<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Intersect.html>}

6 to find the intersections of two geodesic lines or line segments.

7

8 L{Intersectool} and L{Intersector} methods C{All}, C{Closest}, C{Next} and C{Segment} produce

9

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

11 C{Segment} return or yield L{Intersectool5Tuple} or L{Intersector5Tuple}s with the lat-, longitude

12 azimuth of each intersection as an extended or C{Position}-like L{GDict}.

13

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

15<https://GeographicLib.SourceForge.io/C++/doc/IntersectTool.1.html>} utility like class L{Intersector},

16 but intended I{for testing purposes only}.

17

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

19

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

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

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

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

24'''

25# make sure int/int division yields float quotient

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

27

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

29 _xinstanceof, _xor

30 from 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

33 from pygeodesy.ellipsoids import _EWGS84, Fmt, unstr

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

35 _xgeodesics, _xkwds_get, _xkwds_kwds

36# from pygeodesy.errors import exception_chaining # _MODS

37 from pygeodesy.fmath import euclid, fdot

38 from pygeodesy.fsums import Fsum, fsum1_, _ceil

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

40 _HASH_, _M_, _not_, _SPACE_, _too_

41 from pygeodesy.interns import _m_ # PYCHOK used!

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

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

44 _getenv, _PYGEODESY_INTERSECTTOOL_

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

46 from pygeodesy.props import deprecated_method, Property, \

47 Property_RO, property_RO

48 from pygeodesy.solveBase import _SolveCapsBase, pairs

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

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

51 from pygeodesy.units import Degrees, Float, Int, Lat, Lon, \

52 Meter, Meter_

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

54

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

56

57__all__ = _ALL_LAZY.geodesici

58__version__ = '24.07.09'

59

60_0t = 0, # int

61_1_1t = -1, +1

62_1_0_1t = -1, 0, +1

63_c__ = '-c' # PYCHOK used!

64_C__ = '-C' # PYCHOK used!

65_EPS3 = EPS * _3_0

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

67_i__ = '-i' # PYCHOK used!

68_latA_ = 'latA'

69_lonA_ = 'lonA'

70_n__ = '-n' # PYCHOK used!

71_o__ = '-o' # PYCHOK used!

72_R__ = '-R'

73_sAB_ = 'sAB'

74_sX0_ = 'sX0'

75_TRIPS = 128

76

77

78 class Azi(Degrees):

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

80 '''

81 pass

82

83

84 class XDict(ADict):

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

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

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

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

89 or 0 otherwise.

90

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

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

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

94

95 Segment indicators C{kA} and C{kB} are C{0} if the segments

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

97 the start, respectively I{after} the end of the segment, like

98 L{Intersection3Tuple<Intersection3Tuple>}. Segment indicator

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

100 '''

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

102

103 def __add__(self, other):

104 X = _copy(self)

105 X += other

106 return X

107

108 def __eq__(self, other):

109 return not self.__ne__(other)

110

111 def __iadd__(self, other):

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

113 a, b = other

114 else:

115 # _xinstanceof(XDict, other=other)

116 a = other.sA

117 b = other.sB

118 if other.c:

119 self.c = other.c

120 self.sA += a # PYCHOK sA

121 self.sB += b # PYCHOK sB

122 return self

123

124 def __le__(self, other):

125 # _xinstanceof(XDict, other=other)

126 return self == other or self < other

127

128 def __lt__(self, other):

129 # _xinstanceof(XDict, other=other)

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

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

132

133 def __ne__(self, other):

134 # _xinstanceof(XDict, other=other)

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

136

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

138 # yield all corners further than C{T2}

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

140 for x in (0, sA):

141 for y in (0, sB):

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

143 yield XDict_(x, y)

144

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

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

147 c = c0 or X.c

148 if c:

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

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

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

152 return X

153

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

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

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

157

158 def _g(): # intersection in smallest gap

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

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

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

162 s = a + b

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

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

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

166 s = sB + (a - b)

167 return (sB if t else sA) - s

168

169 ta = -a

170 tb = sA - a

171 tc = -c * b

172 td = -c * (b - sB)

173

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

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

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

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

178

179 # test opposite rectangle sides first

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

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

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

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

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

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

186 self += s, s * c

187

188 @property_RO

189 def _is00(self):

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

191

192 def L1(self, other=None):

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

194 '''

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

196 if other is not None:

197 # _xinstanceof(XDict, other=other)

198 a -= other.sA

199 b -= other.sB

200 return _L1(a, b)

201

202 def _nD1(self, D1):

203 # yield the C{Closest} starts

204 D_ = 0, D1, -D1

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

206 (0, 0, 0, 1, -1)):

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

208

209 def _nD2(self, D2):

210 # yield the C{Next} starts

211 D22 = D2 * _2_0

212 D_ = 0, D2, D22, -D22, -D2

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

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

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

216

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

218 # yield the C{All} starts

219 yield self

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

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

222 if i or j: # skip self

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

224 (i - j) * D3)

225

226 def _outSide(self, sA, sB):

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

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

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

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

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

232 return bool(kA or kB)

233

234 def _skip(self, S_, T1_Delta):

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

236 for j, S in _enumereverse(S_):

237 if S.L1(self) < T1_Delta:

238 S_.pop(j)

239

240

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

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

243 '''

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

245

246_X000 = XDict_() # PYCHOK origin

247_XINF = XDict_(INF)

248

249

250 class _IntersectBase(_NamedBase):

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

252 '''

253 # _g = None

254

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

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

257 self._g = geodesic

258 if name:

259 self.name = name

260

261 @Property_RO

262 def a(self):

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

264 '''

265 return self.ellipsoid.a

266

267 equatoradius = a # = Requatorial

268

269 @Property_RO

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

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

272

273 @Property_RO

274 def _cMax(self): # outer circumference

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

276

277 @property_RO

278 def datum(self):

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

280 '''

281 return self.geodesic.datum

282

283 @Property_RO

284 def ellipsoid(self):

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

286 '''

287 return self.geodesic.datum.ellipsoid

288

289 @Property_RO

290 def f(self):

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

292 '''

293 return self.ellipsoid.f

294

295 flattening = f

296

297 @property_RO

298 def geodesic(self):

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

300 '''

301 return self._g

302

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

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

305 B.lat2, B.lon2)

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

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

308

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

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

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

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

313 else:

314 t += Azi(azi1=azi1_lat2),

315 return t

316

317 @deprecated_method

318 def Next5s(self, glA, glB, X0=_X000, aMax=1801, sMax=0, **unused):

319 '''DEPRECATED on 2024.07.02, use method L{All5}.'''

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

321

322 @Property_RO

323 def R(self):

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

325 '''

326 return self.ellipsoid.R2

327

328 def _sMaX0_C2(self, aMaX0, **sMaX0_C):

329 _g = _xkwds_get

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

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

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

333 if a: # degrees to meter, approx.

334 s = max(s, self.R * radians(a)) # ellipsoid.degrees2m(a)

335 s = _g(sMaX0_C, _R=s)

336 if s < _EPS3:

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

338 return s, _xkwds_kwds(sMaX0_C, _C=False)

339

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

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

342 if eps1 is not None:

343 a = glA.lat1 - glB.lat1

344 b = glA.lon1 - glB.lon1

345 if euclid(a, b) > eps1:

346 raise GeodesicError(lat=a, lon=b, eps1=eps1)

347 return _xkwds_kwds(eps_C, _C=False)

348

349

350 class Intersectool(_IntersectBase, _SolveCapsBase):

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

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

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

354

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

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

357

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

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

360 '''

361 _Error = GeodesicError

362 _linelimit = 1200 # line printer width X 10

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

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

365 _Xable_name = 'IntersectTool'

366 _Xable_path = _getenv(_PYGEODESY_INTERSECTTOOL_, _PYGEODESY_INTERSECTTOOL_)

367

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

369 '''New L{IntersectTool}.

370

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

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

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

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

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

376 is in C{meter}, ignored otherwise.

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

378

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

380 large or no initial convergence.

381

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

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

384 '''

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

386 self._GeodesicExact(a_geodesic, f))

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

388

389 def All(self, glA, glB, X0=_X000, eps1=_0_0, aMaX0=0, **sMaX0_C):

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

391

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

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

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

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

396

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

398 '''

399 def _xz2(**gl):

400 try:

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

402 try:

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

404 except (AttributeError, KeyError):

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

406 except Exception as x:

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

408

409 _t, a = _xz2(glA=glA)

410 _x, b = _xz2(glB=glB)

411 if _x is not _t:

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

413

414 A = glA.lat1, glA.lon1

415 B = glB.lat1, glB.lon1

416 if _x is self._c_alt:

417 if X0 is _X000 or X0._is00:

418 if eps1 is not None and \

419 euclid(glA.lat1 - glB.lat1,

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

421 _x, B = self._n_alt, ()

422 else: # non-zero offset

423 _x = self._o_alt

424 b += X0.sA, X0.sB

425

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

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

428 _R=sMaX0, **_C):

429 yield X.set_(c=int(X.c))

430

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

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

433

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

435 '''

436 for X in self.All(glA, glB, X0=X0, _C=True, **aMaX0_sMaX0):

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

438

439 @Property_RO

440 def _C_option(self):

441 return (_C__,)

442

443 @Property_RO

444 def _cmdBasic(self):

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

446 '''

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

448 self._E_option +

449 self._p_option)

450

451 @Property_RO

452 def _c_alt(self):

453 '''(INTERNAL) Get the C{Closest} format (C{tuple}).

454 '''

455 return _c__, # latA lonA aziA latB lonB aziB

456

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

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

459

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

461

462 @return: An L{XDict}.

463 '''

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

465 glB.lat1, glB.lon1, glB.azi1

466 if X0 is _X000 or X0._is000:

467 _x = self._c_alt

468 else:

469 _x = self._o_alt

470 args += X0.sA, X0.sB

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

472

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

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

475

476 @return: An L{Intersectool5Tuple}.

477 '''

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

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

480

481 @property_RO

482 def _GeodesicExact(self):

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

484 '''

485 Intersectool._GeodesicExact = G = _MODS.geodesicx.GeodesicExact # overwrite property_RO

486 return G

487

488 @Property_RO

489 def _i_alt(self):

490 '''(INTERNAL) Get the C{Segment} format (C{tuple}).

491 '''

492 return _i__, # latA1 lonA1 latA2 lonA2 latB1 lonB1 latB2 lonB2

493

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

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

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

497 if k2:

498 A.set_(k2=X.kA)

499 B.set_(k2=X.kB)

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

501 sAB = _xkwds_get(X, sAB=s)

502 if a and s and s != sAB:

503 a *= sAB / s # adjust a

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

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

506

507 @Property

508 def IntersectTool(self):

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

510 executable (C{filename}).

511 '''

512 return self._Xable_path

513

514 @IntersectTool.setter # PYCHOK setter!

515 def IntersectTool(self, path):

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

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

518

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

520 C{IntersectTool} executable.

521 '''

522 self._setXable(path)

523

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

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

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

527

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

529 '''

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

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

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

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

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

535 if name:

536 gl.name= name

537 return gl

538

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

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

541

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

543

544 @return: An L{XDict}.

545 '''

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

547 return X

548

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

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

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

552

553 def _smi4(**gl):

554 try:

555 n, gl = gl.popitem()

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

557 except Exception as x:

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

559 s = il.s13

560 m = s * _0_5

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

562

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

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

565 X = XDict_(mA, mB) # center

566 _ = X._outSide(sA, sB)

567 if _C: # _Names_ABs

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

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

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

571 return X, A, iA, B, iB

572

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

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

575

576 @return: A L{Middle5Tuple}.

577 '''

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

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

580 return Middle5Tuple(_llz2G(A, glA, iA), _llz2G(B, glB, iB), s, a, c)

581

582 @Property_RO

583 def _n_alt(self):

584 '''(INTERNAL) Get the C{Next} format (C{tuple}).

585 '''

586 return _n__, # latA lonA aziA aziB

587

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

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

590

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

592

593 @return: An L{XDict}.

594 '''

595 if eps1 or _C:

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

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

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

599 **_C) # _R=None

600

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

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

603

604 @return: An L{Intersectool5Tuple}.

605 '''

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

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

608

609 @Property_RO

610 def _o_alt(self):

611 '''(INTERNAL) Get the C{Offset} format (C{tuple}).

612 '''

613 return _o__, # latA lonA aziA latB lonB aziB x0 y0

614

615 def _R_option(self, _R=None):

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

617 '''

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

619

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

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

622

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

624

625 @return: An L{XDict}.

626 '''

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

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

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

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

631 k = int(X.k)

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

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

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

635 return X.set_(k=k, kA=kA, kB=kB)

636 raise GeodesicError(k=k, glA=glA, glB=glB, txt=repr(X))

637

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

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

640

641 @return: An L{Intersectool5Tuple}.

642 '''

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

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

645

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

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

648

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

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

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

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

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

654

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

656 '''

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

658 status=self.status,

659 IntersectTool=self.IntersectTool)

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

661

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

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

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

665 '''

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

667 # self._i_alt: 8,

668 # self._n_alt: 4,

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

670 cmd = self._cmdBasic

671 Names = self._Names_XDict # has _c_ always

672 if _k_sX0:

673 Names += _k_sX0,

674 if _C:

675 cmd += self._C_option

676 Names += self._Names_ABs

677 if _R:

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

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

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

681

682

683 class Intersector(_IntersectBase):

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

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

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

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

688

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

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

691 '''

692 # _D1 = 0

693 # _D2 = 0

694 # _T1 = 0

695 # _T5 = 0

696

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

698 '''New L{Intersector}.

699

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

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

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

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

704

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

706 large or no initial convergence.

707

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

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

710 '''

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

712 E = self.ellipsoid

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

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

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

716 if self.f > 0:

717 t3 = self._obliqDist4()[0]

718 t4 = t1

719 else: # PYCHOK no cover

720 t1, t2, t3 = t2, t1, t5

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

722

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

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

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

726 self._T1 = t1 # min distance between intersects

727 self._T2 = t2 = t1 * _2_0

728# self._T5 = t5 # not used

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

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

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

732

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

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

735

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

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

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

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

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

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

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

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

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

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

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

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

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

749

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

751

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

753 invalid, incompatible or ill-configured.

754

755 @raise IntersectionError: No convergence.

756 '''

757 self._xLines(glA, glB)

758 if X0 is None:

759 try: # determine X0

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

761 except GeodesicError: # no .Distance

762 X0 = _X000

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

764

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

766 xMaX0 = sMaX0 + D

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

768 d3 = xMaX0 / m

769 T2d3D = self._T2d3Delta(d3)

770

771 C_ = _List(D) # closest coincident

772 X_ = _List(D) # intersections found

773 c0 = 0

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

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

776 while S_:

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

778 if Q in X_:

779 continue

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

781 _X0fx = X0._fixCoincident

782 Q = _X0fx(Q) # Q = Q'

783 if c0 and Q in C_:

784 continue

785 C_.addend(Q)

786 # elimate all existing intersections

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

788 c0 = Q.c

789 for j, X in _enumereverse(X_):

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

791 X_.pop(j)

792

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

794 args = self._m12_M12_M21(glA, s0)

795 _cjD = self._conjDist

796 for s in (-_D, _D):

797 s += s0

798 sa = 0

799 while True:

800 i += 1

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

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

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

804 break

805

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

807 continue

808 else:

809 a = len(X_)

810

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

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

813 X._skip(S_, T2d3D)

814

815 return X_.sorter(sMaX0, self._C, glA, glB, **_C) # generator

816

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

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

819

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

821 for each intersection found.

822

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

824 '''

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

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

827

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

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

830 '''

831 X = _copy(S)

832 for _ in range(_TRIPS):

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

834 X += S

835 i += 1

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

837 return self._Delto(X), i

838

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

840

841 def _C(self, X, glA, glB, _C=False, _AB=True):

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

843 if _C:

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

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

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

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

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

849 if _AB:

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

851 else:

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

853 return X

854

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

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

857

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

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

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

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

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

863 and C{aSB} between the intersections.

864

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

866

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

868 invalid, incompatible or ill-configured.

869

870 @raise IntersectionError: No convergence.

871 '''

872 self._xLines(glA, glB)

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

874 while S_:

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

876 X = X0._fixCoincident(X)

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

878 d0 = X.L1(X0)

879 if d0 < self._T1:

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

881 break

882 if d0 < d or Q is X0:

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

884 X._skip(S_, self._T2D1Delta)

885

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

887

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

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

890

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

892 or C{None} if none found.

893

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

895 '''

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

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

898

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

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

901 # if semi:

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

903 # else:

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

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

906 for _ in range(_TRIPS):

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

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

909 m23 = m13 * M12

910 M32 = M31 * M12

911 if m12: # PYCHOK no cover

912 m23 -= m12 * M13

913 if m13:

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

915 if semi:

916 M23 = M13 * M21

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

918 if m12 and m13:

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

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

921 else:

922 d = -m23 / M32

923 s, d = _S2(d)

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

925 break

926 return s

927

928 _gl3 = None

929

930 @Property

931 def _conjDist3s(self):

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

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

934

935 @_conjDist3s.setter # PYCHOK setter!

936 def _conjDist3(self, gl):

937 # _XLines(gl, gl)

938 self._gl3 = gl

939

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

941 for s in self._conjDist3s:

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

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

944

945 def _conjDist5(self, azi):

946 gl = self._Line(azi1=azi)

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

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

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

950

951 @Property_RO

952 def Delta(self):

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

954 '''

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

956

957 def _Delto(self, X):

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

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

960 return X

961

962 @Property_RO

963 def _EPS3R(self):

964 return _EPS3 * self.R

965

966 @Property_RO

967 def _faPI_4(self):

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

969

970 @Property_RO

971 def _GeodesicLines(self):

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

973 '''

974 return type(self._Line()),

975

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

977 # Return an intersection as C{Intersector5Tuple}.

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

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

980 if k2:

981 A.set_(k2=X.kA)

982 B.set_(k2=X.kB)

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

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

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

986

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

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

989

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

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

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

993

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

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

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

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

998 the second point (C{degrees}).

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

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

1001

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

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

1004 '''

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

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

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

1008 if name:

1009 gl.name= name

1010 return gl

1011

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

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

1014

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

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

1017

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

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

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

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

1022 (C{bool}).

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

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

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

1026

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

1028 B{C{_C}} above.

1029

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

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

1032 Line} or other C{InverseLine}.

1033 '''

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

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

1036 return self._C(M, glA, glB, _AB=False, **_C) if _C else M

1037

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

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

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

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

1042 sA = glA.Distance()

1043 sB = glB.Distance()

1044 X = XDict_(sA * _0_5, sB * _0_5) # center

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

1046

1047 def Middle5(self, glA, glB, **raiser):

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

1049

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

1051

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

1053 '''

1054 M = self.Middle(glA, glB, _C=True, **raiser)

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

1056 return Middle5Tuple(_llz2G(A, glA, glA), _llz2G(B, glB, glB), s, a, c)

1057

1058 def _m12_M12_M21(self, gl, s):

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

1060 return P.m12, P.M12, P.M21

1061

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

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

1064

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

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

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

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

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

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

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

1072 between the intersections.

1073

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

1075

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

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

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

1079

1080 @raise IntersectionError: No convergence.

1081

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

1083 '''

1084 self._xLines(glA, glB)

1085 if eps1 or _C: # eps

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

1087

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

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

1090 while S_:

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

1092 X = X0._fixCoincident(X)

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

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

1095 if z:

1096 if not c:

1097 continue

1098 Tt_ = self._conjDist3Tt_(c, X0)

1099 else:

1100 Tt_ = (X, t),

1101

1102 for T, t in Tt_:

1103 if t < d or Q is _XINF:

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

1105 i += 1

1106

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

1108 if c else _0t):

1109 T = X

1110 if s and c:

1111 s *= self._D2

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

1113 T._skip(S_, self._T2D2Delta)

1114

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

1116

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

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

1119

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

1121 if none found.

1122

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

1124 '''

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

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

1127

1128 def _obliqDist4(self):

1129 zx = 45.0

1130 if self.f:

1131 _abs, _cjD5 = fabs, self._conjDist5

1132

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

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

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

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

1137 for _ in range(16):

1138 if ds1 == ds0:

1139 break

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

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

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

1143 if _abs(ds1) < dsx:

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

1145 if not dsx:

1146 break

1147 else:

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

1149 return sx, zx, sAx, sBx

1150

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

1152 latx = 64.0

1153 lat = _90_0 - latx

1154 if self.f:

1155 _d, _pD2 = fdot, self._polarDist2

1156

1157 s0, lat0 = _pD2(latx - _1_0)

1158 s1, lat1 = _pD2(latx + _1_0)

1159 s2, lat2 = \

1160 sx, latx = _pD2(latx)

1161 prolate = self.f < 0

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

1163 for _ in range(_TRIPS):

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

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

1166 if not d: # or isnan(d)

1167 break

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

1169 (lat0 + lat2) * s1,

1170 (lat2 + lat1) * s0) / d

1171 s0, lat0 = s1, lat1

1172 s1, lat1 = s2, lat2

1173 s2, lat2 = _pD2(lat)

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

1175 sx, latx = s2, lat2

1176 if lats:

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

1178 sx += sx

1179 else:

1180 sx = self._cHalf

1181 return sx, latx, lat

1182

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

1184 gl = self._Line(lat1=lat1)

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

1186 if lat2:

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

1188 return s, lat1

1189

1190 def _Position(self, gl, s):

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

1192

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

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

1195

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

1197 segments intersect, the intersection is the

1198 one closest to the mid-points of segments.

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

1200 intersections on the segments and specify

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

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

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

1204 are an 4-C{args} L{Line<Intersector.Line>} or

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

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

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

1208 and C{aSB} between the intersections.

1209

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

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

1212 intersection is found, C{None}.

1213

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

1215 invalid, incompatible, ill-configured or

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

1217

1218 @raise IntersectionError: No convergence.

1219

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

1221 '''

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

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

1224 if Q is not None:

1225 if Q.c: # anti-/parallel

1226 Q._fixSegment(sA, sB)

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

1228 d0 = X0.L1(Q)

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

1230 i = Q.iteration

1231 for T in X0._corners(sA, sB, self._T2):

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

1233 X = T._fixCoincident(X)

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

1235 d = X0.L1(X)

1236 if d < d0 or proven is None:

1237 Q, d0 = X, d

1238 if proven is None:

1239 break

1240 Q.set_(iteration=i)

1241

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

1243 return Q

1244

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

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

1247

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

1249 or C{None} if none found.

1250

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

1252 '''

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

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

1255

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

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

1258 '''

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

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

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

1262 D = self._Inverse(A, B)

1263

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

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

1266 c, dc = _diff182(a, b)

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

1268 a, da = -a, -da

1269 b, db = -b, -db

1270 sa, ca = _sincos2de(a, da)

1271 sb, cb = _sincos2de(b, db)

1272

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

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

1275 sA = sB = 0 # at intersection

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

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

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

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

1280 sB = -cb * z * _0_5

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

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

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

1284 else: # general case

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

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

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

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

1289 # order to avoid some convergence failures in _Basic2.

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

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

1292 c = 0

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

1294

1295 @Property_RO

1296 def _T2D1Delta(self):

1297 return self._T2d3Delta(self._D1)

1298

1299 @Property_RO

1300 def _T2D2Delta(self):

1301 return self._T2d3Delta(self._D2)

1302

1303 def _T2d3Delta(self, d3):

1304 return self._T2 - d3 - self.Delta

1305

1306 @Property_RO

1307 def _Tol(self): # convergence tolerance

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

1309

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

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

1312

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

1314 for further details.

1315

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

1317 '''

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

1319

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

1321 # check two geodesic lines vs this geodesic

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

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

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

1325 try:

1326 _xgeodesics(gl.geodesic, self.geodesic)

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

1328 t = gl.geodesic.InverseLine.__name__

1329 raise TypeError(_not_(_an(t)))

1330 c = gl.caps & C

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

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

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

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

1335 except Exception as x:

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

1337

1338

1339 class Intersectool5Tuple(_NamedTuple):

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

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

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

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

1344

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

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

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

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

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

1350 L{XDict} for more details.

1351 '''

1352 _Names_ = (_A_, _B_, _sAB_, 'aAB', _c_)

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

1354

1355

1356 class Intersector5Tuple(Intersectool5Tuple):

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

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

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

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

1361

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

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

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

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

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

1367 '''

1368 _Names_ = Intersectool5Tuple._Names_

1369

1370

1371 class Middle5Tuple(Intersectool5Tuple):

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

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

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

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

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

1377 for more details.

1378 '''

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

1380

1381

1382 class _List(list):

1383

1384 _Delta = 0 # equality margin

1385

1386 def __init__(self, Delta):

1387 self._Delta = Delta

1388# list.__init__(self)

1389

1390 def __contains__(self, other):

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

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

1393 D, _D1 = self._Delta, _L1

1394 for X in self:

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

1396 return True

1397 return False

1398

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

1400 # append an item, updated

1401 if d0_i:

1402 d0, i = d0_i

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

1404 self.append(X)

1405 return X.sX0

1406

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

1408 # trim and sort the X items

1409

1410 def _key(X):

1411 return X.sX0 # rank of X

1412

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

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

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

1416

1417

1418 def _L1(a, b):

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

1420 '''

1421 return fabs(a) + fabs(b)

1422

1423

1424 def _llz2G(G, gl, il):

1425 '''(INTERNAL) Set C{InverseLine} attrs in C{G}, a C{GDict}.

1426 '''

1427 return G.set_(lat1=gl.lat1, lon1=gl.lon1, azi1=il.azi1, a12=il.a13, # .Arc()

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

1429

1430

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

1432

1433 from pygeodesy import printf

1434

1435 def _main(args):

1436

1437 from pygeodesy import GeodesicExact

1438 from pygeodesy.internals import _plural, _usage

1439 from pygeodesy.interns import _COLONSPACE_, _DASH_, _DOT_, \

1440 _EQUAL_, _i_, _n_, _version_, _X_

1441 import re

1442

1443 class XY0(Float):

1444 pass

1445

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

1447 ll4 = ' latA1 lonA1'

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

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

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

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

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

1453 alts=((_c_ + llz2),

1454 (_i_ + ll4),

1455 (_m_ + ll4),

1456 (_M_ + ll4),

1457 (_n_ + llz + ' aziB'),

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

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

1460

1461 def _starts(opt, arg, n=0):

1462 return opt.startswith(arg) and len(arg.lstrip(_DASH_)) > n

1463

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

1465

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

1467 M = m = _R = None

1468 _T = _V = _h = _C = False

1469

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

1471 arg = args.pop(0)

1472 if arg == _C__ or _starts('--Check', arg):

1473 _C = True

1474 elif arg == _c__:

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

1476 elif arg == '-h' or _starts('--help', arg):

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

1478 elif arg == _i__:

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

1480 elif arg == '-m':

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

1482 _R = None

1483 elif arg == '-M':

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

1485 _R = _C = None

1486 elif arg == _n__:

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

1488 elif arg == _o__:

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

1490 elif arg == _R__ and args:

1491 _R = args.pop(0)

1492 elif arg == '-T' or _starts('--Tool', arg):

1493 I = Intersectool() # PYCHOK I

1494 if _V:

1495 I.verbose = True

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

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

1498 M, _T = None, True

1499 elif arg == '-V' or _starts('--Verbose', arg):

1500 _V = True

1501 if _T:

1502 I.verbose = True

1503 elif arg == '-v' or _starts('--version', arg):

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

1505 (__version__, repr(I)))))

1506 else:

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

1508

1509 if _h or M is None:

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

1511# exit(0)

1512 else:

1513 n = len(args)

1514 if n < m:

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

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

1517 args[:] = args[:m]

1518

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

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

1521 # get latA lonA aziA latA lonA aziB

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

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

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

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

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

1527 if _R:

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

1529 kwds.update(sMaX0=m)

1530 M = I.All

1531

1532 n = len(args) // 2

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

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

1535

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

1537 if _V:

1538 X = _SPACE_(_X_, _EQUAL_, m)

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

1540

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

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

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

1544 else:

1545 for i, X in enumerate(X):

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

1547

1548 from sys import argv, stderr

1549 try:

1550 if len(argv) == 2 and argv[1] == _HASH_:

1551 from pygeodesy.internals import _usage_argv

1552

1553 s = _SPACE_(*_usage_argv(__file__))

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

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

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

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

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

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

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

1561 '# -631414 5988887 0 -3',

1562 '# -4.05187 -4.00000 -4.05187 -4.00000 0',

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

1564 '-M 0 0 10 10 50 -4 50S 4W',

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

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

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

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

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

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

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

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

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

1574 if t.startswith(_HASH_):

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

1576 else:

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

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

1579 _main(t)

1580 else:

1581 _main(argv[1:])

1582

1583 except Exception as x:

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

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

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

1587 raise

1588 exit(1)

1589

1590# % env PYGEODESY_INTERSECTTOOL=... python3 -m pygeodesy.geodesici -h

1591#

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

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

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

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

1596# -M latA1 lonA1 latA2 lonA2 latB1 lonB1 latB2 lonB2 |

1597# -n latA lonA aziA aziB |

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

1599

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

1601#

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

1603

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

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

1606

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

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

1609

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

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

1612

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

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

1615

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

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

1618

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

1620# -631414 5988887 0 -3

1621# -4.05187 -4.00000 -4.05187 -4.00000 0

1622

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

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

1625

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

1627# Intersector.Middle5: Middle5Tuple(A=GDict(a12=14.106434, a1M=7.053396, azi1=44.75191, azi2=45.629037, aziM=44.969535, lat1=0.0, lat2=10.0, latM=5.019509, lon1=0.0, lon2=10.0, lonM=4.961883, s12=1565109.099218, s1M=782554.549609), B=GDict(a12=99.810444, a1M=49.869061, azi1=180.0, azi2=180.0, aziM=180.0, lat1=50.0, lat2=-50.0, latM=0.036282, lon1=-4.0, lon2=-4.0, lonM=-4.0, s12=11073670.322999, s1M=5536835.161499), sMM=1138574.546746, aMM=10.262308, c=0)

1628

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

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

1631

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

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

1634

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

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

1637

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

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

1640

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

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

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

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

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

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

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

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

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

1650

1651# % python3 -m ....pygeodesy.geodesici -C -R 4e7 -c 50 -4 -147.7 0 0 90

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

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

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

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

1656

1657# % env PYGEODESY_INTERSECTTOOL=... python3 -m ....pygeodesy.geodesici -T -C -R 4e7 -c 50 -4 -147.7 0 0 90

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

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

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

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

1662

1663# % python3 -m ....pygeodesy.geodesici -R 4e7 -i 0 0 10 10 50 -4 -50 -4

1664# Intersector.All[0]: XDict(c=0, sA=-631414.26877, sB=5988887.278435, sX0=1866020.935315)

1665# Intersector.All[1]: XDict(c=0, sA=19422725.117572, sB=-14062417.105648, sX0=38239422.83511)

1666# Intersector.All[2]: XDict(c=0, sA=19422725.117572, sB=25945445.811603, sX0=39048781.218067)

1667# Intersector.All[3]: XDict(c=0, sA=39476927.464575, sB=5894074.699478, sX0=39051612.452944)

1668

1669# % env PYGEODESY_INTERSECTTOOL=... python3 -m ....pygeodesy.geodesici -T -R 4e7 -i 0 0 10 10 50 -4 -50 -4

1670# Intersectool.All[0]: XDict(c=0, sA=-631414.26877, sB=5988887.278435, sX0=1862009.05513)

1671# Intersectool.All[1]: XDict(c=0, sA=19422725.117572, sB=-14062417.105648, sX0=38243434.715295)

1672# Intersectool.All[2]: XDict(c=0, sA=19422725.117572, sB=25945445.811603, sX0=39044769.337882)

1673# Intersectool.All[3]: XDict(c=0, sA=39476927.464575, sB=5894074.699478, sX0=39047600.57276)

1674

1675

1676# **) MIT License

1677#

1678# Copyright (C) 2024-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1679#

1680# Permission is hereby granted, free of charge, to any person obtaining a

1681# copy of this software and associated documentation files (the "Software"),

1682# to deal in the Software without restriction, including without limitation

1683# the rights to use, copy, modify, merge, publish, distribute, sublicense,

1684# and/or sell copies of the Software, and to permit persons to whom the

1685# Software is furnished to do so, subject to the following conditions:

1686#

1687# The above copyright notice and this permission notice shall be included

1688# in all copies or substantial portions of the Software.

1689#

1690# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

1691# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

1692# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

1693# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

1694# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

1695# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

1696# OTHER DEALINGS IN THE SOFTWARE.