Coverage for pygeodesy/ellipsoidalBaseDI.py: 96%

310 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-23 16:38 -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.23' 

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 self._lat = fmean_(lat, *lats) if lats else lat 

407 self._r = max(EPS, E.rocMean(self._lat)) 

408 self._m = max(EPS, tol_m) 

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

410 

411 @property_RO 

412 def degrees(self): 

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

414 ''' 

415 return self._deg 

416 

417 def degrees2m(self, deg): 

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

419 ''' 

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

421 

422 def degError(self, Error=_ValueError): 

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

424 ''' 

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

426 

427 def done(self, deg): 

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

429 ''' 

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

431 return True 

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

433 self._prev = deg 

434 return False 

435 

436 @property_RO 

437 def lat(self): 

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

439 ''' 

440 return self._lat 

441 

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

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

444 ''' 

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

446 if m2km(m) > self.meter: 

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

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

449 

450 @property_RO 

451 def meter(self): 

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

453 ''' 

454 return self._m 

455 

456 @property_RO 

457 def radius(self): 

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

459 ''' 

460 return self._r 

461 

462 def reset(self): 

463 '''Reset tolerances. 

464 ''' 

465 self._min = MAX 

466 self._prev = None 

467 

468 

469def _Equidistant00(equidistant, p1): 

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

471 ''' 

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

473 equidistant = p1.Equidistant 

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

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

476 raise _IsnotError(*t, equidistant=equidistant) 

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

478 

479 

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

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

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

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

484 ''' 

485 _LLS = _MODS.sphericalTrigonometry.LatLon 

486 _si = _MODS.sphericalTrigonometry._intersect 

487 _vi3 = _MODS.vector3d._intersect3d3 

488 

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

490 # compute opposing and distance 

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

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

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

494 return b, t.distance 

495 

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

497 # compute an end point along the initial bearing 

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

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

500 # radians in .sphericalTrigonometry._int3d2 and bearing 

501 # comparison in .sphericalTrigonometry._intb 

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

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

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

505 e = s.destination(d, e) 

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

507 

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

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

510 ll = not isscalar(e) 

511 if ll: 

512 e = s.others(**end) 

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

514 e = _unrollon(s, e) 

515 return e, ll 

516 

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

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

519 

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

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

522 if not o: # intersection may be on start 

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

524 return o 

525 return -n if b else n 

526 

527 E = s1.ellipsoids(s2) 

528 

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

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

531 

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

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

534 

535 # get the azimuthal equidistant projection 

536 A = _Equidistant00(equidistant, s1) 

537 

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

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

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

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

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

543 

544 if not ll1: 

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

546 if not ll2: 

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

548 

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

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

551 for i in range(1, _TRIPS): 

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

553 # convert start and end points to projection 

554 # space and compute an intersection there 

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

556 eps=e.meter, useZ=False) 

557 # convert intersection back to geodetic 

558 t, d = A._reverse2(v) 

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

560 t._iteration = i 

561 break 

562 else: 

563 raise e.degError(Error=IntersectionError) 

564 

565 # like .sphericalTrigonometry._intersect, if this intersection 

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

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

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

569 if b1: 

570 t = t.antipodal() 

571 b1 = not b1 

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

573 

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

575 iteration=t._iteration, name=n) 

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

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

578 

579 

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

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

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

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

584 point and an initial bearing from North. 

585 ''' 

586 s1 = _xellipsoidal(start1=start1) 

587 s2 = s1.others(start2=start2) 

588 try: 

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

590 equidistant=equidistant, tol=tol, 

591 LatLon=LatLon, **LatLon_kwds) 

592 except (TypeError, ValueError) as x: 

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

594 

595 

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

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

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

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

600 ''' 

601 c1 = _xellipsoidal(center1=center1) 

602 c2 = c1.others(center2=center2) 

603 try: 

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

605 equidistant=equidistant, tol=tol, 

606 LatLon=LatLon, **LatLon_kwds) 

607 except (TypeError, ValueError) as x: 

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

609 center2=center2, radius2=radius2) 

610 

611 

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

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

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

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

