Coverage for pygeodesy/ellipsoidalBaseDI.py: 93%

314 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-07 07:28 -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, fsumf_ 

19from pygeodesy.formy import opposing, _radical2 

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

21 _low_, _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, Radius_, Scalar 

29from pygeodesy.utily import m2km, unroll180, _unrollon, _unrollon3, \ 

30 _Wrap, wrap360 

31 

32from math import degrees, radians 

33 

34__all__ = _ALL_LAZY.ellipsoidalBaseDI 

35__version__ = '23.05.15' 

36 

37_polar__ = 'polar?' 

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: If C{True}, wrap or I{normalize} and unroll the 

55 B{C{other}} point (C{bool}). 

56 

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

58 

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

60 

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

62 ellipsoids are not compatible. 

63 ''' 

64 r = self._Inverse(other, wrap) 

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

66 

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

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

69 the given distance from this point along a geodesic given 

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

71 method L{destination2} for more details. 

72 

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

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

75 @kwarg height: Optional height, overriding the default 

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

77 

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

79 ''' 

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

81 

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

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

84 azimuth) after having travelled for the given distance from 

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

86 using this C{Direct} method. 

87 

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

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

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

91 

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

93 are in compass C{degrees360}. 

94 

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

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

97 

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

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

100 @kwarg height: Optional height, overriding the default 

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

102 

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

104 ''' 

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

106 return self._xnamed(r) 

107 

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

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

110 

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

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

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

114 ''' 

115 g = self.geodesic 

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

117 if LL: 

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

119 return r 

120 

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

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

123 ''' 

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

125 d = LL(*_Wrap.latlon(r.lat, r.lon), height=h, **_xkwds_not(None, 

126 datum=self.datum, name=self.name, 

127 epoch=self.epoch, reframe=self.reframe)) 

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

129 

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

131 '''Compute the distance between this and an other point along 

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

133 L{distanceTo3} for more details. 

134 

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

136 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the 

137 B{C{other}} point (C{bool}). 

138 

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

140 

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

142 

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

144 ellipsoids are not compatible. 

145 ''' 

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

147 

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

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

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

151 C{Inverse} method. 

152 

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

154 conventionally meter. The distance is measured on the surface 

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

156 

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

158 are in compass C{degrees360} from North. 

159 

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

161 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the 

162 B{C{other}} point (C{bool}). 

163 

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

165 

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

167 

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

169 ellipsoids are not compatible. 

170 ''' 

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

172 

173 def finalBearingOn(self, distance, bearing): 

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

175 travelled for the given distance along a geodesic given by 

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

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

178 

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

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

181 

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

183 ''' 

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

185 

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

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

188 travelled along a geodesic from this point to an other 

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

190 L{distanceTo3} for more details. 

191 

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

193 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll 

194 the B{C{other}} point (C{bool}). 

195 

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

197 

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

199 

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

201 ellipsoids are not compatible. 

202 ''' 

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

204 

205 @Property_RO 

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

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

208 ''' 

209 return None # PYCHOK no cover 

210 

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

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

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

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

215 for more details. 

216 

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

218 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll 

219 the B{C{other}} point (C{bool}). 

220 

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

222 

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

224 

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

226 ellipsoids are not compatible. 

227 ''' 

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

229 

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

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

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

233 methods. 

234 

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

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

237 at this and 1.0 at the other point. 

238 @kwarg height: Optional height, overriding the fractional 

239 height (C{meter}). 

240 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the 

241 B{C{other}} point (C{bool}). 

242 

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

244 

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

246 

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

248 

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

250 ellipsoids are not compatible. 

251 

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

253 C{rhumbMidpointTo}. 

254 ''' 

255 f = Scalar(fraction=fraction) 

256 if isnear0(f): 

257 r = self 

258 elif isnear1(f) and not wrap: 

259 r = self.others(other) 

260 else: # negative fraction OK 

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

262 h = self._havg(other, f=f, h=height) 

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

264 return r 

265 

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

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

268 

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

270 ''' 

271 _ = self.ellipsoids(other) 

272 g = self.geodesic 

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

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

275 

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

277 equidistant=None, tol=_TOL_M): 

278 '''I{Iteratively} locate the point on a path or polygon closest 

279 to this point. 

280 

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

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

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

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

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

286 is taken into account for distances. 

287 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the 

288 B{C{points}} (C{bool}). 

289 

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

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

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

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

294 

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

296 

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

298 

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

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

301 

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

303 ''' 

304 _d3 = self.distanceTo3 # Distance3Tuple 

305 _n3 = _nearestOn3 

306 try: 

307 Ps = self.PointsIter(points, loop=1, wrap=wrap) 

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

309 _ = self.ellipsoids(p1) 

310 c3 = _d3(c, wrap=wrap) # XXX wrap=False? 

311 

312 except (TypeError, ValueError) as x: 

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

314 closed=closed, height=height, wrap=wrap) 

315 

316 # get the azimuthal equidistant projection, once 

317 A = _Equidistant00(equidistant, c) 

318 b = _Box(c, c3.distance) 

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

320 

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

322 LatLon=self.classof, # this LatLon 

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

