Coverage for pygeodesy/ellipsoidalBaseDI.py: 96%

312 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-11 14:35 -0400

1 

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

3 

4u'''(INTERNAL) Private, ellipsoidal Direct/Inverse geodesy base 

5class C{LatLonEllipsoidalBaseDI} and functions. 

6''' 

7# make sure int/int division yields float quotient, see .basics 

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

9 

10from pygeodesy.basics import isscalar, issubclassof 

11from pygeodesy.constants import EPS, MAX, PI, PI2, PI_4, isnear0, isnear1, \ 

12 _0_0, _0_5, _1_5, _3_0 

13from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase, Property_RO, \ 

14 property_RO, _TOL_M 

15from pygeodesy.errors import _AssertionError, IntersectionError, _IsnotError, \ 

16 _or, _ValueError, _xellipsoidal, _xError, _xkwds_not 

17from pygeodesy.fmath import favg, fmean_ 

18from pygeodesy.fsums import Fmt, fsum_ 

19from pygeodesy.formy import opposing, _radical2 

20from pygeodesy.interns import _antipodal_, _concentric_, _exceed_PI_radians_, \ 

21 _near_, _SPACE_, _too_ 

22from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS 

23from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \ 

24 Intersection3Tuple, NearestOn2Tuple, \ 

25 NearestOn8Tuple, _LL4Tuple 

26# from pygeodesy.props import Property_RO, property_RO # from .ellipsoidalBase 

27# from pygeodesy.streprs import Fmt # from .fsums 

28from pygeodesy.units import _fi_j2, Height, Radius_, Scalar 

29from pygeodesy.utily import m2km, unroll180, _unrollon, wrap90, wrap180, wrap360 

30 

31from math import degrees, radians 

32 

33__all__ = _ALL_LAZY.ellipsoidalBaseDI 

34__version__ = '23.04.11' 

35 

36_polar__ = 'polar?' 

37_too_low_ = _too_('low') 

38_B2END = _1_5 # _intersect3 bearing to pseudo-end point factor 

39_TRIPS = 33 # _intersect3, _intersects2, _nearestOn interations, 6..9 sufficient? 

40 

41 

42class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase): 

43 '''(INTERNAL) Base class for C{ellipsoidal*.LatLon} classes 

44 with I{overloaded} C{Direct} and C{Inverse} methods. 

45 ''' 

46 

47 def bearingTo2(self, other, wrap=False): 

48 '''Compute the initial and final bearing (forward and reverse 

49 azimuth) from this to an other point, using this C{Inverse} 

50 method. See methods L{initialBearingTo} and L{finalBearingTo} 

51 for more details. 

52 

53 @arg other: The other point (C{LatLon}). 

54 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

55 

56 @return: A L{Bearing2Tuple}C{(initial, final)}. 

57 

58 @raise TypeError: The B{C{other}} point is not C{LatLon}. 

59 

60 @raise ValueError: If this and the B{C{other}} point's L{Datum} 

61 ellipsoids are not compatible. 

62 ''' 

63 r = self._Inverse(other, wrap) 

64 return Bearing2Tuple(r.initial, r.final, name=self.name) 

65 

66 def destination(self, distance, bearing, height=None): 

67 '''Compute the destination point after having travelled for 

68 the given distance from this point along a geodesic given 

69 by an initial bearing, using this C{Direct} method. See 

70 method L{destination2} for more details. 

71 

72 @arg distance: Distance (C{meter}). 

73 @arg bearing: Initial bearing in (compass C{degrees360}). 

74 @kwarg height: Optional height, overriding the default 

75 height (C{meter}, same units as C{distance}). 

76 

77 @return: The destination point (C{LatLon}). 

78 ''' 

79 return self._Direct(distance, bearing, self.classof, height).destination 

80 

81 def destination2(self, distance, bearing, height=None): 

82 '''Compute the destination point and the final bearing (reverse 

83 azimuth) after having travelled for the given distance from 

84 this point along a geodesic given by an initial bearing, 

85 using this C{Direct} method. 

86 

87 The distance must be in the same units as this point's datum 

88 axes, conventionally C{meter}. The distance is measured on 

89 the surface of the ellipsoid, ignoring this point's height. 

90 

91 The initial and final bearing (forward and reverse azimuth) 

92 are in compass C{degrees360}. 

93 

94 The destination point's height and datum are set to this 

95 point's height and datum, unless the former is overridden. 

96 

97 @arg distance: Distance (C{meter}). 

98 @arg bearing: Initial bearing (compass C{degrees360}). 

99 @kwarg height: Optional height, overriding the default 

100 height (C{meter}, same units as C{distance}). 

101 

102 @return: A L{Destination2Tuple}C{(destination, final)}. 

103 ''' 

