Coverage for pygeodesy/ellipsoidalBaseDI.py: 91%

330 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-02 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 isLatLon, issubclassof, map2 

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

12 _EPSqrt as _TOL, _0_0, _0_5, _1_5, _3_0 

13# from pygeodesy.dms import F_DMS # _MODS 

14from pygeodesy.ellipsoidalBase import LatLonEllipsoidalBase, _TOL_M, property_RO 

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 _isequalTo, opposing, _radical2 

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

21 _ellipsoidal_, _exceed_PI_radians_, _low_, \ 

22 _near_, _SPACE_, _too_ 

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

24from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \ 

25 Intersection3Tuple, NearestOn2Tuple, \ 

26 NearestOn8Tuple, _LL4Tuple 

27# from pygeodesy.props import property_RO # from .ellipsoidalBase 

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

29from pygeodesy.units import _fi_j2, _isDegrees, _isHeight, _isRadius, \ 

30 Radius_, Scalar 

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

32 _Wrap, wrap360 

33 

34from math import degrees, radians 

35 

36__all__ = _ALL_LAZY.ellipsoidalBaseDI 

37__version__ = '24.02.14' 

38 

39_polar__ = 'polar?' 

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

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

42 

43 

44class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase): 

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

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

47 ''' 

48 

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

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

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

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

53 for more details. 

54 

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

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

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

58 

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

60 

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

62 

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

64 ellipsoids are not compatible. 

65 ''' 

66 r = self._Inverse(other, wrap) 

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

68 

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

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

71 the given distance from this point along a geodesic given 

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

73 method L{destination2} for more details. 

74 

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

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

77 @kwarg height: Optional height, overriding the default 

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

79 

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

81 ''' 

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

83 

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

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

86 azimuth) after having travelled for the given distance from 

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

88 using this C{Direct} method. 

89 

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

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

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

93 

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

95 are in compass C{degrees360}. 

96 

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

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

99 

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

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

102 @kwarg height: Optional height, overriding the default 

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

104 

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

106 ''' 

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

108 return self._xnamed(r) 

109 

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

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

112 

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

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

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

116 ''' 

117 g = self.geodesic 

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

119 if LL: 

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

121 return r 

122 

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

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

125 ''' 

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

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

128 datum=self.datum, name=self.name, 

129 epoch=self.epoch, reframe=self.reframe)) 

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

131 

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

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

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

135 L{distanceTo3} for more details. 

136 

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

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

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

140 

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

142 

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

144 

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

146 ellipsoids are not compatible. 

147 ''' 

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

149 

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

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

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

153 C{Inverse} method. 

154 

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

156 conventionally meter. The distance is measured on the surface 

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

158 

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

160 are in compass C{degrees360} from North. 

161 

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

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

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

165 

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

167 

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

169 

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

171 ellipsoids are not compatible. 

172 ''' 

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

174 

175 def finalBearingOn(self, distance, bearing): 

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

177 travelled for the given distance along a geodesic given by 

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

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

180 

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

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

183 

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

185 ''' 

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

187 

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

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

190 travelled along a geodesic from this point to an other 

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

192 L{distanceTo3} for more details. 

193 

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

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

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

197 

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

199 

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

201 

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

203 ellipsoids are not compatible. 

204 ''' 

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

206 

207 @property_RO 

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

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

210 ''' 

211 return None # PYCHOK no cover 

212 

213 def _g_gl_p3(self, point, other, exact, wrap): 

214 '''(INTERNAL) Helper for methods C{.intersecant2} and C{.plumbTo}. 

215 ''' 

216 p = _unrollon(self, self.others(point=point), wrap=wrap) 

217 g = self.datum.ellipsoid.geodesic_(exact=exact) 

218 gl = g._DirectLine( p, other) if _isDegrees(other) else \ 

219 g._InverseLine(p, self.others(other), wrap) 

220 return g, gl, p 

221 

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

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

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

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

226 for more details. 

227 

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

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

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

231 

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

233 

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

235 

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

237 ellipsoids are not compatible. 

238 ''' 

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

240 

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

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

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

244 methods. 

245 

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

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

248 at this and 1.0 at the other point. 

249 @kwarg height: Optional height, overriding the fractional 

250 height (C{meter}). 

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

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

253 

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

255 

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

257 

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

259 

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

261 ellipsoids are not compatible. 

262 

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

264 C{rhumbMidpointTo}. 