616 ''' 

617 _LLS = _MODS.sphericalTrigonometry.LatLon 

618 _si2 = _MODS.sphericalTrigonometry._intersects2 

619 _vi2 = _MODS.vector3d._intersects2 

620 

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

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

623 iteration=t.iteration, name=n) 

624 

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

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

627 

628 r1 = Radius_(radius1=radius1) 

629 r2 = Radius_(radius2=radius2) 

630 

631 E = c1.ellipsoids(c2) 

632 # get the azimuthal equidistant projection 

633 A = _Equidistant00(equidistant, c1) 

634 

635 if r1 < r2: 

636 c1, c2 = c2, c1 

637 r1, r2 = r2, r1 

638 

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

640 raise _ValueError(_exceed_PI_radians_) 

641 

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

643 c2 = _unrollon(c1, c2) 

644 

645 # distance between centers and radii are 

646 # measured along the ellipsoid's surface 

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

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

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

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

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

652 

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

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

655 

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

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

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

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

660 

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

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

663 ts, ta = [], None 

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

665 for i in range(1, _TRIPS): 

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

667 # convert centers to projection space and 

668 # compute the intersections there 

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

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

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

672 sphere=False, too_d=m) 

673 # convert intersections back to geodetic 

674 if v1 is v2: # abutting 

675 t, d = A._reverse2(v1) 

676 else: # consider the closer intersection 

677 t1, d1 = A._reverse2(v1) 

678 t2, d2 = A._reverse2(v2) 

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

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

681 t._iteration = i # _NamedTuple._iteration 

682 ts.append(t) 

683 if v1 is v2: # abutting 

684 ta = t 

685 break 

686 else: 

687 raise e.degError(Error=IntersectionError) 

688 e.reset() 

689 

690 if ta: # abutting circles 

691 pass # PYCHOK no cover 

692 elif len(ts) == 2: 

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

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

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

696 ta = ts[0] # assume abutting 

697 else: # PYCHOK no cover 

698 raise _AssertionError(ts=ts) 

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

700 return r, r 

701 

702 

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

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

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

706 separated to allow callers to embellish any exceptions. 

707 ''' 

708 p1 = p.others(point1=point1) 

709 p2 = p.others(point2=point2) 

710 

711 _ = p.ellipsoids(p1) 

712# E = p.ellipsoids(p2) # done in _nearestOn3 

713 

714 # get the azimuthal equidistant projection 

715 A = _Equidistant00(equidistant, p) 

716 

717 if wrap: 

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

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

720 p2 = _unrollon(p1, p2) 

721 

722 r, f, _ = _nearestOn3(p, p1, p2, A, within=within, height=height, 

723 tol=tol, **LatLon_and_kwds) 

724 return NearestOn2Tuple(r, f) 

725 

726 

727def _nearestOn3(p, p1, p2, A, within=True, height=None, tol=_TOL_M, 

728 LatLon=None, **LatLon_kwds): 

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

730 _LLS = _MODS.sphericalNvector.LatLon 

731 _V3d = _MODS.vector3d.Vector3d 

732 _vn2 = _MODS.vector3d._nearestOn2 

733 

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

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

736 

737 E = p.ellipsoids(p2) 

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

739 

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

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

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

743 height=height) 

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

745 if height is None: 

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

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

748 h0 = favg(h1, h2) 

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

750 h0 = h1 = h2 = _0_0 

751 

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

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

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

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

756 for i in range(1, _TRIPS): 

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

758 # convert points to projection space and compute 

759 # the nearest one (and its height) there 

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

761 v, f = _vn2(vp, _V3d(s.x, s.y, h1), 

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

763 # convert nearest one back to geodetic 

764 t, d = A._reverse2(v) 

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

766 t._iteration = i # _NamedTuple._iteration 

767 break 

768 else: 

769 raise e.degError() 

770 

771 if height is None: 

772 h = v.z # nearest 

773 elif isscalar(height): 

774 h = height 

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

776 iteration=t.iteration, name=n) 

777 return r, f, e # fraction or None 

778 

779 

780__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI) 

781 

782# **) MIT License 

783# 

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

785# 

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

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

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

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

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

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

792# 

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

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

795# 

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

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

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

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

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

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

802# OTHER DEALINGS IN THE SOFTWARE.