104 r = self._Direct(distance, bearing, self.classof, height) 

105 return self._xnamed(r) 

106 

107 def _Direct(self, distance, bearing, LL, height): # overloaded by I{Vincenty} 

108 '''(INTERNAL) I{Karney}'s C{Direct} method. 

109 

110 @return: A L{Destination2Tuple}C{(destination, final)} or 

111 a L{Destination3Tuple}C{(lat, lon, final)} if 

112 B{C{LL}} is C{None}. 

113 ''' 

114 g = self.geodesic 

115 r = g.Direct3(self.lat, self.lon, bearing, distance) 

116 if LL: 

117 r = self._Direct2Tuple(LL, height, r) 

118 return r 

119 

120 def _Direct2Tuple(self, LL, height, r): 

121 '''(INTERNAL) Helper for C{._Direct} result L{Destination2Tuple}. 

122 ''' 

123 h = self.height if height is None else height 

124 d = LL(wrap90(r.lat), wrap180(r.lon), height=h, datum=self.datum, name=self.name, 

125 **_xkwds_not(None, epoch=self.epoch, reframe=self.reframe)) 

126 return Destination2Tuple(d, wrap360(r.final)) 

127 

128 def distanceTo(self, other, wrap=False, **unused): # ignore radius=R_M 

129 '''Compute the distance between this and an other point 

130 along a geodesic, using this C{Inverse} method. See method 

131 L{distanceTo3} for more details. 

132 

133 @arg other: The other point (C{LatLon}). 

134 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

135 

136 @return: Distance (C{meter}). 

137 

138 @raise TypeError: The B{C{other}} point is not C{LatLon}. 

139 

140 @raise ValueError: If this and the B{C{other}} point's L{Datum} 

141 ellipsoids are not compatible. 

142 ''' 

143 return self._Inverse(other, wrap, azis=False).distance 

144 

145 def distanceTo3(self, other, wrap=False): 

146 '''Compute the distance, the initial and final bearing along 

147 a geodesic between this and an other point, using this 

148 C{Inverse} method. 

149 

150 The distance is in the same units as this point's datum axes, 

151 conventionally meter. The distance is measured on the surface 

152 of the ellipsoid, ignoring this point's height. 

153 

154 The initial and final bearing (forward and reverse azimuth) 

155 are in compass C{degrees360} from North. 

156 

157 @arg other: Destination point (C{LatLon}). 

158 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

159 

160 @return: A L{Distance3Tuple}C{(distance, initial, final)}. 

161 

162 @raise TypeError: The B{C{other}} point is not C{LatLon}. 

163 

164 @raise ValueError: If this and the B{C{other}} point's L{Datum} 

165 ellipsoids are not compatible. 

166 ''' 

167 return self._xnamed(self._Inverse(other, wrap)) 

168 

169 def finalBearingOn(self, distance, bearing): 

170 '''Compute the final bearing (reverse azimuth) after having 

171 travelled for the given distance along a geodesic given by 

172 an initial bearing from this point, using this C{Direct} 

173 method. See method L{destination2} for more details. 

174 

175 @arg distance: Distance (C{meter}). 

176 @arg bearing: Initial bearing (compass C{degrees360}). 

177 

178 @return: Final bearing (compass C{degrees360}). 

179 ''' 

180 return self._Direct(distance, bearing, None, None).final 

181 

182 def finalBearingTo(self, other, wrap=False): 

183 '''Compute the final bearing (reverse azimuth) after having 

184 travelled along a geodesic from this point to an other 

185 point, using this C{Inverse} method. See method 

186 L{distanceTo3} for more details. 

187 

188 @arg other: The other point (C{LatLon}). 

189 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

190 

191 @return: Final bearing (compass C{degrees360}). 

192 

193 @raise TypeError: The B{C{other}} point is not C{LatLon}. 

194 

195 @raise ValueError: If this and the B{C{other}} point's L{Datum} 

196 ellipsoids are not compatible. 

197 ''' 

