Coverage for pygeodesy/ellipsoidalBaseDI.py: 91%

337 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-12-02 13:46 -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, isscalar, 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_ 

19from pygeodesy.formy import opposing, _radical2 

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, Radius_, Scalar 

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

30 _Wrap, wrap360 

31 

32from math import degrees, radians 

33 

34__all__ = _ALL_LAZY.ellipsoidalBaseDI 

35__version__ = '23.11.16' 

36 

37_polar__ = 'polar?' 

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

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

40 

41 

42class LatLonEllipsoidalBaseDI(LatLonEllipsoidalBase): 

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

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

45 ''' 

46 

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

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

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

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

51 for more details. 

52 

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

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

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

56 

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

58 

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

60 

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

62 ellipsoids are not compatible. 

63 ''' 

64 r = self._Inverse(other, wrap) 

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

66 

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

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

69 the given distance from this point along a geodesic given 

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

71 method L{destination2} for more details. 

72 

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

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

75 @kwarg height: Optional height, overriding the default 

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

77 

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

79 ''' 

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

81 

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

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

84 azimuth) after having travelled for the given distance from 

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

86 using this C{Direct} method. 

87 

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

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

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

91 

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

93 are in compass C{degrees360}. 

94 

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

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

97 

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

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

100 @kwarg height: Optional height, overriding the default 

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

102 

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

104 ''' 

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

106 return self._xnamed(r) 

107 

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

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

110 

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

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

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

114 ''' 

115 g = self.geodesic 

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

117 if LL: 

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

119 return r 

120 

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

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

123 ''' 

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

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

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

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

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

129 

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

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

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

133 L{distanceTo3} for more details. 

134 

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

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

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

138 

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

140 

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

142 

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

144 ellipsoids are not compatible. 

145 ''' 

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

147 

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

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

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

151 C{Inverse} method. 

152 

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

154 conventionally meter. The distance is measured on the surface 

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

156 

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

158 are in compass C{degrees360} from North. 

159 

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

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

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

163 

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

165 

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

167 

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

169 ellipsoids are not compatible. 

170 ''' 

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

172 

173 def finalBearingOn(self, distance, bearing): 

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

175 travelled for the given distance along a geodesic given by 

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

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

178 

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

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

181 

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

183 ''' 

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

185 

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

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

188 travelled along a geodesic from this point to an other 

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

190 L{distanceTo3} for more details. 

191 

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

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

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

195 

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

197 

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

199 

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

201 ellipsoids are not compatible. 

202 ''' 

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

204 

205 @property_RO 

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

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

208 ''' 

209 return None # PYCHOK no cover 

210 

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

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

213 ''' 

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

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

216 gl = g._DirectLine( p, other) if isscalar(other) else \ 

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

218 return g, gl, p 

219 

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

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

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

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

224 for more details. 

225 

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

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

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

229 

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

231 

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

233 

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

235 ellipsoids are not compatible. 

236 ''' 

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

238 

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

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

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

242 methods. 

243 

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

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

246 at this and 1.0 at the other point. 

247 @kwarg height: Optional height, overriding the fractional 

248 height (C{meter}). 

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

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

251 

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

253 

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

255 

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

257 

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

259 ellipsoids are not compatible. 

260 

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

262 C{rhumbMidpointTo}. 

263 ''' 

264 f = Scalar(fraction=fraction) 

265 if isnear0(f): 

266 r = self 

267 elif isnear1(f) and not wrap: 

268 r = self.others(other) 

269 else: # negative fraction OK 

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

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

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

273 return r 

274 

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

276 wrap=False, tol=_TOL): 

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

278 two points or as a point and bearing. 

279 

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

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

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

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

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

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

286 method L{Ellipsoid.geodesic_}. 

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

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

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

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

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

292 

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

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

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

296 

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

298 

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

300 or B{C{other}} invalid. 

301 

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

303 B{C{height}}. 

304 

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

306 ''' 

307 try: 

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

309 r = Radius_(circle=circle) if isscalar(circle) else \ 

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

311 

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

313 form=_MODS.dms.F_DMS) 

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

315 self.intersecant2) 

316 except (TypeError, ValueError) as x: 

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

318 exact=exact, wrap=wrap) 

319 

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

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

322 

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

324 ''' 

325 _ = self.ellipsoids(other) 

326 g = self.geodesic 

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

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

329 

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

331 equidistant=None, tol=_TOL_M): 

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

333 to this point. 

334 

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

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

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

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

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

340 is taken into account for distances. 

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

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

343 

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

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

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

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

348 

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

350 

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

352 

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

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

355 

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

