Coverage for pygeodesy/ellipsoidalBaseDI.py: 91%

330 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-15 16:36 -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 

20# from pygeodesy.internals import _dunder_nameof # from .lazily 

21from pygeodesy.interns import _antipodal_, _concentric_, _ellipsoidal_, \ 

22 _exceed_PI_radians_, _low_, _near_, \ 

23 _SPACE_, _too_ 

24from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, \ 

25 _dunder_nameof 

26from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \ 

27 Intersection3Tuple, NearestOn2Tuple, \ 

28 NearestOn8Tuple, _LL4Tuple 

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

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

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

32 Radius_, Scalar 

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

34 _Wrap, wrap360 

35 

36from math import degrees, radians 

37 

38__all__ = _ALL_LAZY.ellipsoidalBaseDI 

39__version__ = '24.05.13' 

40 

41_polar__ = 'polar?' 

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

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

44 

45 

46class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase): 

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

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

49 ''' 

50 

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

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

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

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

55 for more details. 

56 

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

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

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

60 

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

62 

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

64 

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

66 ellipsoids are not compatible. 

67 ''' 

68 r = self._Inverse(other, wrap) 

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

70 

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

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

73 the given distance from this point along a geodesic given 

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

75 method L{destination2} for more details. 

76 

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

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

79 @kwarg height: Optional height, overriding the default 

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

81 

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

83 ''' 

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

85 

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

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

88 azimuth) after having travelled for the given distance from 

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

90 using this C{Direct} method. 

91 

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

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

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

95 

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

97 are in compass C{degrees360}. 

98 

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

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

101 

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

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

104 @kwarg height: Optional height, overriding the default 

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

106 

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

108 ''' 

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

110 return self._xnamed(r) 

111 

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

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

114 

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

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

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

118 ''' 

119 g = self.geodesic 

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

121 if LL: 

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

123 return r 

124 

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

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

127 ''' 

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

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

130 datum=self.datum, name=self.name, 

131 epoch=self.epoch, reframe=self.reframe)) 

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

133 

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

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

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

137 L{distanceTo3} for more details. 

138 

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

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

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

142 

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

144 

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

146 

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

148 ellipsoids are not compatible. 

149 ''' 

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

151 

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

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

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

155 C{Inverse} method. 

156 

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

158 conventionally meter. The distance is measured on the surface 

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

160 

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

162 are in compass C{degrees360} from North. 

163 

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

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

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

167 

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

169 

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

171 

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

173 ellipsoids are not compatible. 

174 ''' 

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

176 

177 def finalBearingOn(self, distance, bearing): 

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

179 travelled for the given distance along a geodesic given by 

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

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

182 

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

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

185 

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

187 ''' 

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

189 

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

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

192 travelled along a geodesic from this point to an other 

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

194 L{distanceTo3} for more details. 

195 

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

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

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

199 

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

201 

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

203 

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

205 ellipsoids are not compatible. 

206 ''' 

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

208 

209 @property_RO 

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

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

212 ''' 

213 return None # PYCHOK no cover 

214 

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

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

217 ''' 

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

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

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

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

222 return g, gl, p 

223 

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

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

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

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

228 for more details. 

229 

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

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

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

233 

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

235 

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

237 

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

239 ellipsoids are not compatible. 

240 ''' 

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

242 

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

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

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

246 methods. 

247 

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

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

250 at this and 1.0 at the other point. 

251 @kwarg height: Optional height, overriding the fractional 

252 height (C{meter}). 

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

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

255 

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

257 

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

259 

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

261 

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

263 ellipsoids are not compatible. 

264 

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

266 C{rhumbMidpointTo}. 

267 ''' 

268 f = Scalar(fraction=fraction) 

269 if isnear0(f): 

270 r = self 

271 elif isnear1(f) and not wrap: 

272 r = self.others(other) 

273 else: # negative fraction OK 

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

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

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

277 return r 

278 

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