198 return self._Inverse(other, wrap).final 

199 

200 @Property_RO 

201 def geodesic(self): # overloaded by I{Karney}'s, N/A for I{Vincenty} 

202 '''N/A, invalid (C{None} I{always}). 

203 ''' 

204 return None # PYCHOK no cover 

205 

206 def initialBearingTo(self, other, wrap=False): 

207 '''Compute the initial bearing (forward azimuth) to travel 

208 along a geodesic from this point to an other point, 

209 using this C{Inverse} method. See method L{distanceTo3} 

210 for more details. 

211 

212 @arg other: The other point (C{LatLon}). 

213 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

214 

215 @return: Initial bearing (compass C{degrees360}). 

216 

217 @raise TypeError: The B{C{other}} point is not C{LatLon}. 

218 

219 @raise ValueError: If this and the B{C{other}} point's L{Datum} 

220 ellipsoids are not compatible. 

221 ''' 

222 return self._Inverse(other, wrap).initial 

223 

224 def intermediateTo(self, other, fraction, height=None, wrap=False): 

225 '''Return the point at given fraction along the geodesic between 

226 this and an other point, using this C{Direct} and C{Inverse} 

227 methods. 

228 

229 @arg other: The other point (C{LatLon}). 

230 @arg fraction: Fraction between both points (C{scalar}, 

231 0.0 at this and 1.0 at the other point. 

232 @kwarg height: Optional height, overriding the fractional 

233 height (C{meter}). 

234 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

235 

236 @return: Intermediate point (C{LatLon}). 

237 

238 @raise TypeError: The B{C{other}} point is not C{LatLon}. 

239 

240 @raise UnitError: Invalid B{C{fraction}} or B{C{height}}. 

241 

242 @raise ValueError: If this and the B{C{other}} point's L{Datum} 

243 ellipsoids are not compatible. 

244 

245 @see: Methods L{distanceTo3}, L{destination}, C{midpointTo} and 

246 C{rhumbMidpointTo}. 

247 ''' 

248 f = Scalar(fraction=fraction) 

249 if isnear0(f): 

250 r = self 

251 elif isnear1(f) and not wrap: 

252 r = self.others(other) 

253 else: # negative fraction OK 

254 t = self.distanceTo3(other, wrap=wrap) 

255 h = self._havg(other, f=f) if height is None else Height(height) 

256 r = self.destination(t.distance * f, t.initial, height=h) 

257 return r 

258 

259 def _Inverse(self, other, wrap, **unused): # azis=False, overloaded by I{Vincenty} 

260 '''(INTERNAL) I{Karney}'s C{Inverse} method. 

261 

262 @return: A L{Distance3Tuple}C{(distance, initial, final)}. 

263 

264 @raise TypeError: The B{C{other}} point is not L{LatLon}. 

265 

266 @raise ValueError: If this and the B{C{other}} point's 

267 L{Datum} ellipsoids are not compatible. 

268 ''' 

269 _ = self.ellipsoids(other) 

270 g = self.geodesic 

271 _, lon = unroll180(self.lon, other.lon, wrap=wrap) 

272 return g.Inverse3(self.lat, self.lon, other.lat, lon) 

273 

274 def nearestOn8(self, points, closed=False, height=None, wrap=False, 

275 equidistant=None, tol=_TOL_M): 

276 '''Iteratively locate the point on a path or polygon closest 

277 to this point. 

278 

279 @arg points: The path or polygon points (C{LatLon}[]). 

280 @kwarg closed: Optionally, close the polygon (C{bool}). 

281 @kwarg height: Optional height, overriding the height of this and 

282 all other points (C{meter}, conventionally). If 

283 B{C{height}} is C{None}, the height of each point 

284 is taken into account for distances. 

285 

286 @return: A L{NearestOn8Tuple}C{(closest, distance, fi, j, start, 

287 end, initial, final)} with C{distance} in C{meter}, 

288 conventionally and with the C{closest}, the C{start} 

289 the C{end} point each an instance of this C{LatLon}. 

290 

291 @raise PointsError: Insufficient number of B{C{points}}. 

292 

293 @raise TypeError: Some B{C{points}} or B{C{equidistant}} invalid. 

294 

295 @raise ValueError: Some B{C{points}}' datum or ellipsoid incompatible 

296 or no convergence for the given B{C{tol}}. 

297 

298 @see: Function L{pygeodesy.nearestOn6} and method C{nearestOn6}. 

299 ''' 

