Coverage for pygeodesy/ellipsoidalBaseDI.py: 91%

330 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-06-10 14:08 -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_, _ellipsoidal_, \ 

21 _exceed_PI_radians_, _low_, _near_, \ 

22 _SPACE_, _too_ 

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

24 _dunder_nameof 

25from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \ 

26 Intersection3Tuple, NearestOn2Tuple, \ 

27 NearestOn8Tuple, _LL4Tuple 

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

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

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

31 Radius_, Scalar 

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

33 _Wrap, wrap360 

34 

35from math import degrees, radians 

36 

37__all__ = _ALL_LAZY.ellipsoidalBaseDI 

38__version__ = '24.05.19' 

39 

40_polar__ = 'polar?' 

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

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

43 

44 

45class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase): 

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

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

48 ''' 

49 

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

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

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

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

54 for more details. 

55 

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

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

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

59 

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

61 

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

63 

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

65 ellipsoids are not compatible. 

66 ''' 

67 r = self._Inverse(other, wrap) 

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

69 

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

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

72 the given distance from this point along a geodesic given 

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

74 method L{destination2} for more details. 

75 

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

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

78 @kwarg height: Optional height, overriding the default 

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

80 

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

82 ''' 

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

84 

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

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

87 azimuth) after having travelled for the given distance from 

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

89 using this C{Direct} method. 

90 

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

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

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

94 

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

96 are in compass C{degrees360}. 

97 

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

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

100 

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

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

103 @kwarg height: Optional height, overriding the default 

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

105 

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

107 ''' 

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

109 return self._xnamed(r) 

110 

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

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

113 

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

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

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

117 ''' 

118 g = self.geodesic 

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

120 if LL: 

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

122 return r 

123 

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

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

126 ''' 

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

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

129 datum=self.datum, name=self.name, 

130 epoch=self.epoch, reframe=self.reframe)) 

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

132 

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

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

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

136 L{distanceTo3} for more details. 

137 

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

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

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

141 

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

143 

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

145 

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

147 ellipsoids are not compatible. 

148 ''' 

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

150 

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

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

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

154 C{Inverse} method. 

155 

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

157 conventionally meter. The distance is measured on the surface 

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

159 

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

161 are in compass C{degrees360} from North. 

162 

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

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

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

166 

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

168 

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

170 

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

172 ellipsoids are not compatible. 

173 ''' 

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

175 

176 def finalBearingOn(self, distance, bearing): 

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

178 travelled for the given distance along a geodesic given by 

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

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

181 

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

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

184 

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

186 ''' 

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

188 

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

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

191 travelled along a geodesic from this point to an other 

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

193 L{distanceTo3} for more details. 

194 

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

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

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

198 

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

200 

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

202 

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

204 ellipsoids are not compatible. 

205 ''' 

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

207 

208 @property_RO 

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

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

211 ''' 

212 return None # PYCHOK no cover 

213 

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

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

216 ''' 

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

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

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

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

221 return g, gl, p 

222 

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

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

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

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

227 for more details. 

228 

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

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

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

232 

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

234 

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

236 

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

238 ellipsoids are not compatible. 

239 ''' 

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

241 

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

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

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

245 methods. 

246 

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

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

249 at this and 1.0 at the other point. 

250 @kwarg height: Optional height, overriding the fractional 

251 height (C{meter}). 

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

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

254 

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

256 

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

258 

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

260 

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

262 ellipsoids are not compatible. 

263 

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

265 C{rhumbMidpointTo}. 

266 ''' 

267 f = Scalar(fraction=fraction) 

268 if isnear0(f): 

269 r = self 

270 elif isnear1(f) and not wrap: 

271 r = self.others(other) 

272 else: # negative fraction OK 

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

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

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

276 return r 

277 

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

279 wrap=False, tol=_TOL): 

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

281 two points or as a point and bearing. 

282 

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

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

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

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

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

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