265 ''' 

266 f = Scalar(fraction=fraction) 

267 if isnear0(f): 

268 r = self 

269 elif isnear1(f) and not wrap: 

270 r = self.others(other) 

271 else: # negative fraction OK 

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

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

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

275 return r 

276 

277 def intersecant2(self, circle, point, other, exact=False, height=None, # PYCHOK signature 

278 wrap=False, tol=_TOL): 

279 '''Compute the intersections of a circle and a geodesic (line) given as 

280 two points or as a point and bearing. 

281 

282 @arg circle: Radius of the circle centered at this location (C{meter}, 

283 conventionally) or a point on the circle (this C{LatLon}). 

284 @arg point: A location on the geodesic (this C{LatLon}). 

285 @arg other: An other point on the geodesic (this C{LatLon}) or the 

286 (forward) bearing at the B{C{point}} (compass C{degrees}). 

287 @kwarg exact: Exact C{geodesic...} to use (C{bool} or C{Geodesic...}), see 

288 method L{Ellipsoid.geodesic_}. 

289 @kwarg height: Optional height for the intersection points (C{meter}, 

290 conventionally) or C{None} for interpolated heights. 

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

292 B{C{circle}} and/or B{C{other}} (C{bool}). 

293 @kwarg tol: Convergence tolerance (C{scalar}). 

294 

295 @return: 2-Tuple of the intersection points (representing a geodesic chord), 

296 each an instance of this C{LatLon}. Both points are the same 

297 instance if the geodesic (line) is tangential to the circle. 

298 

299 @raise IntersectionError: The circle and geodesic do not intersect. 

300 

301 @raise TypeError: If B{C{point}} is not this C{LatLon} or B{C{circle}} 

302 or B{C{other}} invalid. 

303 

304 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{exact}} or 

305 B{C{height}}. 

306 

307 @see: Method L{rhumbIntersecant2<LatLonBase.rhumbIntersecant2>}. 

308 ''' 

309 try: 

310 g, gl, p = self._g_gl_p3(point, other, exact, wrap) 

311 r = Radius_(circle=circle) if _isRadius(circle) else \ 

312 g._Inverse(self, self.others(circle=circle), wrap).s12 

313 

314 P, Q = _MODS.geodesicw._Intersecant2(gl, self.lat, self.lon, r, tol=tol, 

315 form=_MODS.dms.F_DMS) 

316 return self._intersecend2(p, other, wrap, height, g, P, Q, 

317 self.intersecant2) 

318 except (TypeError, ValueError) as x: 

319 raise _xError(x, center=self, circle=circle, point=point, other=other, 

320 exact=exact, wrap=wrap) 

321 

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

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

324 

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

326 ''' 

327 _ = self.ellipsoids(other) 

328 g = self.geodesic 

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

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

331 

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

333 equidistant=None, tol=_TOL_M): 

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

335 to this point. 

336 

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

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

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

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

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

342 is taken into account for distances. 

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

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

345 

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

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

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

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

350 

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

352 

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

354 

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

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

357 

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

359 ''' 

360 _d3 = self.distanceTo3 # Distance3Tuple 

361 _n3 = _nearestOn3 

362 try: 

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

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

365 _ = self.ellipsoids(p1) 

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

367 

368 except (TypeError, ValueError) as x: 

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

370 closed=closed, height=height, wrap=wrap) 

371 

372 # get the azimuthal equidistant projection, once 

373 A = _Equidistant00(equidistant, c) 

374 b = _Box(c, c3.distance) 

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

376 

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

378 LatLon=self.classof, # this LatLon 

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

380 try: 

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

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

383 p2 = _unrollon(p1, p2) 

384 # skip edge if no overlap with box around closest 

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

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

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

388 if d3.distance < c3.distance: 

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

390 b = _Box(c, c3.distance) 

391 m = max(m, c.iteration) 

392 p1, i = p2, j 

393 

394 except (TypeError, ValueError) as x: 

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

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

397 closed=closed, height=height, wrap=wrap) 

398 

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

400 

401 n = self.nearestOn8.__name__ 

402 c.rename(n) 

403 if s is not c: 

404 s = s.copy(name=n) 

405 if e is not c: 

406 e = e.copy(name=n) 

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

408 iteration=m) # ._iteration for tests 

409 

410 def plumbTo(self, point, other, exact=False, height=None, # PYCHOK signature 

411 wrap=False, tol=_TOL): 

412 '''Compute the I{perpendicular} intersection of a geodesic (line) with 