300 D3 = self.distanceTo3 # Distance3Tuple 

301 

302 try: 

303 Ps = self.PointsIter(points, loop=1) 

304 p1 = c = s = e = Ps[0] 

305 _ = self.ellipsoids(p1) 

306 c3 = D3(c, wrap=wrap) # XXX wrap=False? 

307 

308 except (TypeError, ValueError) as x: 

309 raise _xError(x, Fmt.SQUARE(points=0), p1, this=self, tol=tol, 

310 closed=closed, height=height, wrap=wrap) 

311 

312 # get the azimuthal equidistant projection, once 

313 A = _Equidistant00(equidistant, c) 

314 b = _Box(c, c3.distance) 

315 m = f = i = 0 # p1..p2 == points[i]..[j] 

316 

317 kwds = dict(within=True, height=height, tol=tol, 

318 LatLon=self.classof, # this LatLon 

319 datum=self.datum, epoch=self.epoch, reframe=self.reframe) 

320 try: 

321 for j, p2 in Ps.enumerate(closed=closed): 

322 if wrap and j != 0: 

323 p2 = _unrollon(p1, p2) 

324 # skip edge if no overlap with box around closest 

325 if j < 4 or b.overlaps(p1.lat, p1.lon, p2.lat, p2.lon): 

326 p, t, _ = _nearestOn3_(self, p1, p2, A, **kwds) 

327 d3 = D3(p, wrap=False) # already unrolled 

328 if d3.distance < c3.distance: 

329 c3, c, s, e, f = d3, p, p1, p2, (i + t) 

330 b = _Box(c, c3.distance) 

331 m = max(m, c.iteration) 

332 p1, i = p2, j 

333 

334 except (TypeError, ValueError) as x: 

335 raise _xError(x, Fmt.SQUARE(points=i), p1, 

336 Fmt.SQUARE(points=j), p2, this=self, tol=tol, 

337 closed=closed, height=height, wrap=wrap) 

338 

339 f, j = _fi_j2(f, len(Ps)) # like .vector3d.nearestOn6 

340 

341 n = self.nearestOn8.__name__ 

342 c.rename(n) 

343 if s is not c: 

344 s = s.copy(name=n) 

345 if e is not c: 

346 e = e.copy(name=n) 

347 return NearestOn8Tuple(c, c3.distance, f, j, s, e, c3.initial, c3.final, 

348 iteration=m) # ._iteration for tests 

349 

350 

351class _Box(object): 

352 '''Bounding box around a C{LatLon} point. 

353 

354 @see: Function C{_box4} in .clipy.py. 

355 ''' 

356 def __init__(self, center, distance): 

357 '''New L{_Box} around point. 

358 

359 @arg center: The center point (C{LatLon}). 

360 @arg distance: Radius, half-size of the box 

361 (C{meter}, conventionally) 

362 ''' 

363 E = center.ellipsoid() 

364 d = degrees(distance / max(E.a, E.b)) + _0_5 # some margin 

365 self._N = center.lat + d 

366 self._S = center.lat - d 

367 self._E = center.lon + d 

368 self._W = center.lon - d 

369 

370 def overlaps(self, lat1, lon1, lat2, lon2): 

371 '''Check whether this box overlaps a line between 2 points. 

372 

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

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

375 @arg lat2: Latitude of second point (C{degrees}). 

376 @arg lon2: Longitude of second point (C{degrees}). 

377 

378 @return: C{False} if there is certainly no overlap, 

379 C{True} otherwise (C{bool}). 

380 ''' 

381 non_ = ((lat1 > self._N or lat2 < self._S) if lat1 < lat2 else 

382 (lat2 > self._N or lat1 < self._S)) or \ 

383 ((lon1 > self._E or lon2 < self._W) if lon1 < lon2 else 

384 (lon2 > self._E or lon1 < self._W)) 

385 return not non_ 

386 

387 

388class _Tol(object): 

389 '''Handle a tolerance in C{meter} as C{degrees} and C{meter}. 

390 ''' 

391 _deg = 0 # tol in degrees 

392 _lat = 0 

393 _m = 0 # tol in meter 

394 _min = MAX # degrees 

395 _prev = None 

396 _r = 0 

397 