280 wrap=False, tol=_TOL): 

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

282 two points or as a point and bearing. 

283 

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

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

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

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

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

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

290 method L{Ellipsoid.geodesic_}. 

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

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

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

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

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

296 

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

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

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

300 

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

302 

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

304 or B{C{other}} invalid. 

305 

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

307 B{C{height}}. 

308 

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

310 ''' 

311 try: 

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

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

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

315 

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

317 form=_MODS.dms.F_DMS) 

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

319 self.intersecant2) 

320 except (TypeError, ValueError) as x: 

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

322 exact=exact, wrap=wrap) 

323 

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

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

326 

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

328 ''' 

329 _ = self.ellipsoids(other) 

330 g = self.geodesic 

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

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

333 

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

335 equidistant=None, tol=_TOL_M): 

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

337 to this point. 

338 

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

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

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

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

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

344 is taken into account for distances. 

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

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

347 

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

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

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

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

352 

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

354 

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

356 

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

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

359 

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

361 ''' 

362 _d3 = self.distanceTo3 # Distance3Tuple 

363 _n3 = _nearestOn3 

364 try: 

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

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

367 _ = self.ellipsoids(p1) 

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

369 

370 except (TypeError, ValueError) as x: 

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

372 closed=closed, height=height, wrap=wrap) 

373 

374 # get the azimuthal equidistant projection, once 

375 A = _Equidistant00(equidistant, c) 

376 b = _Box(c, c3.distance) 

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

378 

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

380 LatLon=self.classof, # this LatLon 

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

382 try: 

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

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

385 p2 = _unrollon(p1, p2) 

386 # skip edge if no overlap with box around closest 

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

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

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

390 if d3.distance < c3.distance: 

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

392 b = _Box(c, c3.distance) 

393 m = max(m, c.iteration) 

394 p1, i = p2, j 

395 

396 except (TypeError, ValueError) as x: 

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

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

399 closed=closed, height=height, wrap=wrap) 

400 

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

402 

403 n = self.nearestOn8.__name__ 

404 c.rename(n) 

405 if s is not c: 

406 s = s.copy(name=n) 

407 if e is not c: 

408 e = e.copy(name=n) 

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

410 iteration=m) # ._iteration for tests 

411 

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

413 wrap=False, tol=_TOL): 

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

415 a geodesic from this point. 

416 

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

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

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

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

421 see method L{Ellipsoid.geodesic_}. 

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

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

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

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

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

427 

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

429 

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

431 

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

433 ''' 

434 try: 

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

436 

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

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

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

440 p._iteration = P.iteration 

441 except (TypeError, ValueError) as x: 

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

443 exact=exact, wrap=wrap) 

444 return p 

445 

446 

447class _Box(object): 

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

449 

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

451 ''' 

452 def __init__(self, center, distance): 

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

454 

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

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

457 (C{meter}, conventionally) 

458 ''' 

459 m = Radius_(distance=distance) 

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

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

462 self._N = center.lat + d 

463 self._S = center.lat - d 

464 self._E = center.lon + d 

465 self._W = center.lon - d 

466 

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

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

469 

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

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

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

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

474 

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

476 otherwise (C{bool}). 

477 ''' 

478 if lat1 > lat2: 

479 lat1, lat2 = lat2, lat1 

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

481 return False 

482 if lon1 > lon2: 

483 lon1, lon2 = lon2, lon1 

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

485 return False 

486 return True 

487 

488 

489class _Tol(object): 

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

491 ''' 

492 _deg = 0 # tol in degrees 

493 _lat = 0 

494 _m = 0 # tol in meter 

495 _min = MAX # degrees 

496 _prev = None 

497 _r = 0 

498 

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

500 '''New L{_Tol}. 

501 

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

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

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

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

506 ''' 

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

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

509 self._m = max(EPS, tol_m) 

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

511 

512 @property_RO 

513 def degrees(self): 

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