413 a geodesic from this point. 

414 

415 @arg point: A location (C{LatLon}) on the geodesic (line). 

416 @arg other: An other point (C{LatLon}) on the geodesic (line) or the 

417 (forward) bearing at the B{C{point}} (compass C{degrees}). 

418 @kwarg exact: Exact C{geodesic...} to use (C{bool} or C{Geodesic...}), 

419 see method L{Ellipsoid.geodesic_}. 

420 @kwarg height: Optional height for the intersection point (C{meter}, 

421 conventionally) or C{None} for an interpolated height. 

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

423 and/or B{C{other}} (C{bool}). 

424 @kwarg tol: Convergence tolerance (C{scalar}). 

425 

426 @return: The intersection point, an instance of this C{LatLon}. 

427 

428 @raise TypeError: If B{C{point}} or B{C{other}} not this C{LatLon}. 

429 

430 @raise UnitError: Invalid B{C{other}}, B{C{exact}} or B{C{height}}. 

431 ''' 

432 try: 

433 g, gl, p = self._g_gl_p3(point, other, exact, wrap) 

434 

435 P = _MODS.geodesicw._PlumbTo(gl, self.lat, self.lon, tol=tol) 

436 h = self._havg(p, h=height) 

437 p = self.classof(P.lat2, P.lon2, datum=self.datum, height=h) # name=n 

438 p._iteration = P.iteration 

439 except (TypeError, ValueError) as x: 

440 raise _xError(x, plumb=self, point=point, other=other, 

441 exact=exact, wrap=wrap) 

442 return p 

443 

444 

445class _Box(object): 

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

447 

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

449 ''' 

450 def __init__(self, center, distance): 

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

452 

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

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

455 (C{meter}, conventionally) 

456 ''' 

457 m = Radius_(distance=distance) 

458 E = center.ellipsoid() # XXX see above 

459 d = E.m2degrees(m) + _0_5 # some margin 

460 self._N = center.lat + d 

461 self._S = center.lat - d 

462 self._E = center.lon + d 

463 self._W = center.lon - d 

464 

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

466 '''Check whether this box overlaps an other box. 

467 

468 @arg lat1: Latitude of a box corner (C{degrees}). 

469 @arg lon1: Longitude of a box corner (C{degrees}). 

470 @arg lat2: Latitude of the opposing corner (C{degrees}). 

471 @arg lon2: Longitude of the opposing corner (C{degrees}). 

472 

473 @return: C{True} if there is some overlap, C{False} 

474 otherwise (C{bool}). 

475 ''' 

476 if lat1 > lat2: 

477 lat1, lat2 = lat2, lat1 

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

479 return False 

480 if lon1 > lon2: 

481 lon1, lon2 = lon2, lon1 

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

483 return False 

484 return True 

485 

486 

487class _Tol(object): 

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

489 ''' 

490 _deg = 0 # tol in degrees 

491 _lat = 0 

492 _m = 0 # tol in meter 

493 _min = MAX # degrees 

494 _prev = None 

495 _r = 0 

496 

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

498 '''New L{_Tol}. 

499 

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

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

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

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

504 ''' 

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

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

507 self._m = max(EPS, tol_m) 

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

509 

510 @property_RO 

511 def degrees(self): 

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

513 ''' 

514 return self._deg 

515 

516 def degrees2m(self, deg): 

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

518 ''' 

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

520 

521 def degError(self, Error=_ValueError): 

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

523 ''' 

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

525 

526 def done(self, deg): 

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

528 ''' 

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

530 return True 

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

532 self._prev = deg 

533 return False 

534 

535 @property_RO 

536 def lat(self): 

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

538 ''' 

539 return self._lat 

540 

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

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

543 ''' 

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

545 if m2km(m) > self.meter: 

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

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

548 

549 @property_RO 

550 def meter(self): 

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

552 ''' 

553 return self._m 

554 

555 @property_RO 

556 def radius(self): 

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

558 ''' 

559 return self._r 

560 

561 def reset(self): 

562 '''Reset tolerances. 

563 ''' 

564 self._min = MAX 

565 self._prev = None 

566 

567 

568def _Equidistant00(equidistant, p1): 

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