398 def __init__(self, tol_m, E, lat, *lats): 

399 '''New L{_Tol}. 

400 

401 @arg tol_m: Tolerance (C{meter}, only). 

402 @arg E: Earth ellipsoid (L{Ellipsoid}). 

403 @arg lat: Latitude (C{degrees}). 

404 @arg lats: Additional latitudes (C{degrees}). 

405 ''' 

406 if lats: 

407 lat = fmean_(lat, *lats) 

408 self._lat = lat 

409 self._r = max(EPS, E.rocMean(lat)) 

410 self._m = max(EPS, tol_m) 

411 self._deg = max(EPS, degrees(self._m / self._r)) # avoid m2degrees! 

412 

413 @property_RO 

414 def degrees(self): 

415 '''Get this tolerance in C{degrees}. 

416 ''' 

417 return self._deg 

418 

419 def degrees2m(self, deg): 

420 '''Convert B{C{deg}} to meter at the same C{lat} and earth radius. 

421 ''' 

422 return self.radius * radians(deg) / PI2 # avoid degrees2m! 

423 

424 def degError(self, Error=_ValueError): 

425 '''Compose an error with the C{deg}rees minimum. 

426 ''' 

427 return self.mError(self.degrees2m(self._min), Error=Error) 

428 

429 def done(self, deg): 

430 '''Check C{deg} vs. tolerance and previous value. 

431 ''' 

432 if deg < self._deg or deg == self._prev: 

433 return True 

434 self._min = min(self._min, deg) 

435 self._prev = deg 

436 return False 

437 

438 @property_RO 

439 def lat(self): 

440 '''Get the mean latitude in C{degrees}. 

441 ''' 

442 return self._lat 

443 

444 def mError(self, m, Error=_ValueError): 

445 '''Compose an error with B{C{m}}eter minimum. 

446 ''' 

447 t = _SPACE_(Fmt.tolerance(self.meter), _too_low_) 

448 if m2km(m) > self.meter: 

449 t = _or(t, _antipodal_, _near_(_polar__)) 

450 return Error(Fmt.no_convergence(m), txt=t) 

451 

452 @property_RO 

453 def meter(self): 

454 '''Get this tolerance in C{meter}. 

455 ''' 

456 return self._m 

457 

458 @property_RO 

459 def radius(self): 

460 '''Get the earth radius in C{meter}. 

461 ''' 

462 return self._r 

463 

464 def reset(self): 

465 '''Reset tolerances. 

466 ''' 

467 self._min = MAX 

468 self._prev = None 

469 

470 

471def _Equidistant00(equidistant, p1): 

472 '''(INTERNAL) Get an C{Equidistant*(0, 0, ...)} instance. 

473 ''' 

474 if equidistant is None or not callable(equidistant): 

475 equidistant = p1.Equidistant 

476 elif not issubclassof(equidistant, *_MODS.azimuthal._Equidistants): # PYCHOK no cover 

477 t = tuple(_.__name__ for _ in _MODS.azimuthal._Equidistants) 

478 raise _IsnotError(*t, equidistant=equidistant) 

479 return equidistant(0, 0, p1.datum) 

480 

481 

482def _intersect3(s1, end1, s2, end2, height=None, wrap=True, # MCCABE 16 

483 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds): 

484 '''(INTERNAL) Intersect two (ellipsoidal) lines, see ellipsoidal method 

485 L{intersection3}, separated to allow callers to embellish any exceptions. 

486 ''' 

487 _LLS = _MODS.sphericalTrigonometry.LatLon 

488 _si = _MODS.sphericalTrigonometry._intersect 

489 _vi3 = _MODS.vector3d._intersect3d3 

490 

491 def _b_d(s, e, w, t, h=_0_0): 

492 # compute opposing and distance 

493 t = s.classof(t.lat, t.lon, height=h, name=t.name) 

494 t = s.distanceTo3(t, wrap=w) # Distance3Tuple 

495 b = opposing(e, t.initial) # "before" start 

496 return b, t.distance 

497 

498 def _b_e(s, e, w, t): 

499 # compute an end point along the initial bearing 

500 # about 1.5 times the distance to the gu-/estimate, at 

501 # least 1/8 and at most 3/8 of the earth perimeter like 

502 # radians in .sphericalTrigonometry._int3d2 and bearing 

503 # comparison in .sphericalTrigonometry._intb 