324 try: 

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

326 if wrap and j != 0: # not Ps.looped 

327 p2 = _unrollon(p1, p2) 

328 # skip edge if no overlap with box around closest 

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

330 p, t, _ = _n3(self, p1, p2, A, **kwds) 

331 d3 = _d3(p, wrap=False) # already unrolled 

332 if d3.distance < c3.distance: 

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

334 b = _Box(c, c3.distance) 

335 m = max(m, c.iteration) 

336 p1, i = p2, j 

337 

338 except (TypeError, ValueError) as x: 

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

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

341 closed=closed, height=height, wrap=wrap) 

342 

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

344 

345 n = self.nearestOn8.__name__ 

346 c.rename(n) 

347 if s is not c: 

348 s = s.copy(name=n) 

349 if e is not c: 

350 e = e.copy(name=n) 

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

352 iteration=m) # ._iteration for tests 

353 

354 

355class _Box(object): 

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

357 

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

359 ''' 

360 def __init__(self, center, distance): 

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

362 

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

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

365 (C{meter}, conventionally) 

366 ''' 

367 E = center.ellipsoid() 

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

369 self._N = center.lat + d 

370 self._S = center.lat - d 

371 self._E = center.lon + d 

372 self._W = center.lon - d 

373 

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

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

376 

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

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

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

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

381 

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

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

384 ''' 

385 if lat1 > lat2: 

386 lat1, lat2 = lat2, lat1 

387 if lat1 > self._N or lat2 < self._S: 

388 return False 

389 if lon1 > lon2: 

390 lon1, lon2 = lon2, lon1 

391 if lon1 > self._E or lon2 < self._W: 

392 return False 

393 return True 

394 

395 

396class _Tol(object): 

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

398 ''' 

399 _deg = 0 # tol in degrees 

400 _lat = 0 

401 _m = 0 # tol in meter 

402 _min = MAX # degrees 

403 _prev = None 

404 _r = 0 

405 

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

407 '''New L{_Tol}. 

408 

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

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

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

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

413 ''' 

414 self._lat = fmean_(lat, *lats) if lats else lat 

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

416 self._m = max(EPS, tol_m) 

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

418 

419 @property_RO 

420 def degrees(self): 

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

422 ''' 

423 return self._deg 

424 

425 def degrees2m(self, deg): 

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

427 ''' 

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

429 

430 def degError(self, Error=_ValueError): 

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

432 ''' 

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

434 

435 def done(self, deg): 

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

437 ''' 

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

439 return True 

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

441 self._prev = deg 

442 return False 

443 

444 @property_RO 

445 def lat(self): 

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

447 ''' 

448 return self._lat 

449 

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

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

452 ''' 

453 t = _SPACE_(Fmt.tolerance(self.meter), _too_(_low_)) 

454 if m2km(m) > self.meter: 

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

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

457 

458 @property_RO 

459 def meter(self): 

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

461 ''' 

462 return self._m 

463 

464 @property_RO 

465 def radius(self): 

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

467 ''' 

468 return self._r 

469 

470 def reset(self): 

471 '''Reset tolerances. 

472 ''' 

473 self._min = MAX 

474 self._prev = None 

475 

476 

477def _Equidistant00(equidistant, p1): 

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

479 ''' 

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

481 equidistant = p1.Equidistant 

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

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

484 raise _IsnotError(*t, equidistant=equidistant) 

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

486 

487 

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

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

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

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

492 ''' 

493 _LLS = _MODS.sphericalTrigonometry.LatLon 

494 _si = _MODS.sphericalTrigonometry._intersect 

495 _vi3 = _MODS.vector3d._intersect3d3 

496 

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

498 # compute opposing and distance 

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

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

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

502 return b, t.distance 

503 

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

505 # compute an end point along the initial bearing 

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

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

508 # radians in .sphericalTrigonometry._int3d2 and bearing 

509 # comparison in .sphericalTrigonometry._intb 

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

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

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

513 e = s.destination(d, e) 

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

515 

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

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

518 ll = not isscalar(e) 

519 if ll: 

520 e = s.others(**end) 

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

522 e = _unrollon(s, e) 

523 return e, ll 

524 

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

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

527 

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

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

530 if not o: # intersection may be on start 

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

532 return o 

533 return -n if b else n 

534 

535 E = s1.ellipsoids(s2) 

536 

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

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

539 

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

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

542 

543 # get the azimuthal equidistant projection 

544 A = _Equidistant00(equidistant, s1) 

545 

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

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

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

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

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

551 

552 if not ll1: 

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

554 if not ll2: 

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

556 

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

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

559 for i in range(1, _TRIPS): 

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

561 # convert start and end points to projection 

562 # space and compute an intersection there 

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

564 eps=e.meter, useZ=False) 

565 # convert intersection back to geodetic 

566 t, d = A._reverse2(v) 

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

568 t._iteration = i 

569 break 

570 else: 

571 raise e.degError(Error=IntersectionError) 

572 

573 # like .sphericalTrigonometry._intersect, if this intersection 

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

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

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

577 if b1: 

578 t = t.antipodal() 

579 b1 = not b1 

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

581 

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