570 ''' 

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

572 equidistant = p1.Equidistant 

573 else: 

574 Eqs = _MODS.azimuthal._Equidistants 

575 if not issubclassof(equidistant, *Eqs): # PYCHOK no cover 

576 t = map2(_dunder_nameof, Eqs) 

577 raise _IsnotError(*t, equidistant=equidistant) 

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

579 

580 

581def intersecant2(center, circle, point, other, **exact_height_wrap_tol): 

582 '''Compute the intersections of a circle and a geodesic given as two points 

583 or as a point and (forward) bearing. 

584 

585 @arg center: Center of the circle (C{LatLon}). 

586 @arg circle: Radius of the circle (C{meter}, conventionally) or a point on 

587 the circle (C{LatLon}, as B{C{center}}). 

588 @arg point: A point of the geodesic (C{LatLon}, as B{C{center}}). 

589 @arg other: An other point of the geodesic (C{LatLon}, as B{C{center}}) or 

590 the (forward) bearing at the B{C{point}} (compass C{degrees}). 

591 @kwarg exact_height_wrap_tol: Optional keyword arguments, see below. 

592 

593 @raise NotImplementedError: Method C{intersecant2} not available. 

594 

595 @raise TypeError: If B{C{center}}, B{C{point}} or B{C{circle}} or B{C{other}} 

596 points not ellipsoidal or not compatible with B{C{center}}. 

597 

598 @see: Method C{LatLon.intersecant2} of class L{ellipsoidalExact.LatLon}, 

599 L{ellipsoidalKarney.LatLon} or L{ellipsoidalVincenty.LatLon}. 

600 ''' 

601 if not isLatLon(center, ellipsoidal=True): # isinstance(center, LatLonEllipsoidalBase) 

602 raise _IsnotError(_ellipsoidal_, center=center) 

603 return center.intersecant2(circle, point, other, **exact_height_wrap_tol) 

604 

605 

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

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

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

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

610 ''' 

611 _LLS = _MODS.sphericalTrigonometry.LatLon 

612 _si = _MODS.sphericalTrigonometry._intersect 

613 _vi3 = _MODS.vector3d._intersect3d3 

614 

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

616 # compute opposing and distance 

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

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

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

620 return b, t.distance 

621 

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

623 # compute an end point along the initial bearing about 

624 # 1.5 times the distance to the gu-/estimate, at least 

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

626 # radians in .sphericalTrigonometry._int3d2 and bearing 

627 # comparison in .sphericaltrigonometry._intb 

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

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

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

631 e = s.destination(d, e) 

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

633 

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

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

636 ll = not _isDegrees(e) 

637 if ll: 

638 e = s.others(**end) 

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

640 e = _unrollon(s, e) 

641 return e, ll 

642 

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

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

645 if not o: # intersection may be on start 

646 if _isequalTo(s, t, eps=e.degrees): 

647 return o 

648 return -n if b else n 

649 

650 E = s1.ellipsoids(s2) 

651 

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

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

654 

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

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

657 

658 # get the azimuthal equidistant projection 

659 A = _Equidistant00(equidistant, s1) 

660 

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

662 t = _si(_LLS(s1), (_LLS(e1) if ll1 else e1), 

663 _LLS(s2), (_LLS(e2) if ll2 else e2), 

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

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

666 

667 if not ll1: 

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

669 if not ll2: 

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

671 

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

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

674 for i in range(1, _TRIPS): 

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

676 # convert start and end points to projection 

677 # space and compute an intersection there 

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

679 eps=e.meter, useZ=False) 

680 # convert intersection back to geodetic 

681 t, d = A._reverse2(v) 

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

683 break 

684 else: 

685 raise e.degError(Error=IntersectionError) 

686 

687 # like .sphericalTrigonometry._intersect, if this intersection 

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

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

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

691 if b1: 

692 t = t.antipodal() 

693 b1 = not b1 

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

695 

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

697 iteration=i, name=n) 

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

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

700 

701 

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

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

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

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

706 point and an initial bearing from North. 

707 ''' 

708 s1 = _xellipsoidal(start1=start1) 

709 s2 = s1.others(start2=start2) 

710 try: 

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

712 equidistant=equidistant, tol=tol, 

713 LatLon=LatLon, **LatLon_kwds) 

714 except (TypeError, ValueError) as x: 

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

716 

717 

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

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

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

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

722 ''' 

723 c1 = _xellipsoidal(center1=center1) 

724 c2 = c1.others(center2=center2) 

725 try: 

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