289 method L{Ellipsoid.geodesic_}. 

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

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

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

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

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

295 

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

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

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

299 

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

301 

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

303 or B{C{other}} invalid. 

304 

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

306 B{C{height}}. 

307 

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

309 ''' 

310 try: 

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

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

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

314 

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

316 form=_MODS.dms.F_DMS) 

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

318 self.intersecant2) 

319 except (TypeError, ValueError) as x: 

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

321 exact=exact, wrap=wrap) 

322 

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

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

325 

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

327 ''' 

328 _ = self.ellipsoids(other) 

329 g = self.geodesic 

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

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

332 

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

334 equidistant=None, tol=_TOL_M): 

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

336 to this point. 

337 

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

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

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

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

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

343 is taken into account for distances. 

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

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

346 

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

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

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

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

351 

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

353 

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

355 

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

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

358 

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

360 ''' 

361 _d3 = self.distanceTo3 # Distance3Tuple 

362 _n3 = _nearestOn3 

363 try: 

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

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

366 _ = self.ellipsoids(p1) 

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

368 

369 except (TypeError, ValueError) as x: 

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

371 closed=closed, height=height, wrap=wrap) 

372 

373 # get the azimuthal equidistant projection, once 

374 A = _Equidistant00(equidistant, c) 

375 b = _Box(c, c3.distance) 

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

377 

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

379 LatLon=self.classof, # this LatLon 

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

381 try: 

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

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

384 p2 = _unrollon(p1, p2) 

385 # skip edge if no overlap with box around closest 

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

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

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

389 if d3.distance < c3.distance: 

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

391 b = _Box(c, c3.distance) 

392 m = max(m, c.iteration) 

393 p1, i = p2, j 

394 

395 except (TypeError, ValueError) as x: 

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

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

398 closed=closed, height=height, wrap=wrap) 

399 

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

401 

402 n = self.nearestOn8.__name__ # _dunder_nameof 

403 c.rename(n) 

404 if s is not c: 

405 s = s.copy(name=n) 

406 if e is not c: 

407 e = e.copy(name=n) # name__=self.nearestOn8 

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

409 iteration=m) # ._iteration for tests 

410 

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

412 wrap=False, tol=_TOL): 

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

414 a geodesic from this point. 

415 

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

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

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

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

420 see method L{Ellipsoid.geodesic_}. 

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

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

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

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

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

426 

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

428 

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

430 

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

432 ''' 

433 try: 

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

435 

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

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

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

439 p._iteration = P.iteration 

440 except (TypeError, ValueError) as x: 

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

442 exact=exact, wrap=wrap) 

443 return p 

444 

445 

446class _Box(object): 

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

448 

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

450 ''' 

451 def __init__(self, center, distance): 

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

453 

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

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

456 (C{meter}, conventionally) 

457 ''' 

458 m = Radius_(distance=distance) 

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

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

461 self._N = center.lat + d 

462 self._S = center.lat - d 

463 self._E = center.lon + d 

464 self._W = center.lon - d 

465 

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

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

468 

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

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

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

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

473 

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

475 otherwise (C{bool}). 

476 ''' 

477 if lat1 > lat2: 

478 lat1, lat2 = lat2, lat1 

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

480 return False 

481 if lon1 > lon2: 

482 lon1, lon2 = lon2, lon1 

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

484 return False 

485 return True 

486 

487 

488class _Tol(object): 

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

490 ''' 

491 _deg = 0 # tol in degrees 

492 _lat = 0 

493 _m = 0 # tol in meter 

494 _min = MAX # degrees 

495 _prev = None 

496 _r = 0 

497 

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

499 '''New L{_Tol}. 

500 

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

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

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

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

505 ''' 

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

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

508 self._m = max(EPS, tol_m) 

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

510 

511 @property_RO 

512 def degrees(self): 

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