583 iteration=t._iteration, name=n) 

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

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

586 

587 

588def _intersection3(start1, end1, start2, end2, height=None, wrap=False, # was=True 

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

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

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

592 point and an initial bearing from North. 

593 ''' 

594 s1 = _xellipsoidal(start1=start1) 

595 s2 = s1.others(start2=start2) 

596 try: 

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

598 equidistant=equidistant, tol=tol, 

599 LatLon=LatLon, **LatLon_kwds) 

600 except (TypeError, ValueError) as x: 

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

602 

603 

604def _intersections2(center1, radius1, center2, radius2, height=None, wrap=False, # was=True 

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

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

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

608 ''' 

609 c1 = _xellipsoidal(center1=center1) 

610 c2 = c1.others(center2=center2) 

611 try: 

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

613 equidistant=equidistant, tol=tol, 

614 LatLon=LatLon, **LatLon_kwds) 

615 except (TypeError, ValueError) as x: 

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

617 center2=center2, radius2=radius2) 

618 

619 

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

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

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

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

624 ''' 

625 _LLS = _MODS.sphericalTrigonometry.LatLon 

626 _si2 = _MODS.sphericalTrigonometry._intersects2 

627 _vi2 = _MODS.vector3d._intersects2 

628 

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

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

631 iteration=t.iteration, name=n) 

632 

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

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

635 

636 r1 = Radius_(radius1=radius1) 

637 r2 = Radius_(radius2=radius2) 

638 

639 E = c1.ellipsoids(c2) 

640 # get the azimuthal equidistant projection 

641 A = _Equidistant00(equidistant, c1) 

642 

643 if r1 < r2: 

644 c1, c2 = c2, c1 

645 r1, r2 = r2, r1 

646 

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

648 raise _ValueError(_exceed_PI_radians_) 

649 

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

651 c2 = _unrollon(c1, c2) 

652 

653 # distance between centers and radii are 

654 # measured along the ellipsoid's surface 

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

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

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

658 if fsumf_(r1, r2, -m) < 0: 

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

660 

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

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

663 

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

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

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

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

668 

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

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

671 ts, ta = [], None 

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

673 for i in range(1, _TRIPS): 

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

675 # convert centers to projection space and 

676 # compute the intersections there 

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

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

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

680 sphere=False, too_d=m) 

681 # convert intersections back to geodetic 

682 if v1 is v2: # abutting 

683 t, d = A._reverse2(v1) 

684 else: # consider the closer intersection 

685 t1, d1 = A._reverse2(v1) 

686 t2, d2 = A._reverse2(v2) 

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

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

689 t._iteration = i # _NamedTuple._iteration 

690 ts.append(t) 

691 if v1 is v2: # abutting 

692 ta = t 

693 break 

694 else: 

695 raise e.degError(Error=IntersectionError) 

696 e.reset() 

697 

698 if ta: # abutting circles 

699 pass # PYCHOK no cover 

700 elif len(ts) == 2: 

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

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

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

704 ta = ts[0] # assume abutting 

705 else: # PYCHOK no cover 

706 raise _AssertionError(ts=ts) 

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

708 return r, r 

709 

710 

711def _nearestOn2(p, point1, point2, within=True, height=None, wrap=False, # was=True 

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

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

714 separated to allow callers to embellish any exceptions. 

715 ''' 

716 p1 = p.others(point1=point1) 

717 p2 = p.others(point2=point2) 

718 

719 _ = p.ellipsoids(p1) 

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

721 

722 # get the azimuthal equidistant projection 

723 A = _Equidistant00(equidistant, p) 

724 

725 p1, p2, _ = _unrollon3(p, p1, p2, wrap) # XXX don't unroll? 

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

727 tol=tol, **LatLon_and_kwds) 

728 return NearestOn2Tuple(r, f) 

729 

730 

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

732 LatLon=None, **LatLon_kwds): 

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

734 _LLS = _MODS.sphericalNvector.LatLon 

735 _V3d = _MODS.vector3d.Vector3d 

736 _vn2 = _MODS.vector3d._nearestOn2 

737 

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

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

740 

741 E = p.ellipsoids(p2) 

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

743 

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

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

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

747 height=height) 

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

749 if height is None: 

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

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

752 h0 = favg(h1, h2) 

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

754 h0 = h1 = h2 = _0_0 

755 

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

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

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

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

760 for i in range(1, _TRIPS): 

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

762 # convert points to projection space and compute 

763 # the nearest one (and its height) there 

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

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

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

767 # convert nearest one back to geodetic 

768 t, d = A._reverse2(v) 

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

770 t._iteration = i # _NamedTuple._iteration 

771 break 

772 else: 

773 raise e.degError() 

774 

775 if height is None: 

776 h = v.z # nearest 

777 elif isscalar(height): 

778 h = height 

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

780 iteration=t.iteration, name=n) 

781 return r, f, e # fraction or None 

782 

783 

784__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI) 

785 

786# **) MIT License 

787# 

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

789# 

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

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

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

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

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

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

796# 

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

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

799# 

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

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

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

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

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

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

806# OTHER DEALINGS IN THE SOFTWARE.