727 equidistant=equidistant, tol=tol, 

728 LatLon=LatLon, **LatLon_kwds) 

729 except (TypeError, ValueError) as x: 

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

731 center2=center2, radius2=radius2) 

732 

733 

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

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

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

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

738 ''' 

739 _LLS = _MODS.sphericalTrigonometry.LatLon 

740 _si2 = _MODS.sphericalTrigonometry._intersects2 

741 _vi2 = _MODS.vector3d._intersects2 

742 

743 def _ll4(t, h, n, c): 

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

745 iteration=t.iteration, name=n) 

746 

747 r1 = Radius_(radius1=radius1) 

748 r2 = Radius_(radius2=radius2) 

749 

750 E = c1.ellipsoids(c2) 

751 # get the azimuthal equidistant projection 

752 A = _Equidistant00(equidistant, c1) 

753 

754 if r1 < r2: 

755 c1, c2 = c2, c1 

756 r1, r2 = r2, r1 

757 

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

759 raise _ValueError(_exceed_PI_radians_) 

760 

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

762 c2 = _unrollon(c1, c2) 

763 

764 # distance between centers and radii are 

765 # measured along the ellipsoid's surface 

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

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

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

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

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

771 

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

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

774 

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

776 t1, t2 = _si2(_LLS(c1), r1, _LLS(c2), r2, radius=e.radius, 

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

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

779 

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

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

782 ts, ta = [], None 

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

784 for i in range(1, _TRIPS): 

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

786 # convert centers to projection space and 

787 # compute the intersections there 

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

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

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

791 sphere=False, too_d=m) 

792 # convert intersections back to geodetic 

793 if v1 is v2: # abutting 

794 t, d = A._reverse2(v1) 

795 else: # consider the closer intersection 

796 t1, d1 = A._reverse2(v1) 

797 t2, d2 = A._reverse2(v2) 

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

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

800 t._iteration = i # _NamedTuple._iteration 

801 ts.append(t) 

802 if v1 is v2: # abutting 

803 ta = t 

804 break 

805 else: 

806 raise e.degError(Error=IntersectionError) 

807 e.reset() 

808 

809 if ta: # abutting circles 

810 pass # PYCHOK no cover 

811 elif len(ts) == 2: 

812 return (_ll4(ts[0], h, n, c1), 

813 _ll4(ts[1], h, n, c2)) 

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

815 ta = ts[0] # assume abutting 

816 else: # PYCHOK no cover 

817 raise _AssertionError(ts=ts) 

818 r = _ll4(ta, h, n, c1) 

819 return r, r 

820 

821 

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

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

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

825 separated to allow callers to embellish any exceptions. 

826 ''' 

827 p1 = p.others(point1=point1) 

828 p2 = p.others(point2=point2) 

829 

830 _ = p.ellipsoids(p1) 

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

832 

833 # get the azimuthal equidistant projection 

834 A = _Equidistant00(equidistant, p) 

835 

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

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

838 tol=tol, **LatLon_and_kwds) 

839 return NearestOn2Tuple(r, f) 

840 

841 

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

843 LatLon=None, **LatLon_kwds): 

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

845 _LLS = _MODS.sphericalNvector.LatLon 

846 _V3d = _MODS.vector3d.Vector3d 

847 _vn2 = _MODS.vector3d._nearestOn2 

848 

849 E = p.ellipsoids(p2) 

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

851 

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

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

854 t = _LLS(p).nearestOn(_LLS(p1), _LLS(p2), within=within, 

855 height=height) 

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

857 if height is None: 

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

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

860 h0 = favg(h1, h2) 

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

862 h0 = h1 = h2 = _0_0 

863 

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

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

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

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

868 for i in range(1, _TRIPS): 

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

870 # convert points to projection space and compute 

871 # the nearest one (and its height) there 

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

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

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

875 # convert nearest one back to geodetic 

876 t, d = A._reverse2(v) 

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

878 break 

879 else: 

880 raise e.degError() 

881 

882 if height is None: 

883 h = v.z # nearest 

884 elif _isHeight(height): 

885 h = height 

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

887 iteration=i, name=n) 

888 return r, f, e # fraction or None 

889 

890 

891__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2) 

892 

893# **) MIT License 

894# 

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

896# 

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

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

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

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

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

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

903# 

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

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

906# 

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

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

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

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

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

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

913# OTHER DEALINGS IN THE SOFTWARE.