Coverage for pygeodesy/ellipsoidalBaseDI.py: 91%

337 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-12-24 11:18 -0500

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 

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, property_RO, _TOL_M 

15from pygeodesy.errors import _AssertionError, IntersectionError, _IsnotError, \ 

16 _or, _ValueError, _xellipsoidal, _xError, _xkwds_not 

17from pygeodesy.fmath import favg, fmean_ 

18from pygeodesy.fsums import Fmt, fsumf_ 

19# from pygeodesy.formy import _isequalTo, opposing, _radical2 # _MODS 

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

21 _exceed_PI_radians_, _low_, _near_, _SPACE_, _too_ 

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

23from pygeodesy.namedTuples import Bearing2Tuple, Destination2Tuple, \ 

24 Intersection3Tuple, NearestOn2Tuple, \ 

25 NearestOn8Tuple, _LL4Tuple 

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

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

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

29 Radius_, Scalar 

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

31 _Wrap, wrap360 

32 

33from math import degrees, radians 

34 

35__all__ = _ALL_LAZY.ellipsoidalBaseDI 

36__version__ = '23.12.12' 

37 

38_polar__ = 'polar?' 

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

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

41 

42 

43class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase): 

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

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

46 ''' 

47 

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

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

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

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

52 for more details. 

53 

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

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

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

57 

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

59 

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

61 

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

63 ellipsoids are not compatible. 

64 ''' 

65 r = self._Inverse(other, wrap) 

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

67 

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

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

70 the given distance from this point along a geodesic given 

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

72 method L{destination2} for more details. 

73 

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

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

76 @kwarg height: Optional height, overriding the default 

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

78 

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

80 ''' 

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

82 

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

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

85 azimuth) after having travelled for the given distance from 

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

87 using this C{Direct} method. 

88 

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

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

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

92 

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

94 are in compass C{degrees360}. 

95 

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

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

98 

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

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

101 @kwarg height: Optional height, overriding the default 

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

103 

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

105 ''' 

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

107 return self._xnamed(r) 

108 

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

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

111 

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

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

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

115 ''' 

116 g = self.geodesic 

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

118 if LL: 

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

120 return r 

121 

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

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

124 ''' 

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

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

127 datum=self.datum, name=self.name, 

128 epoch=self.epoch, reframe=self.reframe)) 

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

130 

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

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

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

134 L{distanceTo3} for more details. 

135 

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

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

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

139 

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

141 

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

143 

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

145 ellipsoids are not compatible. 

146 ''' 

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

148 

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

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

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

152 C{Inverse} method. 

153 

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

155 conventionally meter. The distance is measured on the surface 

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

157 

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

159 are in compass C{degrees360} from North. 

160 

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

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

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

164 

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

166 

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

168 

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

170 ellipsoids are not compatible. 

171 ''' 

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

173 

174 def finalBearingOn(self, distance, bearing): 

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

176 travelled for the given distance along a geodesic given by 

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

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

179 

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

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

182 

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

184 ''' 

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

186 

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

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

189 travelled along a geodesic from this point to an other 

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

191 L{distanceTo3} for more details. 

192 

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

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

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

196 

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

198 

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

200 

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

202 ellipsoids are not compatible. 

203 ''' 

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

205 

206 @property_RO 

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

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

209 ''' 

210 return None # PYCHOK no cover 

211 

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

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

214 ''' 

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

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

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

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

219 return g, gl, p 

220 

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

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

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

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

225 for more details. 

226 

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

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

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

230 

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

232 

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

234 

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

236 ellipsoids are not compatible. 

237 ''' 

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

239 

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

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

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

243 methods. 

244 

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

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

247 at this and 1.0 at the other point. 

248 @kwarg height: Optional height, overriding the fractional 

249 height (C{meter}). 

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

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

252 

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

254 

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

256 

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

258 

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

260 ellipsoids are not compatible. 

261 

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

263 C{rhumbMidpointTo}. 

264 ''' 

265 f = Scalar(fraction=fraction) 

266 if isnear0(f): 

267 r = self 

268 elif isnear1(f) and not wrap: 

269 r = self.others(other) 

270 else: # negative fraction OK 

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

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

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

274 return r 

275 

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

277 wrap=False, tol=_TOL): 

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

279 two points or as a point and bearing. 

280 

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

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

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

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

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

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

287 method L{Ellipsoid.geodesic_}. 

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

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

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

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

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

293 

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

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

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

297 

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

299 

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

301 or B{C{other}} invalid. 

302 

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

304 B{C{height}}. 

305 

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

307 ''' 

308 try: 

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

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

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

312 

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

314 form=_MODS.dms.F_DMS) 

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

316 self.intersecant2) 

317 except (TypeError, ValueError) as x: 

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

319 exact=exact, wrap=wrap) 

320 

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

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

323 

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

325 ''' 