357 ''' 

358 _d3 = self.distanceTo3 # Distance3Tuple 

359 _n3 = _nearestOn3 

360 try: 

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

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

363 _ = self.ellipsoids(p1) 

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

365 

366 except (TypeError, ValueError) as x: 

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

368 closed=closed, height=height, wrap=wrap) 

369 

370 # get the azimuthal equidistant projection, once 

371 A = _Equidistant00(equidistant, c) 

372 b = _Box(c, c3.distance) 

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

374 

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

376 LatLon=self.classof, # this LatLon 

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

378 try: 

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

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

381 p2 = _unrollon(p1, p2) 

382 # skip edge if no overlap with box around closest 

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

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

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

386 if d3.distance < c3.distance: 

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

388 b = _Box(c, c3.distance) 

389 m = max(m, c.iteration) 

390 p1, i = p2, j 

391 

392 except (TypeError, ValueError) as x: 

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

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

395 closed=closed, height=height, wrap=wrap) 

396 

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

398 

399 n = self.nearestOn8.__name__ 

400 c.rename(n) 

401 if s is not c: 

402 s = s.copy(name=n) 

403 if e is not c: 

404 e = e.copy(name=n) 

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

406 iteration=m) # ._iteration for tests 

407 

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

409 wrap=False, tol=_TOL): 

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

411 a geodesic from this point. 

412 

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

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

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

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

417 see method L{Ellipsoid.geodesic_}. 

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

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

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

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

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

423 

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

425 

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

427 

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

429 ''' 

430 try: 

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

432 

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

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

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

436 p._iteration = P.iteration 

437 except (TypeError, ValueError) as x: 

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

439 exact=exact, wrap=wrap) 

440 return p 

441 

442 

443class _Box(object): 

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

445 

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

447 ''' 

448 def __init__(self, center, distance): 

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

450 

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

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

453 (C{meter}, conventionally) 

454 ''' 

455 m = Radius_(distance=distance) 

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

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

458 self._N = center.lat + d 

459 self._S = center.lat - d 

460 self._E = center.lon + d 

461 self._W = center.lon - d 

462 

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

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

465 

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

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

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

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

470 

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

472 otherwise (C{bool}). 

473 ''' 

474 if lat1 > lat2: 

475 lat1, lat2 = lat2, lat1 

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

477 return False 

478 if lon1 > lon2: 

479 lon1, lon2 = lon2, lon1 

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

481 return False 

482 return True 

483 

484 

485class _Tol(object): 

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

487 ''' 

488 _deg = 0 # tol in degrees 

489 _lat = 0 

490 _m = 0 # tol in meter 

491 _min = MAX # degrees 

492 _prev = None 

493 _r = 0 

494 

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

496 '''New L{_Tol}. 

497 

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

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

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

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

502 ''' 

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

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

505 self._m = max(EPS, tol_m) 

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

507 

508 @property_RO 

509 def degrees(self): 

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

511 ''' 

512 return self._deg 

513 

514 def degrees2m(self, deg): 

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

516 ''' 

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

518 

519 def degError(self, Error=_ValueError): 

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

521 ''' 

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

523 

524 def done(self, deg): 

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

526 ''' 

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

528 return True 

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

530 self._prev = deg 

531 return False 

532 

533 @property_RO 

534 def lat(self): 

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

536 ''' 

537 return self._lat 

538 

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

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

541 ''' 

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

543 if m2km(m) > self.meter: 

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

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

546 

547 @property_RO 

548 def meter(self): 

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

550 ''' 

551 return self._m 

552 

553 @property_RO 

554 def radius(self): 

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

556 ''' 

557 return self._r 

558 

559 def reset(self): 

560 '''Reset tolerances. 

561 ''' 

562 self._min = MAX 

563 self._prev = None 

564 

565 

566def _Equidistant00(equidistant, p1): 

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

568 ''' 

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

570 equidistant = p1.Equidistant 

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

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

573 raise _IsnotError(*t, equidistant=equidistant) 

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

575 

576 

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

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

579 or as a point and (forward) bearing. 

580 

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

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

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

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

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

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

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

588 

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

590 

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

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

593 

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

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

596 ''' 

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

598 raise _IsnotError(_ellipsoidal_, center=center) 

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

600 

601 

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

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

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

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

606 ''' 

607 _LLS = _MODS.sphericalTrigonometry.LatLon 

608 _si = _MODS.sphericalTrigonometry._intersect 

609 _vi3 = _MODS.vector3d._intersect3d3 

610 

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

612 # compute opposing and distance 

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

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

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

616 return b, t.distance 

617 

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

619 # compute an end point along the initial bearing 

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

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

622 # radians in .sphericalTrigonometry._int3d2 and bearing 

623 # comparison in .sphericalTrigonometry._intb 

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

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

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

627 e = s.destination(d, e) 

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

629 

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

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

632 ll = not isscalar(e) 

633 if ll: 

634 e = s.others(**end) 

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

636 e = _unrollon(s, e) 

637 return e, ll 

638 

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

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

641 

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

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

644 if not o: # intersection may be on start 

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

646 return o 

647 return -n if b else n 

648 

649 E = s1.ellipsoids(s2) 

650 

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

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

653 

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

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

656 

657 # get the azimuthal equidistant projection 

658 A = _Equidistant00(equidistant, s1) 

659 

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

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

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

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

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

665 

666 if not ll1: 

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

668 if not ll2: 

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

670 

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

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

673 for i in range(1, _TRIPS): 

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

675 # convert start and end points to projection 

676 # space and compute an intersection there 

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

678 eps=e.meter, useZ=False) 

679 # convert intersection back to geodetic 

680 t, d = A._reverse2(v) 

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

682 t._iteration = i 

683 break 

684 else: 

685 raise e.degError(Error=IntersectionError) 

686 

687 # like .sphericalTrigonometry._intersect, if this intersection 

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

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

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

691 if b1: 

692 t = t.antipodal() 

693 b1 = not b1 

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

695 

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

697 iteration=t._iteration, name=n) 

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

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

700 

701 

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

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

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

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

706 point and an initial bearing from North. 

707 ''' 