504 b, d = _b_d(s, e, w, t, h=t.height) 

505 m = s.ellipsoid().R2x * PI_4 # authalic exact 

506 d = min(max(d * _B2END, m), m * _3_0) 

507 e = s.destination(d, e) 

508 return b, (_unrollon(s, e) if w else e) 

509 

510 def _e_ll(s, e, w, **end): 

511 # return 2-tuple (end, False if bearing else True) 

512 ll = not isscalar(e) 

513 if ll: 

514 e = s.others(**end) 

515 if w: # unroll180 == .karney._unroll2 

516 e = _unrollon(s, e) 

517 return e, ll 

518 

519 def _llS(ll): # return an C{_LLS} instance 

520 return _LLS(ll.lat, ll.lon, height=ll.height) 

521 

522 def _o(o, b, n, s, t, e): 

523 # determine C{o}utside before, on or after start point 

524 if not o: # intersection may be on start 

525 if _MODS.latlonBase._isequalTo(s, t, eps=e.degrees): 

526 return o 

527 return -n if b else n 

528 

529 E = s1.ellipsoids(s2) 

530 

531 e1, ll1 = _e_ll(s1, end1, wrap, end1=end1) 

532 e2, ll2 = _e_ll(s2, end2, wrap, end2=end2) 

533 

534 e = _Tol(tol, E, s1.lat, (e1.lat if ll1 else s1.lat), 

535 s2.lat, (e2.lat if ll2 else s2.lat)) 

536 

537 # get the azimuthal equidistant projection 

538 A = _Equidistant00(equidistant, s1) 

539 

540 # gu-/estimate initial intersection, spherically ... 

541 t = _si(_llS(s1), (_llS(e1) if ll1 else e1), 

542 _llS(s2), (_llS(e2) if ll2 else e2), 

543 height=height, wrap=False, LatLon=_LLS) # unrolled already 

544 h, n = t.height, t.name 

545 

546 if not ll1: 

547 b1, e1 = _b_e(s1, e1, wrap, t) 

548 if not ll2: 

549 b2, e2 = _b_e(s2, e2, wrap, t) 

550 

551 # ... and iterate as Karney describes, for references 

552 # @see: Function L{ellipsoidalKarney.intersection3}. 

553 for i in range(1, _TRIPS): 

554 A.reset(t.lat, t.lon) # gu-/estimate as origin 

555 # convert start and end points to projection 

556 # space and compute an intersection there 

557 v, o1, o2 = _vi3(*A._forwards(s1, e1, s2, e2), 

558 eps=e.meter, useZ=False) 

559 # convert intersection back to geodetic 

560 t, d = A._reverse2(v) 

561 if e.done(d): # below tol or unchanged? 

562 t._iteration = i 

563 break 

564 else: 

565 raise e.degError(Error=IntersectionError) 

566 

567 # like .sphericalTrigonometry._intersect, if this intersection 

568 # is "before" the first point, use the antipodal intersection 

569 if not (ll1 or ll2): # end1 and end2 are an initial bearing 

570 b1, _ = _b_d(s1, end1, wrap, t) 

571 if b1: 

572 t = t.antipodal() 

573 b1 = not b1 

574 b2, _ = _b_d(s2, end2, wrap, t) 

575 

576 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=s1, 

577 iteration=t._iteration, name=n) 

578 return Intersection3Tuple(r, (o1 if ll1 else _o(o1, b1, 1, s1, t, e)), 

579 (o2 if ll2 else _o(o2, b2, 2, s2, t, e))) 

580 

581 

582def _intersection3(start1, end1, start2, end2, height=None, wrap=True, 

583 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds): 

584 '''(INTERNAL) Iteratively compute the intersection point of two lines, 

585 each defined by two (ellipsoidal) points or an (ellipsoidal) start 

586 point and an initial bearing from North. 

587 ''' 

588 s1 = _xellipsoidal(start1=start1) 

589 s2 = s1.others(start2=start2) 

590 try: 

591 return _intersect3(s1, end1, s2, end2, height=height, wrap=wrap, 

592 equidistant=equidistant, tol=tol, 

593 LatLon=LatLon, **LatLon_kwds) 

594 except (TypeError, ValueError) as x: 

595 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2) 

596 

597 

598def _intersections2(center1, radius1, center2, radius2, height=None, wrap=True, 

599 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds): 