514 ''' 

515 return self._deg 

516 

517 def degrees2m(self, deg): 

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

519 ''' 

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

521 

522 def degError(self, Error=_ValueError): 

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

524 ''' 

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

526 

527 def done(self, deg): 

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

529 ''' 

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

531 return True 

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

533 self._prev = deg 

534 return False 

535 

536 @property_RO 

537 def lat(self): 

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

539 ''' 

540 return self._lat 

541 

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

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

544 ''' 

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

546 if m2km(m) > self.meter: 

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

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

549 

550 @property_RO 

551 def meter(self): 

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

553 ''' 

554 return self._m 

555 

556 @property_RO 

557 def radius(self): 

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

559 ''' 

560 return self._r 

561 

562 def reset(self): 

563 '''Reset tolerances. 

564 ''' 

565 self._min = MAX 

566 self._prev = None 

567 

568 

569def _Equidistant00(equidistant, p1): 

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

571 ''' 

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

573 equidistant = p1.Equidistant 

574 else: 

575 Eqs = _MODS.azimuthal._Equidistants 

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

577 t = map2(_dunder_nameof, Eqs) 

578 raise _IsnotError(*t, equidistant=equidistant) 

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

580 

581 

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

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

584 or as a point and (forward) bearing. 

585 

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

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

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

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

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

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

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

593 

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

595 

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

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

598 

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

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

601 ''' 

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

603 raise _IsnotError(_ellipsoidal_, center=center) 

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

605 

606 

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

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

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

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

611 ''' 

612 _LLS = _MODS.sphericalTrigonometry.LatLon 

613 _si = _MODS.sphericalTrigonometry._intersect 

614 _vi3 = _MODS.vector3d._intersect3d3 

615 

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

617 # compute opposing and distance 

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

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

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

621 return b, t.distance 

622 

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

624 # compute an end point along the initial bearing about 

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

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

627 # radians in .sphericalTrigonometry._int3d2 and bearing 

628 # comparison in .sphericaltrigonometry._intb 

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

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

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

632 e = s.destination(d, e) 

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

634 

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

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

637 ll = not _isDegrees(e) 

638 if ll: 

639 e = s.others(**end) 

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

641 e = _unrollon(s, e) 

642 return e, ll 

643 

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

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

646 if not o: # intersection may be on start 

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

648 return o 

649 return -n if b else n 

650 

651 E = s1.ellipsoids(s2) 

652 

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

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

655 

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

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

658 

659 # get the azimuthal equidistant projection 

660 A = _Equidistant00(equidistant, s1) 

661 

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

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

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

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

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

667 

668 if not ll1: 

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

670 if not ll2: 

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

672 

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

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

675 for i in range(1, _TRIPS): 

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

677 # convert start and end points to projection 

678 # space and compute an intersection there 

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

680 eps=e.meter, useZ=False) 

681 # convert intersection back to geodetic 

682 t, d = A._reverse2(v) 

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

684 break 

685 else: 

686 raise e.degError(Error=IntersectionError) 

687 

688 # like .sphericalTrigonometry._intersect, if this intersection 

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

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

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

692 if b1: 

693 t = t.antipodal() 

694 b1 = not b1 

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

696 

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

698 iteration=i, name=n) 

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

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

701 

702 

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

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

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

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

707 point and an initial bearing from North. 

708 ''' 

709 s1 = _xellipsoidal(start1=start1) 

710 s2 = s1.others(start2=start2) 

711 try: 

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

713 equidistant=equidistant, tol=tol, 

714 LatLon=LatLon, **LatLon_kwds) 

715 except (TypeError, ValueError) as x: 

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

717 

718 

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

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

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

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

723 ''' 

724 c1 = _xellipsoidal(center1=center1) 

725 c2 = c1.others(center2=center2) 

726 try: 

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

728 equidistant=equidistant, tol=tol, 

729 LatLon=LatLon, **LatLon_kwds) 

730 except (TypeError, ValueError) as x: 

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

732 center2=center2, radius2=radius2) 

733 

734 

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

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

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

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

739 ''' 