515 ''' 

516 return self._deg 

517 

518 def degrees2m(self, deg): 

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

520 ''' 

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

522 

523 def degError(self, Error=_ValueError): 

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

525 ''' 

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

527 

528 def done(self, deg): 

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

530 ''' 

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

532 return True 

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

534 self._prev = deg 

535 return False 

536 

537 @property_RO 

538 def lat(self): 

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

540 ''' 

541 return self._lat 

542 

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

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

545 ''' 

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

547 if m2km(m) > self.meter: 

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

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

550 

551 @property_RO 

552 def meter(self): 

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

554 ''' 

555 return self._m 

556 

557 @property_RO 

558 def radius(self): 

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

560 ''' 

561 return self._r 

562 

563 def reset(self): 

564 '''Reset tolerances. 

565 ''' 

566 self._min = MAX 

567 self._prev = None 

568 

569 

570def _Equidistant00(equidistant, p1): 

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

572 ''' 

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

574 equidistant = p1.Equidistant 

575 else: 

576 Eqs = _MODS.azimuthal._Equidistants 

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

578 t = map2(_dunder_nameof, Eqs) 

579 raise _IsnotError(*t, equidistant=equidistant) 

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

581 

582 

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

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

585 or as a point and (forward) bearing. 

586 

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

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

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

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

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

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

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

594 

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

596 

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

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

599 

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

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

602 ''' 

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

604 raise _IsnotError(_ellipsoidal_, center=center) 

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

606 

607 

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

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

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

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

612 ''' 

613 _LLS = _MODS.sphericalTrigonometry.LatLon 

614 _si = _MODS.sphericalTrigonometry._intersect 

615 _vi3 = _MODS.vector3d._intersect3d3 

616 

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

618 # compute opposing and distance 

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

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

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

622 return b, t.distance 

623 

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

625 # compute an end point along the initial bearing about 

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

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

628 # radians in .sphericalTrigonometry._int3d2 and bearing 

629 # comparison in .sphericaltrigonometry._intb 

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

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

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

633 e = s.destination(d, e) 

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

635 

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

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

638 ll = not _isDegrees(e) 

639 if ll: 

640 e = s.others(**end) 

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

642 e = _unrollon(s, e) 

643 return e, ll 

644 

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

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

647 if not o: # intersection may be on start 

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

649 return o 

650 return -n if b else n 

651 

652 E = s1.ellipsoids(s2) 

653 

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

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

656 

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

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

659 

660 # get the azimuthal equidistant projection 

661 A = _Equidistant00(equidistant, s1) 

662 

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

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

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

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

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

668 

669 if not ll1: 

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

671 if not ll2: 

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

673 

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

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

676 for i in range(1, _TRIPS): 

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

678 # convert start and end points to projection 

679 # space and compute an intersection there 

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

681 eps=e.meter, useZ=False) 

682 # convert intersection back to geodetic 

683 t, d = A._reverse2(v) 

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

685 break 

686 else: 

687 raise e.degError(Error=IntersectionError) 

688 

689 # like .sphericalTrigonometry._intersect, if this intersection 

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

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

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

693 if b1: 

694 t = t.antipodal() 

695 b1 = not b1 

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

697 

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

699 iteration=i, name=n) 

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

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

702 

703 

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

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

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

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

708 point and an initial bearing from North. 

709 ''' 

710 s1 = _xellipsoidal(start1=start1) 

711 s2 = s1.others(start2=start2) 

712 try: 

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

714 equidistant=equidistant, tol=tol, 

715 LatLon=LatLon, **LatLon_kwds) 

716 except (TypeError, ValueError) as x: 

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

718 

719 

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

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

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

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

724 ''' 

725 c1 = _xellipsoidal(center1=center1) 

726 c2 = c1.others(center2=center2) 

727 try: 

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

729 equidistant=equidistant, tol=tol, 

730 LatLon=LatLon, **LatLon_kwds) 

731 except (TypeError, ValueError) as x: 

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

733 center2=center2, radius2=radius2) 

734 

735 

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

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

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

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

740 ''' 