600 '''(INTERNAL) Iteratively compute the intersection points of two circles, 

601 each defined by an (ellipsoidal) center point and a radius. 

602 ''' 

603 c1 = _xellipsoidal(center1=center1) 

604 c2 = c1.others(center2=center2) 

605 try: 

606 return _intersects2(c1, radius1, c2, radius2, height=height, wrap=wrap, 

607 equidistant=equidistant, tol=tol, 

608 LatLon=LatLon, **LatLon_kwds) 

609 except (TypeError, ValueError) as x: 

610 raise _xError(x, center1=center1, radius1=radius1, 

611 center2=center2, radius2=radius2) 

612 

613 

614def _intersects2(c1, radius1, c2, radius2, height=None, wrap=True, # MCCABE 16 

615 equidistant=None, tol=_TOL_M, LatLon=None, **LatLon_kwds): 

616 '''(INTERNAL) Intersect two (ellipsoidal) circles, see L{_intersections2} 

617 above, separated to allow callers to embellish any exceptions. 

618 ''' 

619 _LLS = _MODS.sphericalTrigonometry.LatLon 

620 _si2 = _MODS.sphericalTrigonometry._intersects2 

621 _vi2 = _MODS.vector3d._intersects2 

622 

623 def _latlon4(t, h, n, c): 

624 return _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=c, 

625 iteration=t.iteration, name=n) 

626 

627 def _llS(ll): # return an C{_LLS} instance 

628 return _LLS(ll.lat, ll.lon, height=ll.height) 

629 

630 r1 = Radius_(radius1=radius1) 

631 r2 = Radius_(radius2=radius2) 

632 

633 E = c1.ellipsoids(c2) 

634 # get the azimuthal equidistant projection 

635 A = _Equidistant00(equidistant, c1) 

636 

637 if r1 < r2: 

638 c1, c2 = c2, c1 

639 r1, r2 = r2, r1 

640 

641 if r1 > (min(E.b, E.a) * PI): 

642 raise _ValueError(_exceed_PI_radians_) 

643 

644 if wrap: # unroll180 == .karney._unroll2 

645 c2 = _unrollon(c1, c2) 

646 

647 # distance between centers and radii are 

648 # measured along the ellipsoid's surface 

649 m = c1.distanceTo(c2, wrap=False) # meter 

650 if m < max(r1 - r2, EPS): 

651 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError? 

652 if fsum_(r1, r2, -m) < 0: 

653 raise IntersectionError(_too_(Fmt.distant(m))) 

654 

655 f = _radical2(m, r1, r2).ratio # "radical fraction" 

656 e = _Tol(tol, E, favg(c1.lat, c2.lat, f=f)) 

657 

658 # gu-/estimate initial intersections, spherically ... 

659 t1, t2 = _si2(_llS(c1), r1, _llS(c2), r2, radius=e.radius, 

660 height=height, wrap=False, too_d=m) # unrolled already 

661 h, n = t1.height, t1.name 

662 

663 # ... and iterate as Karney describes, for references 

664 # @see: Function L{ellipsoidalKarney.intersections2}. 

665 ts, ta = [], None 

666 for t in ((t1,) if t1 is t2 else (t1, t2)): 

667 for i in range(1, _TRIPS): 

668 A.reset(t.lat, t.lon) # gu-/estimate as origin 

669 # convert centers to projection space and 

670 # compute the intersections there 

671 t1, t2 = A._forwards(c1, c2) 

672 v1, v2 = _vi2(t1, r1, # XXX * t1.scale?, 

673 t2, r2, # XXX * t2.scale?, 

674 sphere=False, too_d=m) 

675 # convert intersections back to geodetic 

676 if v1 is v2: # abutting 

677 t, d = A._reverse2(v1) 

678 else: # consider the closer intersection 

679 t1, d1 = A._reverse2(v1) 

680 t2, d2 = A._reverse2(v2) 

681 t, d = (t1, d1) if d1 < d2 else (t2, d2) 

682 if e.done(d): # below tol or unchanged? 

683 t._iteration = i # _NamedTuple._iteration 

684 ts.append(t) 

685 if v1 is v2: # abutting 

686 ta = t 

687 break 

688 else: 

689 raise e.degError(Error=IntersectionError) 

690 e.reset() 

691 

692 if ta: # abutting circles 

693 pass # PYCHOK no cover 