326 _ = self.ellipsoids(other) 

327 g = self.geodesic 

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

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

330 

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

332 equidistant=None, tol=_TOL_M): 

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

334 to this point. 

335 

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

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

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

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

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

341 is taken into account for distances. 

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

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

344 

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

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

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

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

349 

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

351 

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

353 

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

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

356 

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

358 ''' 

359 _d3 = self.distanceTo3 # Distance3Tuple 

360 _n3 = _nearestOn3 

361 try: 

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

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

364 _ = self.ellipsoids(p1) 

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

366 

367 except (TypeError, ValueError) as x: 

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

369 closed=closed, height=height, wrap=wrap) 

370 

371 # get the azimuthal equidistant projection, once 

372 A = _Equidistant00(equidistant, c) 

373 b = _Box(c, c3.distance) 

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

375 

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

377 LatLon=self.classof, # this LatLon 

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

379 try: 

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

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

382 p2 = _unrollon(p1, p2) 

383 # skip edge if no overlap with box around closest 

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

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

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

387 if d3.distance < c3.distance: 

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

389 b = _Box(c, c3.distance) 

390 m = max(m, c.iteration) 

391 p1, i = p2, j 

392 

393 except (TypeError, ValueError) as x: 

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

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

396 closed=closed, height=height, wrap=wrap) 

397 

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

399 

400 n = self.nearestOn8.__name__ 

401 c.rename(n) 

402 if s is not c: 

403 s = s.copy(name=n) 

404 if e is not c: 

405 e = e.copy(name=n) 

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

407 iteration=m) # ._iteration for tests 

408 

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

410 wrap=False, tol=_TOL): 

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

412 a geodesic from this point. 

413 

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

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

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

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

418 see method L{Ellipsoid.geodesic_}. 

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

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

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

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

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

424 

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

426 

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

428 

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

430 ''' 

431 try: 

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

433 

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

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

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

437 p._iteration = P.iteration 

438 except (TypeError, ValueError) as x: 

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

440 exact=exact, wrap=wrap) 

441 return p 

442 

443 

444class _Box(object): 

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

446 

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

448 ''' 

449 def __init__(self, center, distance): 

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

451 

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

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

454 (C{meter}, conventionally) 

455 ''' 

456 m = Radius_(distance=distance) 

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

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

459 self._N = center.lat + d 

460 self._S = center.lat - d 

461 self._E = center.lon + d 

462 self._W = center.lon - d 

463 

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

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

466 

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

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

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

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

471 

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

473 otherwise (C{bool}). 

474 ''' 

475 if lat1 > lat2: 

476 lat1, lat2 = lat2, lat1 

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

478 return False 

479 if lon1 > lon2: 

480 lon1, lon2 = lon2, lon1 

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

482 return False 

483 return True 

484 

485 

486class _Tol(object): 

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

488 ''' 

489 _deg = 0 # tol in degrees 

490 _lat = 0 

491 _m = 0 # tol in meter 

492 _min = MAX # degrees 

493 _prev = None 

494 _r = 0 

495 

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

497 '''New L{_Tol}. 

498 

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

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

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

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

503 ''' 

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

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

506 self._m = max(EPS, tol_m) 

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

508 

509 @property_RO 

510 def degrees(self): 

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

512 ''' 

513 return self._deg 

514 

515 def degrees2m(self, deg): 

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

517 ''' 

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

519 

520 def degError(self, Error=_ValueError): 

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

522 ''' 

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

524 

525 def done(self, deg): 

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

527 ''' 

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

529 return True 

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

531 self._prev = deg 

532 return False 

533 

534 @property_RO 

535 def lat(self): 

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

537 ''' 

538 return self._lat 

539 

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

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

542 ''' 

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

544 if m2km(m) > self.meter: 

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

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

547 

548 @property_RO 

549 def meter(self): 

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

551 ''' 

552 return self._m 

553 

554 @property_RO 

555 def radius(self): 

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

557 ''' 

558 return self._r 

559 

560 def reset(self): 

561 '''Reset tolerances. 

562 ''' 

563 self._min = MAX 

564 self._prev = None 

565 

566 

567def _Equidistant00(equidistant, p1): 

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

569 ''' 

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

571 equidistant = p1.Equidistant 

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

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

574 raise _IsnotError(*t, equidistant=equidistant) 

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

576 

577 

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

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

580 or as a point and (forward) bearing. 

581 

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

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

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

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

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

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

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

589 

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

591 

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

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

594 

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

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

597 ''' 

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

599 raise _IsnotError(_ellipsoidal_, center=center) 

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