741 _LLS = _MODS.sphericalTrigonometry.LatLon 

742 _si2 = _MODS.sphericalTrigonometry._intersects2 

743 _vi2 = _MODS.vector3d._intersects2 

744 

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

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

747 iteration=t.iteration, name=n) 

748 

749 r1 = Radius_(radius1=radius1) 

750 r2 = Radius_(radius2=radius2) 

751 

752 E = c1.ellipsoids(c2) 

753 # get the azimuthal equidistant projection 

754 A = _Equidistant00(equidistant, c1) 

755 

756 if r1 < r2: 

757 c1, c2 = c2, c1 

758 r1, r2 = r2, r1 

759 

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

761 raise _ValueError(_exceed_PI_radians_) 

762 

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

764 c2 = _unrollon(c1, c2) 

765 

766 # distance between centers and radii are 

767 # measured along the ellipsoid's surface 

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

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

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

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

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

773 

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

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

776 

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

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

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

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

781 

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

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

784 ts, ta = [], None 

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

786 for i in range(1, _TRIPS): 

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

788 # convert centers to projection space and 

789 # compute the intersections there 

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

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

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

793 sphere=False, too_d=m) 

794 # convert intersections back to geodetic 

795 if v1 is v2: # abutting 

796 t, d = A._reverse2(v1) 

797 else: # consider the closer intersection 

798 t1, d1 = A._reverse2(v1) 

799 t2, d2 = A._reverse2(v2) 

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

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

802 t._iteration = i # _NamedTuple._iteration 

803 ts.append(t) 

804 if v1 is v2: # abutting 

805 ta = t 

806 break 

807 else: 

808 raise e.degError(Error=IntersectionError) 

809 e.reset() 

810 

811 if ta: # abutting circles 

812 pass # PYCHOK no cover 

813 elif len(ts) == 2: 

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

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

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

817 ta = ts[0] # assume abutting 

818 else: # PYCHOK no cover 

819 raise _AssertionError(ts=ts) 

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

821 return r, r 

822 

823 

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

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

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

827 separated to allow callers to embellish any exceptions. 

828 ''' 

829 p1 = p.others(point1=point1) 

830 p2 = p.others(point2=point2) 

831 

832 _ = p.ellipsoids(p1) 

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

834 

835 # get the azimuthal equidistant projection 

836 A = _Equidistant00(equidistant, p) 

837 

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

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

840 tol=tol, **LatLon_and_kwds) 

841 return NearestOn2Tuple(r, f) 

842 

843 

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

845 LatLon=None, **LatLon_kwds): 

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

847 _LLS = _MODS.sphericalNvector.LatLon 

848 _V3d = _MODS.vector3d.Vector3d 

849 _vn2 = _MODS.vector3d._nearestOn2 

850 

851 E = p.ellipsoids(p2) 

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

853 

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

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

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

857 height=height) 

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

859 if height is None: 

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

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

862 h0 = favg(h1, h2) 

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

864 h0 = h1 = h2 = _0_0 

865 

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

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

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

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

870 for i in range(1, _TRIPS): 

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

872 # convert points to projection space and compute 

873 # the nearest one (and its height) there 

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

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

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

877 # convert nearest one back to geodetic 

878 t, d = A._reverse2(v) 

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

880 break 

881 else: 

882 raise e.degError() 

883 

884 if height is None: 

885 h = v.z # nearest 

886 elif _isHeight(height): 

887 h = height 

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

889 iteration=i, name=n) 

890 return r, f, e # fraction or None 

891 

892 

893__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2) 

894 

895# **) MIT License 

896# 

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

898# 

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

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

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

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

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

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

905# 

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

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

908# 

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

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

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

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

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

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

915# OTHER DEALINGS IN THE SOFTWARE.