694 elif len(ts) == 2: 

695 return (_latlon4(ts[0], h, n, c1), 

696 _latlon4(ts[1], h, n, c2)) 

697 elif len(ts) == 1: # PYCHOK no cover 

698 ta = ts[0] # assume abutting 

699 else: # PYCHOK no cover 

700 raise _AssertionError(ts=ts) 

701 r = _latlon4(ta, h, n, c1) 

702 return r, r 

703 

704 

705def _nearestOn2(p, point1, point2, within=True, height=None, wrap=True, 

706 equidistant=None, tol=_TOL_M, **LatLon_and_kwds): 

707 '''(INTERNAL) Closest point and fraction, like L{_intersects2} above, 

708 separated to allow callers to embellish any exceptions. 

709 ''' 

710 p1 = p.others(point1=point1) 

711 p2 = p.others(point2=point2) 

712 

713 _ = p.ellipsoids(p1) 

714# E = p.ellipsoids(p2) # done in __nearestOn2__ 

715 

716 # get the azimuthal equidistant projection 

717 A = _Equidistant00(equidistant, p) 

718 

719 if wrap: 

720 p1 = _unrollon(p, p1) # XXX do not unroll? 

721 p2 = _unrollon(p, p2) # XXX do not unroll? 

722 p2 = _unrollon(p1, p2) 

723 

724 r, f, e = _nearestOn3_(p, p1, p2, A, within=within, height=height, 

725 tol=tol, **LatLon_and_kwds) 

726 return NearestOn2Tuple(r, f) 

727 

728 

729def _nearestOn3_(p, p1, p2, A, within=True, height=None, tol=_TOL_M, 

730 LatLon=None, **LatLon_kwds): 

731 # Only in function C{_nearestOn2} and method C{nearestOn8} above 

732 _LLS = _MODS.sphericalNvector.LatLon 

733 _V3d = _MODS.vector3d.Vector3d 

734 _vnOn2 = _MODS.vector3d._nearestOn2 

735 

736 def _llS(ll): # return an C{_LLS} instance 

737 return _LLS(ll.lat, ll.lon, height=ll.height) 

738 

739 E = p.ellipsoids(p2) 

740 e = _Tol(tol, E, p.lat, p1.lat, p2.lat) 

741 

742 # gu-/estimate initial nearestOn, spherically ... wrap=False, only! 

743 # using sphericalNvector.LatLon.nearestOn for within=False support 

744 t = _llS(p).nearestOn(_llS(p1), _llS(p2), within=within, 

745 height=height) 

746 n, h = t.name, t.height 

747 if height is None: 

748 h1 = p1.height # use heights as pseudo-Z in projection space 

749 h2 = p2.height # to be included in the closest function 

750 h0 = favg(h1, h2) 

751 else: # ignore heights in distances, Z=0 

752 h0 = h1 = h2 = _0_0 

753 

754 # ... and iterate to find the closest (to the origin with .z 

755 # to interpolate height) as Karney describes, for references 

756 # @see: Function L{ellipsoidalKarney.nearestOn}. 

757 vp, f = _V3d(_0_0, _0_0, h0), None 

758 for i in range(1, _TRIPS): 

759 A.reset(t.lat, t.lon) # gu-/estimate as origin 

760 # convert points to projection space and compute 

761 # the nearest one (and its height) there 

762 s, t = A._forwards(p1, p2) 

763 v, f = _vnOn2(vp, _V3d(s.x, s.y, h1), 

764 _V3d(t.x, t.y, h2), within=within) 

765 # convert nearest one back to geodetic 

766 t, d = A._reverse2(v) 

767 if e.done(d): # below tol or unchanged? 

768 t._iteration = i # _NamedTuple._iteration 

769 break 

770 else: 

771 raise e.degError() 

772 

773 if height is None: 

774 h = v.z # nearest 

775 elif isscalar(height): 

776 h = height 

777 r = _LL4Tuple(t.lat, t.lon, h, t.datum, LatLon, LatLon_kwds, inst=p, 

778 iteration=t.iteration, name=n) 

779 return r, f, e # fraction or None 

780 

781 

782__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI) 

783 

784# **) MIT License 

785# 

786# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

787# 

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

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

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

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

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

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

794# 

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

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

797# 

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

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

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

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

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

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

804# OTHER DEALINGS IN THE SOFTWARE.