708 s1 = _xellipsoidal(start1=start1) 

709 s2 = s1.others(start2=start2) 

710 try: 

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

712 equidistant=equidistant, tol=tol, 

713 LatLon=LatLon, **LatLon_kwds) 

714 except (TypeError, ValueError) as x: 

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

716 

717 

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

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

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

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

722 ''' 

723 c1 = _xellipsoidal(center1=center1) 

724 c2 = c1.others(center2=center2) 

725 try: 

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

727 equidistant=equidistant, tol=tol, 

728 LatLon=LatLon, **LatLon_kwds) 

729 except (TypeError, ValueError) as x: 

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

731 center2=center2, radius2=radius2) 

732 

733 

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

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

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

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

738 ''' 

739 _LLS = _MODS.sphericalTrigonometry.LatLon 

740 _si2 = _MODS.sphericalTrigonometry._intersects2 

741 _vi2 = _MODS.vector3d._intersects2 

742 

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

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

745 iteration=t.iteration, name=n) 

746 

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

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

749 

750 r1 = Radius_(radius1=radius1) 

751 r2 = Radius_(radius2=radius2) 

752 

753 E = c1.ellipsoids(c2) 

754 # get the azimuthal equidistant projection 

755 A = _Equidistant00(equidistant, c1) 

756 

757 if r1 < r2: 

758 c1, c2 = c2, c1 

759 r1, r2 = r2, r1 

760 

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

762 raise _ValueError(_exceed_PI_radians_) 

763 

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

765 c2 = _unrollon(c1, c2) 

766 

767 # distance between centers and radii are 

768 # measured along the ellipsoid's surface 

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

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

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

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

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

774 

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

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

777 

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

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

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

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

782 

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

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

785 ts, ta = [], None 

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

787 for i in range(1, _TRIPS): 

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

789 # convert centers to projection space and 

790 # compute the intersections there 

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

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

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

794 sphere=False, too_d=m) 

795 # convert intersections back to geodetic 

796 if v1 is v2: # abutting 

797 t, d = A._reverse2(v1) 

798 else: # consider the closer intersection 

799 t1, d1 = A._reverse2(v1) 

800 t2, d2 = A._reverse2(v2) 

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

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

803 t._iteration = i # _NamedTuple._iteration 

804 ts.append(t) 

805 if v1 is v2: # abutting 

806 ta = t 

807 break 

808 else: 

809 raise e.degError(Error=IntersectionError) 

810 e.reset() 

811 

812 if ta: # abutting circles 

813 pass # PYCHOK no cover 

814 elif len(ts) == 2: 

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

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

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

818 ta = ts[0] # assume abutting 

819 else: # PYCHOK no cover 

820 raise _AssertionError(ts=ts) 

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

822 return r, r 

823 

824 

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

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

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

828 separated to allow callers to embellish any exceptions. 

829 ''' 

830 p1 = p.others(point1=point1) 

831 p2 = p.others(point2=point2) 

832 

833 _ = p.ellipsoids(p1) 

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

835 

836 # get the azimuthal equidistant projection 

837 A = _Equidistant00(equidistant, p) 

838 

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

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

841 tol=tol, **LatLon_and_kwds) 

842 return NearestOn2Tuple(r, f) 

843 

844 

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

846 LatLon=None, **LatLon_kwds): 

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

848 _LLS = _MODS.sphericalNvector.LatLon 

849 _V3d = _MODS.vector3d.Vector3d 

850 _vn2 = _MODS.vector3d._nearestOn2 

851 

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

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

854 

855 E = p.ellipsoids(p2) 

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

857 

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

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

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

861 height=height) 

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

863 if height is None: 

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

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

866 h0 = favg(h1, h2) 

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

868 h0 = h1 = h2 = _0_0 

869 

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

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

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

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

874 for i in range(1, _TRIPS): 

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

876 # convert points to projection space and compute 

877 # the nearest one (and its height) there 

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

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

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

881 # convert nearest one back to geodetic 

882 t, d = A._reverse2(v) 

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

884 t._iteration = i # _NamedTuple._iteration 

885 break 

886 else: 

887 raise e.degError() 

888 

889 if height is None: 

890 h = v.z # nearest 

891 elif isscalar(height): 

892 h = height 

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

894 iteration=t.iteration, name=n) 

895 return r, f, e # fraction or None 

896 

897 

898__all__ += _ALL_DOCS(LatLonEllipsoidalBaseDI, intersecant2) 

899 

900# **) MIT License 

901# 

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

903# 

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

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

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

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

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

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

910# 

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

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

913# 

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

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

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

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

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

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

920# OTHER DEALINGS IN THE SOFTWARE.