740 _LLS = _MODS.sphericalTrigonometry.LatLon 

741 _si2 = _MODS.sphericalTrigonometry._intersects2 

742 _vi2 = _MODS.vector3d._intersects2 

743 

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

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

746 iteration=t.iteration, name=n) 

747 

748 r1 = Radius_(radius1=radius1) 

749 r2 = Radius_(radius2=radius2) 

750 

751 E = c1.ellipsoids(c2) 

752 # get the azimuthal equidistant projection 

753 A = _Equidistant00(equidistant, c1) 

754 

755 if r1 < r2: 

756 c1, c2 = c2, c1 

757 r1, r2 = r2, r1 

758 

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

760 raise _ValueError(_exceed_PI_radians_) 

761 

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

763 c2 = _unrollon(c1, c2) 

764 

765 # distance between centers and radii are 

766 # measured along the ellipsoid's surface 

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

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

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

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

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

772 

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

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

775 

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

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

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

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

780 

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

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

783 ts, ta = [], None 

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

785 for i in range(1, _TRIPS): 

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

787 # convert centers to projection space and 

788 # compute the intersections there 

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

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

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

792 sphere=False, too_d=m) 

793 # convert intersections back to geodetic 

794 if v1 is v2: # abutting 

795 t, d = A._reverse2(v1) 

796 else: # consider the closer intersection 

797 t1, d1 = A._reverse2(v1) 

798 t2, d2 = A._reverse2(v2) 

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

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

801 t._iteration = i # _NamedTuple._iteration 

802 ts.append(t) 

803 if v1 is v2: # abutting 

804 ta = t 

805 break 

806 else: 

807 raise e.degError(Error=IntersectionError) 

808 e.reset() 

809 

810 if ta: # abutting circles 

811 pass # PYCHOK no cover 

812 elif len(ts) == 2: 

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

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

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

816 ta = ts[0] # assume abutting 

817 else: # PYCHOK no cover 

818 raise _AssertionError(ts=ts) 

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

820 return r, r 

821 

822 

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

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

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

826 separated to allow callers to embellish any exceptions. 

827 ''' 

828 p1 = p.others(point1=point1) 

829 p2 = p.others(point2=point2) 

830 

831 _ = p.ellipsoids(p1) 

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

833 

834 # get the azimuthal equidistant projection 

835 A = _Equidistant00(equidistant, p) 

836 

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

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

839 tol=tol, **LatLon_and_kwds) 

840 return NearestOn2Tuple(r, f) 

841 

842 

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

844 LatLon=None, **LatLon_kwds): 

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

846 _LLS = _MODS.sphericalNvector.LatLon 

847 _V3d = _MODS.vector3d.Vector3d 

848 _vn2 = _MODS.vector3d._nearestOn2 

849 

850 E = p.ellipsoids(p2) 

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

852 

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

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

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

856 height=height) 

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

858 if height is None: 

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

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

861 h0 = favg(h1, h2) 

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

863 h0 = h1 = h2 = _0_0 

864 

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

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

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

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

869 for i in range(1, _TRIPS): 

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

871 # convert points to projection space and compute 

872 # the nearest one (and its height) there 

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

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

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

876 # convert nearest one back to geodetic 

877 t, d = A._reverse2(v) 

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

879 break 

880 else: 

881 raise e.degError() 

882 

883 if height is None: 

884 h = v.z # nearest 

885 elif _isHeight(height): 

886 h = height 

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

888 iteration=i, name=n) 

889 return r, f, e # fraction or None 

890 

891 

892__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2) 

893 

894# **) MIT License 

895# 

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

897# 

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

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

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

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

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

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

904# 

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

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

907# 

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

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

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

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

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

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

914# OTHER DEALINGS IN THE SOFTWARE.