601 

602 

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

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

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

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

607 ''' 

608 _opp = _MODS.formy.opposing 

609 _LLS = _MODS.sphericalTrigonometry.LatLon 

610 _si = _MODS.sphericalTrigonometry._intersect 

611 _vi3 = _MODS.vector3d._intersect3d3 

612 

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

614 # compute opposing and distance 

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

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

617 b = _opp(e, t.initial) # "before" start 

618 return b, t.distance 

619 

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

621 # compute an end point along the initial bearing 

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

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

624 # radians in .sphericalTrigonometry._int3d2 and bearing 

625 # comparison in .sphericalTrigonometry._intb 

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

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

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

629 e = s.destination(d, e) 

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

631 

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

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

634 ll = not _isDegrees(e) 

635 if ll: 

636 e = s.others(**end) 

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

638 e = _unrollon(s, e) 

639 return e, ll 

640 

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

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

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 _MODS.formy._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 t._iteration = i 

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=t._iteration, 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 _latlon4(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 def _llS(ll): # return an C{_LLS} instance 

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

751 

752 r1 = Radius_(radius1=radius1) 

753 r2 = Radius_(radius2=radius2) 

754 

755 E = c1.ellipsoids(c2) 

756 # get the azimuthal equidistant projection 

757 A = _Equidistant00(equidistant, c1) 

758 

759 if r1 < r2: 

760 c1, c2 = c2, c1 

761 r1, r2 = r2, r1 

762 

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

764 raise _ValueError(_exceed_PI_radians_) 

765 

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

767 c2 = _unrollon(c1, c2) 

768 

769 # distance between centers and radii are 

770 # measured along the ellipsoid's surface 

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

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

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

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

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

776 

777 f = _MODS.formy._radical2(m, r1, r2).ratio # "radical fraction" 

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

779 

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

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

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

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

784 

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

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

787 ts, ta = [], None 

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

789 for i in range(1, _TRIPS): 

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

791 # convert centers to projection space and 

792 # compute the intersections there 

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

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

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

796 sphere=False, too_d=m) 

797 # convert intersections back to geodetic 

798 if v1 is v2: # abutting 

799 t, d = A._reverse2(v1) 

800 else: # consider the closer intersection 

801 t1, d1 = A._reverse2(v1) 

802 t2, d2 = A._reverse2(v2) 

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

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

805 t._iteration = i # _NamedTuple._iteration 

806 ts.append(t) 

807 if v1 is v2: # abutting 

808 ta = t 

809 break 

810 else: 

811 raise e.degError(Error=IntersectionError) 

812 e.reset() 

813 

814 if ta: # abutting circles 

815 pass # PYCHOK no cover 

816 elif len(ts) == 2: 

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

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

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

820 ta = ts[0] # assume abutting 

821 else: # PYCHOK no cover 

822 raise _AssertionError(ts=ts) 

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

824 return r, r 

825 

826 

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

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

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

830 separated to allow callers to embellish any exceptions. 

831 ''' 

832 p1 = p.others(point1=point1) 

833 p2 = p.others(point2=point2) 

834 

835 _ = p.ellipsoids(p1) 

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

837 

838 # get the azimuthal equidistant projection 

839 A = _Equidistant00(equidistant, p) 

840 

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

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

843 tol=tol, **LatLon_and_kwds) 

844 return NearestOn2Tuple(r, f) 

845 

846 

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

848 LatLon=None, **LatLon_kwds): 

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

850 _LLS = _MODS.sphericalNvector.LatLon 

851 _V3d = _MODS.vector3d.Vector3d 

852 _vn2 = _MODS.vector3d._nearestOn2 

853 

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

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

856 

857 E = p.ellipsoids(p2) 

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

859 

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

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

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

863 height=height) 

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

865 if height is None: 

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

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

868 h0 = favg(h1, h2) 

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

870 h0 = h1 = h2 = _0_0 

871 

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

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

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

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

876 for i in range(1, _TRIPS): 

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

878 # convert points to projection space and compute 

879 # the nearest one (and its height) there 

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

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

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

883 # convert nearest one back to geodetic 

884 t, d = A._reverse2(v) 

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

886 t._iteration = i # _NamedTuple._iteration 

887 break 

888 else: 

889 raise e.degError() 

890 

891 if height is None: 

892 h = v.z # nearest 

893 elif _isHeight(height): 

894 h = height 

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

896 iteration=t.iteration, name=n) 

897 return r, f, e # fraction or None 

898 

899 

900__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2) 

901 

902# **) MIT License 

903# 

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

905# 

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

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

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

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

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

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

912# 

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

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

915# 

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

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

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

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

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

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

922# OTHER DEALINGS IN THE SOFTWARE.