Coverage for pygeodesy/sphericalTrigonometry.py: 93%

387 statements  

« prev     ^ index     » next       coverage.py v7.6.1, created at 2024-11-12 16:17 -0500

1 

2# -*- coding: utf-8 -*- 

3 

4u'''Spherical, C{trigonometry}-based geodesy. 

5 

6Trigonometric classes geodetic (lat-/longitude) L{LatLon} and 

7geocentric (ECEF) L{Cartesian} and functions L{areaOf}, L{intersection}, 

8L{intersections2}, L{isPoleEnclosedBy}, L{meanOf}, L{nearestOn3} and 

9L{perimeterOf}, I{all spherical}. 

10 

11Pure Python implementation of geodetic (lat-/longitude) methods using 

12spherical trigonometry, transcoded from JavaScript originals by 

13I{(C) Chris Veness 2011-2024} published under the same MIT Licence**, see 

14U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}. 

15''' 

16# make sure int/int division yields float quotient, see .basics 

17from __future__ import division as _; del _ # PYCHOK semicolon 

18 

19from pygeodesy.basics import copysign0, map1, signOf 

20from pygeodesy.constants import EPS, EPS1, EPS4, PI, PI2, PI_2, PI_4, R_M, \ 

21 isnear0, isnear1, isnon0, _0_0, _0_5, \ 

22 _1_0, _2_0, _90_0 

23from pygeodesy.datums import _ellipsoidal_datum, _mean_radius 

24from pygeodesy.errors import _AssertionError, CrossError, crosserrors, \ 

25 _TypeError, _ValueError, IntersectionError, \ 

26 _xError, _xkwds, _xkwds_get, _xkwds_pop2 

27from pygeodesy.fmath import favg, fdot, fdot_, fmean, hypot 

28from pygeodesy.fsums import Fsum, fsum, fsumf_ 

29from pygeodesy.formy import antipode_, bearing_, _bearingTo2, excessAbc_, \ 

30 excessGirard_, excessLHuilier_, opposing_, _radical2, \ 

31 vincentys_ 

32from pygeodesy.interns import _1_, _2_, _coincident_, _composite_, _colinear_, \ 

33 _concentric_, _convex_, _end_, _infinite_, \ 

34 _invalid_, _line_, _near_, _null_, _parallel_, \ 

35 _point_, _SPACE_, _too_ 

36from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER 

37# from pygeodesy.nvectorBase import NvectorBase, sumOf # _MODE 

38from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, NearestOn3Tuple, \ 

39 Triangle7Tuple, Triangle8Tuple 

40from pygeodesy.points import ispolar, nearestOn5 as _nearestOn5, \ 

41 Fmt as _Fmt # XXX shadowed 

42from pygeodesy.props import deprecated_function, deprecated_method 

43from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \ 

44 _intersecant2, LatLonSphericalBase, \ 

45 _rads3, _radians2m, _trilaterate5 

46# from pygeodesy.streprs import Fmt as _Fmt # from .points XXX shadowed 

47from pygeodesy.units import Bearing_, Height, _isDegrees, _isRadius, Lamd, \ 

48 Phid, Radius_, Scalar 

49from pygeodesy.utily import acos1, asin1, atan1d, atan2d, degrees90, degrees180, \ 

50 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \ 

51 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI 

52from pygeodesy.vector3d import sumOf, Vector3d 

53 

54from math import asin, atan2, cos, degrees, fabs, radians, sin 

55 

56__all__ = _ALL_LAZY.sphericalTrigonometry 

57__version__ = '24.11.07' 

58 

59_PI_EPS4 = PI - EPS4 

60if _PI_EPS4 >= PI: 

61 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4) 

62 

63 

64class Cartesian(CartesianSphericalBase): 

65 '''Extended to convert geocentric, L{Cartesian} points to 

66 spherical, geodetic L{LatLon}. 

67 ''' 

68 

69 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon 

70 '''Convert this cartesian point to a C{spherical} geodetic point. 

71 

72 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword 

73 arguments. Use C{B{LatLon}=...} to override 

74 this L{LatLon} class or specify C{B{LatLon}=None}. 

75 

76 @return: The geodetic point (L{LatLon}) or if C{B{LatLon} is None}, 

77 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} 

78 with C{C} and C{M} if available. 

79 

80 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument. 

81 ''' 

82 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum) 

83 return CartesianSphericalBase.toLatLon(self, **kwds) 

84 

85 

86class LatLon(LatLonSphericalBase): 

87 '''New point on a spherical earth model, based on trigonometry formulae. 

88 ''' 

89 

90 def _ab1_ab2_db5(self, other, wrap): 

91 '''(INTERNAL) Helper for several methods. 

92 ''' 

93 a1, b1 = self.philam 

94 a2, b2 = self.others(other, up=2).philam 

95 if wrap: 

96 a2, b2 = _Wrap.philam(a2, b2) 

97 db, b2 = unrollPI(b1, b2, wrap=wrap) 

98 else: # unrollPI shortcut 

99 db = b2 - b1 

100 return a1, b1, a2, b2, db 

101 

102 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False): 

103 '''Compute the (signed) distance from the start to the closest 

104 point on the great circle line defined by a start and an 

105 end point. 

106 

107 That is, if a perpendicular is drawn from this point to the 

108 great circle line, the along-track distance is the distance 

109 from the start point to the point where the perpendicular 

110 crosses the line. 

111 

112 @arg start: Start point of the great circle line (L{LatLon}). 

113 @arg end: End point of the great circle line (L{LatLon}). 

114 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 

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

116 the B{C{start}} and B{C{end}} point (C{bool}). 

117 

118 @return: Distance along the great circle line (C{radians} 

119 if C{B{radius} is None} or C{meter}, same units 

120 as B{C{radius}}), positive if I{after} the 

121 B{C{start}} toward the B{C{end}} point of the 

122 line, I{negative} if before or C{0} if at the 

123 B{C{start}} point. 

124 

125 @raise TypeError: Invalid B{C{start}} or B{C{end}} point. 

126 

127 @raise ValueError: Invalid B{C{radius}}. 

128 ''' 

129 r, x, b = self._a_x_b3(start, end, radius, wrap) 

130 cx = cos(x) 

131 return _0_0 if isnear0(cx) else \ 

132 _radians2m(copysign0(acos1(cos(r) / cx), cos(b)), radius) 

133 

134 def _a_x_b3(self, start, end, radius, wrap): 

135 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo. 

136 ''' 

137 s = self.others(start=start) 

138 e = self.others(end=end) 

139 s, e, w = _unrollon3(self, s, e, wrap) 

140 

141 r = Radius_(radius) 

142 r = s.distanceTo(self, r, wrap=w) / r 

143 

144 b = radians(s.initialBearingTo(self, wrap=w) 

145 - s.initialBearingTo(e, wrap=w)) 

146 x = asin(sin(r) * sin(b)) 

147 return r, x, -b 

148 

149 @deprecated_method 

150 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover 

151 '''DEPRECATED, use method L{initialBearingTo}. 

152 ''' 

153 return self.initialBearingTo(other, wrap=wrap, raiser=raiser) 

154 

155 def crossingParallels(self, other, lat, wrap=False): 

156 '''Return the pair of meridians at which a great circle defined 

157 by this and an other point crosses the given latitude. 

158 

159 @arg other: The other point defining great circle (L{LatLon}). 

160 @arg lat: Latitude at the crossing (C{degrees}). 

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

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

163 

164 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or 

165 C{None} if the great circle doesn't reach B{C{lat}}. 

166 ''' 

167 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap) 

168 sa, ca, sa1, ca1, \ 

169 sa2, ca2, sdb, cdb = sincos2_(radians(lat), a1, a2, db) 

170 sa1 *= ca2 * ca 

171 

172 x = sa1 * sdb 

173 y = sa1 * cdb - ca1 * sa2 * ca 

174 z = ca1 * sdb * ca2 * sa 

175 

176 h = hypot(x, y) 

177 if h < EPS or fabs(z) > h: # PYCHOK no cover 

178 return None # great circle doesn't reach latitude 

179 

180 m = atan2(-y, x) + b1 # longitude at max latitude 

181 d = acos1(z / h) # delta longitude to intersections 

182 return degrees180(m - d), degrees180(m + d) 

183 

184 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False): 

185 '''Compute the (signed) distance from this point to a great 

186 circle from a start to an end point. 

187 

188 @arg start: Start point of the great circle line (L{LatLon}). 

189 @arg end: End point of the great circle line (L{LatLon}). 

190 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 

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

192 the B{C{start}} and B{C{end}} point (C{bool}). 

193 

194 @return: Distance to the great circle (C{radians} if 

195 B{C{radius}} or C{meter}, same units as 

196 B{C{radius}}), I{negative} if to the left or 

197 I{positive} if to the right of the line. 

198 

199 @raise TypeError: If B{C{start}} or B{C{end}} is not L{LatLon}. 

200 

201 @raise ValueError: Invalid B{C{radius}}. 

202 ''' 

203 _, x, _ = self._a_x_b3(start, end, radius, wrap) 

204 return _radians2m(x, radius) 

205 

206 def destination(self, distance, bearing, radius=R_M, height=None): 

207 '''Locate the destination from this point after having 

208 travelled the given distance on a bearing from North. 

209 

210 @arg distance: Distance travelled (C{meter}, same units as 

211 B{C{radius}}). 

212 @arg bearing: Bearing from this point (compass C{degrees360}). 

213 @kwarg radius: Mean earth radius (C{meter}). 

214 @kwarg height: Optional height at destination (C{meter}, same 

215 units a B{C{radius}}). 

216 

217 @return: Destination point (L{LatLon}). 

218 

219 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}}, 

220 B{C{radius}} or B{C{height}}. 

221 ''' 

222 a, b = self.philam 

223 r, t = _m2radians(distance, radius, low=None), Bearing_(bearing) 

224 

225 a, b = _destination2(a, b, r, t) 

226 h = self._heigHt(height) 

227 return self.classof(degrees90(a), degrees180(b), height=h) 

228 

229 def distanceTo(self, other, radius=R_M, wrap=False): 

230 '''Compute the (angular) distance from this to an other point. 

231 

232 @arg other: The other point (L{LatLon}). 

233 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 

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

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

236 

237 @return: Distance between this and the B{C{other}} point 

238 (C{meter}, same units as B{C{radius}} or 

239 C{radians} if C{B{radius} is None}). 

240 

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

242 

243 @raise ValueError: Invalid B{C{radius}}. 

244 ''' 

245 a1, _, a2, _, db = self._ab1_ab2_db5(other, wrap) 

246 return _radians2m(vincentys_(a2, a1, db), radius) 

247 

248# @Property_RO 

249# def Ecef(self): 

250# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}. 

251# ''' 

252# return _MODS.ecef.EcefKarney 

253 

254 def greatCircle(self, bearing, Vector=Vector3d, **Vector_kwds): 

255 '''Compute the vector normal to great circle obtained by heading 

256 from this point on the bearing from North. 

257 

258 Direction of vector is such that initial bearing vector 

259 b = c × n, where n is an n-vector representing this point. 

260 

261 @arg bearing: Bearing from this point (compass C{degrees360}). 

262 @kwarg Vector: Vector class to return the great circle, 

263 overriding the default L{Vector3d}. 

264 @kwarg Vector_kwds: Optional, additional keyword argunents 

265 for B{C{Vector}}. 

266 

267 @return: Vector representing great circle (C{Vector}). 

268 

269 @raise ValueError: Invalid B{C{bearing}}. 

270 ''' 

271 a, b = self.philam 

272 sa, ca, sb, cb, st, ct = sincos2_(a, b, Bearing_(bearing)) 

273 

274 sa *= st 

275 return Vector(fdot_(sb, ct, -cb, sa), 

276 -fdot_(cb, ct, sb, sa), 

277 ca * st, **Vector_kwds) # XXX .unit()? 

278 

279 def initialBearingTo(self, other, wrap=False, raiser=False): 

280 '''Compute the initial bearing (forward azimuth) from this 

281 to an other point. 

282 

283 @arg other: The other point (spherical L{LatLon}). 

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

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

286 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}), 

287 use C{B{raiser}=True} for behavior like 

288 C{sphericalNvector.LatLon.initialBearingTo}. 

289 

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

291 

292 @raise CrossError: If this and the B{C{other}} point coincide 

293 and if B{C{raiser}} and L{crosserrors 

294 <pygeodesy.crosserrors>} are both C{True}. 

295 

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

297 ''' 

298 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap) 

299 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

300 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(db)) < EPS: 

301 raise CrossError(_point_, self, other=other, wrap=wrap, txt=_coincident_) 

302 

303 return degrees(bearing_(a1, b1, a2, b2, final=False)) 

304 

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

306 '''Locate the point at given fraction between (or along) this 

307 and an other point. 

308 

309 @arg other: The other point (L{LatLon}). 

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

311 0.0 at this and 1.0 at the other point). 

312 @kwarg height: Optional height, overriding the intermediate 

313 height (C{meter}). 

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

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

316 

317 @return: Intermediate point (L{LatLon}). 

318 

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

320 

321 @raise ValueError: Invalid B{C{fraction}} or B{C{height}}. 

322 

323 @see: Methods C{midpointTo} and C{rhumbMidpointTo}. 

324 ''' 

325 p = self 

326 f = Scalar(fraction=fraction) 

327 if not isnear0(f): 

328 p = p.others(other) 

329 if wrap: 

330 p = _Wrap.point(p) 

331 if not isnear1(f): # and not near0 

332 a1, b1 = self.philam 

333 a2, b2 = p.philam 

334 db, b2 = unrollPI(b1, b2, wrap=wrap) 

335 r = vincentys_(a2, a1, db) 

336 sr = sin(r) 

337 if isnon0(sr): 

338 sa1, ca1, sa2, ca2, \ 

339 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2) 

340 

341 t = f * r 

342 a = sin(r - t) # / sr superflous 

343 b = sin( t) # / sr superflous 

344 

345 x = fdot_(a, ca1 * cb1, b, ca2 * cb2) 

346 y = fdot_(a, ca1 * sb1, b, ca2 * sb2) 

347 z = fdot_(a, sa1, b, sa2) 

348 

349 a = atan1d(z, hypot(x, y)) 

350 b = atan2d(y, x) 

351 

352 else: # PYCHOK no cover 

353 a = degrees90( favg(a1, a2, f=f)) # coincident 

354 b = degrees180(favg(b1, b2, f=f)) 

355 

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

357 p = self.classof(a, b, height=h) 

358 return p 

359 

360 def intersection(self, end1, other, end2, height=None, wrap=False): 

361 '''Compute the intersection point of two lines, each defined by 

362 two points or a start point and a bearing from North. 

363 

364 @arg end1: End point of this line (L{LatLon}) or the initial 

365 bearing at this point (compass C{degrees360}). 

366 @arg other: Start point of the other line (L{LatLon}). 

367 @arg end2: End point of the other line (L{LatLon}) or the 

368 initial bearing at the B{C{other}} point (compass 

369 C{degrees360}). 

370 @kwarg height: Optional height for intersection point, 

371 overriding the mean height (C{meter}). 

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

373 B{C{start2}} and both B{C{end*}} points (C{bool}). 

374 

375 @return: The intersection point (L{LatLon}). An alternate 

376 intersection point might be the L{antipode} to 

377 the returned result. 

378 

379 @raise IntersectionError: Ambiguous or infinite intersection 

380 or colinear, parallel or otherwise 

381 non-intersecting lines. 

382 

383 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}} 

384 or B{C{end2}} not C{scalar} nor L{LatLon}. 

385 

386 @raise ValueError: Invalid B{C{height}} or C{null} line. 

387 ''' 

388 try: 

389 s2 = self.others(other) 

390 return _intersect(self, end1, s2, end2, height=height, wrap=wrap, 

391 LatLon=self.classof) 

392 except (TypeError, ValueError) as x: 

393 raise _xError(x, start1=self, end1=end1, 

394 other=other, end2=end2, wrap=wrap) 

395 

396 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0, 

397 height=None, wrap=True): 

398 '''Compute the intersection points of two circles, each defined 

399 by a center point and a radius. 

400 

401 @arg rad1: Radius of the this circle (C{meter} or C{radians}, 

402 see B{C{radius}}). 

403 @arg other: Center point of the other circle (L{LatLon}). 

404 @arg rad2: Radius of the other circle (C{meter} or C{radians}, 

405 see B{C{radius}}). 

406 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}}, 

407 B{C{rad2}} and B{C{eps}} are given in C{radians}). 

408 @kwarg eps: Required overlap (C{meter} or C{radians}, see 

409 B{C{radius}}). 

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

411 conventionally) or C{None} for the I{"radical height"} 

412 at the I{radical line} between both centers. 

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

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

415 

416 @return: 2-Tuple of the intersection points, each a L{LatLon} 

417 instance. For abutting circles, both intersection 

418 points are the same instance, aka the I{radical center}. 

419 

420 @raise IntersectionError: Concentric, antipodal, invalid or 

421 non-intersecting circles. 

422 

423 @raise TypeError: If B{C{other}} is not L{LatLon}. 

424 

425 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}}, 

426 B{C{eps}} or B{C{height}}. 

427 ''' 

428 try: 

429 c2 = self.others(other) 

430 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps, 

431 height=height, wrap=wrap, 

432 LatLon=self.classof) 

433 except (TypeError, ValueError) as x: 

434 raise _xError(x, center=self, rad1=rad1, 

435 other=other, rad2=rad2, wrap=wrap) 

436 

437 @deprecated_method 

438 def isEnclosedBy(self, points): # PYCHOK no cover 

439 '''DEPRECATED, use method C{isenclosedBy}.''' 

440 return self.isenclosedBy(points) 

441 

442 def isenclosedBy(self, points, wrap=False): 

443 '''Check whether a (convex) polygon or composite encloses this point. 

444 

445 @arg points: The polygon points or composite (L{LatLon}[], 

446 L{BooleanFHP} or L{BooleanGH}). 

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

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

449 

450 @return: C{True} if this point is inside the polygon or 

451 composite, C{False} otherwise. 

452 

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

454 

455 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

456 

457 @raise ValueError: Invalid B{C{points}}, non-convex polygon. 

458 

459 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy} 

460 and L{pygeodesy.ispolar} especially if the B{C{points}} may 

461 enclose a pole or wrap around the earth I{longitudinally}. 

462 ''' 

463 if _MODS.booleans.isBoolean(points): 

464 return points._encloses(self.lat, self.lon, wrap=wrap) 

465 

466 Ps = self.PointsIter(points, loop=2, dedup=True, wrap=wrap) 

467 n0 = self._N_vector 

468 

469 v2 = Ps[0]._N_vector 

470 p1 = Ps[1] 

471 v1 = p1._N_vector 

472 # check whether this point on same side of all 

473 # polygon edges (to the left or right depending 

474 # on the anti-/clockwise polygon direction) 

475 gc1 = v2.cross(v1) 

476 t0 = gc1.angleTo(n0) > PI_2 

477 s0 = None 

478 # get great-circle vector for each edge 

479 for i, p2 in Ps.enumerate(closed=True): 

480 if wrap and not Ps.looped: 

481 p2 = _unrollon(p1, p2) 

482 p1 = p2 

483 v2 = p2._N_vector 

484 gc = v1.cross(v2) 

485 t = gc.angleTo(n0) > PI_2 

486 if t != t0: # different sides of edge i 

487 return False # outside 

488 

489 # check for convex polygon: angle between 

490 # gc vectors, signed by direction of n0 

491 # (otherwise the test above is not reliable) 

492 s = signOf(gc1.angleTo(gc, vSign=n0)) 

493 if s != s0: 

494 if s0 is None: 

495 s0 = s 

496 else: 

497 t = _Fmt.SQUARE(points=i) 

498 raise _ValueError(t, p2, wrap=wrap, txt_not_=_convex_) 

499 gc1, v1 = gc, v2 

500 

501 return True # inside 

502 

503 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False): 

504 '''Find the midpoint between this and an other point. 

505 

506 @arg other: The other point (L{LatLon}). 

507 @kwarg height: Optional height for midpoint, overriding 

508 the mean height (C{meter}). 

509 @kwarg fraction: Midpoint location from this point (C{scalar}), 

510 may be negative or greater than 1.0. 

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

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

513 

514 @return: Midpoint (L{LatLon}). 

515 

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

517 

518 @raise ValueError: Invalid B{C{height}}. 

519 

520 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}. 

521 ''' 

522 if fraction is _0_5: 

523 # see <https://MathForum.org/library/drmath/view/51822.html> 

524 a1, b, a2, _, db = self._ab1_ab2_db5(other, wrap) 

525 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db) 

526 

527 x = ca2 * cdb + ca1 

528 y = ca2 * sdb 

529 

530 a = atan1d(sa1 + sa2, hypot(x, y)) 

531 b = degrees180(b + atan2(y, x)) 

532 

533 h = self._havg(other, h=height) 

534 r = self.classof(a, b, height=h) 

535 else: 

536 r = self.intermediateTo(other, fraction, height=height, wrap=wrap) 

537 return r 

538 

539 def nearestOn(self, point1, point2, radius=R_M, **wrap_adjust_limit): 

540 '''Locate the point between two other points closest to this point. 

541 

542 Distances are approximated by function L{pygeodesy.equirectangular4}, 

543 subject to the supplied B{C{options}}. 

544 

545 @arg point1: Start point (L{LatLon}). 

546 @arg point2: End point (L{LatLon}). 

547 @kwarg radius: Mean earth radius (C{meter}). 

548 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

549 L{sphericalTrigonometry.nearestOn3} and 

550 L{pygeodesy.equirectangular4}, 

551 

552 @return: Closest point on the great circle line (L{LatLon}). 

553 

554 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}}, 

555 see function L{pygeodesy.equirectangular4}. 

556 

557 @raise NotImplementedError: Keyword argument C{B{within}=False} 

558 is not (yet) supported. 

559 

560 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}. 

561 

562 @raise ValueError: Invalid B{C{radius}} or B{C{options}}. 

563 

564 @see: Functions L{pygeodesy.equirectangular4} and L{pygeodesy.nearestOn5} 

565 and method L{sphericalTrigonometry.LatLon.nearestOn3}. 

566 ''' 

567 # remove kwarg B{C{within}} if present 

568 w, kwds = _xkwds_pop2(wrap_adjust_limit, within=True) 

569 if not w: 

570 self._notImplemented(within=w) 

571 

572# # UNTESTED - handle C{B{within}=False} and C{B{within}=True} 

573# wrap = _xkwds_get(options, wrap=False) 

574# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap) 

575# if fabs(a) < EPS or (within and a < EPS): 

576# return point1 

577# d = point1.distanceTo(point2, radius=radius, wrap=wrap) 

578# if isnear0(d): 

579# return point1 # or point2 

580# elif fabs(d - a) < EPS or (a + EPS) > d: 

581# return point2 

582# f = a / d 

583# if within: 

584# if f > EPS1: 

585# return point2 

586# elif f < EPS: 

587# return point1 

588# return point1.intermediateTo(point2, f, wrap=wrap) 

589 

590 # without kwarg B{C{within}}, use backward compatible .nearestOn3 

591 return self.nearestOn3([point1, point2], closed=False, radius=radius, 

592 **kwds)[0] 

593 

594 @deprecated_method 

595 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover 

596 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}. 

597 

598 @return: ... 2-Tuple C{(closest, distance)} of the closest 

599 point (L{LatLon}) on the polygon and the distance 

600 to that point from this point in C{meter}, same 

601 units of B{C{radius}}. 

602 ''' 

603 r = self.nearestOn3(points, closed=closed, radius=radius, **options) 

604 return r.closest, r.distance 

605 

606 def nearestOn3(self, points, closed=False, radius=R_M, **wrap_adjust_limit): 

607 '''Locate the point on a polygon closest to this point. 

608 

609 Distances are approximated by function L{pygeodesy.equirectangular4}, 

610 subject to the supplied B{C{options}}. 

611 

612 @arg points: The polygon points (L{LatLon}[]). 

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

614 @kwarg radius: Mean earth radius (C{meter}). 

615 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

616 L{sphericalTrigonometry.nearestOn3} and 

617 L{pygeodesy.equirectangular4}, 

618 

619 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the 

620 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular4} 

621 C{distance} between this and the C{closest} point converted to 

622 C{meter}, same units as B{C{radius}}. The C{angle} from this 

623 to the C{closest} point is in compass C{degrees360}, like 

624 function L{pygeodesy.compassAngle}. 

625 

626 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}}, 

627 see function L{pygeodesy.equirectangular4}. 

628 

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

630 

631 @raise TypeError: Some B{C{points}} are not C{LatLon}. 

632 

633 @raise ValueError: Invalid B{C{radius}} or B{C{options}}. 

634 

635 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular4} 

636 and L{pygeodesy.nearestOn5}. 

637 ''' 

638 return nearestOn3(self, points, closed=closed, radius=radius, 

639 LatLon=self.classof, **wrap_adjust_limit) 

640 

641 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None 

642 '''Convert this point to C{Karney}-based cartesian (ECEF) coordinates. 

643 

644 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}} and other 

645 keyword arguments, ignored if C{B{Cartesian} is 

646 None}. Use C{B{Cartesian}=...} to override this 

647 L{Cartesian} class or specify C{B{Cartesian}=None}. 

648 

649 @return: The cartesian point (L{Cartesian}) or if C{B{Cartesian} is None}, 

650 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with C{C} 

651 and C{M} if available. 

652 

653 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument. 

654 ''' 

655 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum) 

656 return LatLonSphericalBase.toCartesian(self, **kwds) 

657 

658 def triangle7(self, otherB, otherC, radius=R_M, wrap=False): 

659 '''Compute the angles, sides and area of a spherical triangle. 

660 

661 @arg otherB: Second triangle point (C{LatLon}). 

662 @arg otherC: Third triangle point (C{LatLon}). 

663 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, L{Ellipsoid}, 

664 L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) or C{None}. 

665 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll points B{C{otherB}} 

666 and B{C{otherC}} (C{bool}). 

667 

668 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if B{C{radius} is 

669 None}, a L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)}. 

670 

671 @see: Function L{triangle7} and U{Spherical trigonometry 

672 <https://WikiPedia.org/wiki/Spherical_trigonometry>}. 

673 ''' 

674 B = self.others(otherB=otherB) 

675 C = self.others(otherC=otherC) 

676 B, C, _ = _unrollon3(self, B, C, wrap) 

677 

678 r = self.philam + B.philam + C.philam 

679 t = triangle8_(*r, wrap=wrap) 

680 return self._xnamed(_t7Tuple(t, radius)) 

681 

682 def triangulate(self, bearing1, other, bearing2, **height_wrap): 

683 '''Locate a point given this, an other point and a bearing from 

684 North at both points. 

685 

686 @arg bearing1: Bearing at this point (compass C{degrees360}). 

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

688 @arg bearing2: Bearing at the other point (compass C{degrees360}). 

689 @kwarg height_wrap_tol: Optional keyword arguments C{B{height}=None}, 

690 C{B{wrap}=False}, see method L{intersection}. 

691 

692 @return: Triangulated point (C{LatLon}). 

693 

694 @see: Method L{intersection} for further details. 

695 ''' 

696 if _isDegrees(bearing1) and _isDegrees(bearing2): 

697 return self.intersection(bearing1, other, bearing2, **height_wrap) 

698 raise _TypeError(bearing1=bearing1, bearing2=bearing2, **height_wrap) 

699 

700 def trilaterate5(self, distance1, point2, distance2, point3, distance3, 

701 area=True, eps=EPS1, radius=R_M, wrap=False): 

702 '''Trilaterate three points by I{area overlap} or I{perimeter intersection} 

703 of three corresponding circles. 

704 

705 @arg distance1: Distance to this point (C{meter}, same units as B{C{radius}}). 

706 @arg point2: Second center point (C{LatLon}). 

707 @arg distance2: Distance to point2 (C{meter}, same units as B{C{radius}}). 

708 @arg point3: Third center point (C{LatLon}). 

709 @arg distance3: Distance to point3 (C{meter}, same units as B{C{radius}}). 

710 @kwarg area: If C{True}, compute the area overlap, otherwise the perimeter 

711 intersection of the circles (C{bool}). 

712 @kwarg eps: The required I{minimal overlap} for C{B{area}=True} or the 

713 I{intersection margin} if C{B{area}=False} (C{meter}, same 

714 units as B{C{radius}}). 

715 @kwarg radius: Mean earth radius (C{meter}, conventionally). 

716 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{point2}} and 

717 B{C{point3}} (C{bool}). 

718 

719 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)} with 

720 C{min} and C{max} in C{meter}, same units as B{C{eps}}, the 

721 corresponding trilaterated points C{minPoint} and C{maxPoint} 

722 as I{spherical} C{LatLon} and C{n}, the number of trilatered 

723 points found for the given B{C{eps}}. 

724 

725 If only a single trilaterated point is found, C{min I{is} max}, 

726 C{minPoint I{is} maxPoint} and C{n = 1}. 

727 

728 For C{B{area}=True}, C{min} and C{max} are the smallest respectively 

729 largest I{radial} overlap found. 

730 

731 For C{B{area}=False}, C{min} and C{max} represent the nearest 

732 respectively farthest intersection margin. 

733 

734 If C{B{area}=True} and all 3 circles are concentric, C{n=0} and 

735 C{minPoint} and C{maxPoint} are both the B{C{point#}} with the 

736 smallest B{C{distance#}} C{min} and C{max} the largest B{C{distance#}}. 

737 

738 @raise IntersectionError: Trilateration failed for the given B{C{eps}}, 

739 insufficient overlap for C{B{area}=True} or 

740 no intersection or all (near-)concentric if 

741 C{B{area}=False}. 

742 

743 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}. 

744 

745 @raise ValueError: Coincident B{C{point2}} or B{C{point3}} or invalid 

746 B{C{distance1}}, B{C{distance2}}, B{C{distance3}} 

747 or B{C{radius}}. 

748 ''' 

749 return _trilaterate5(self, distance1, 

750 self.others(point2=point2), distance2, 

751 self.others(point3=point3), distance3, 

752 area=area, radius=radius, eps=eps, wrap=wrap) 

753 

754 

755_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon}) 

756 

757 

758def areaOf(points, radius=R_M, wrap=False): # was=True 

759 '''Calculate the area of a (spherical) polygon or composite (with the 

760 points joined by great circle arcs). 

761 

762 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP} 

763 or L{BooleanGH}). 

764 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, 

765 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) 

766 or C{None}. 

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

768 (C{bool}). 

769 

770 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}} 

771 or C{radians} if C{B{radius} is None}). 

772 

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

774 

775 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

776 

777 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge. 

778 

779 @note: The area is based on I{Karney}'s U{'Area of a spherical 

780 polygon'<https://MathOverflow.net/questions/97711/ 

781 the-area-of-spherical-polygons>}, 3rd Answer. 

782 

783 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf}, 

784 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and 

785 L{ellipsoidalKarney.areaOf}. 

786 ''' 

787 if _MODS.booleans.isBoolean(points): 

788 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap) 

789 

790 _at2, _t_2 = atan2, tan_2 

791 _un, _w180 = unrollPI, wrap180 

792 

793 Ps = _T00.PointsIter(points, loop=1, wrap=wrap) 

794 p1 = p2 = Ps[0] 

795 a1, b1 = p1.philam 

796 ta1, z1 = _t_2(a1), None 

797 

798 A = Fsum() # mean phi 

799 R = Fsum() # see L{pygeodesy.excessKarney_} 

800 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360° 

801 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html> 

802 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points 

803 D = Fsum() 

804 for i, p2 in Ps.enumerate(closed=True): 

805 a2, b2 = p2.philam 

806 db, b2 = _un(b1, b2, wrap=wrap and not Ps.looped) 

807 A += a2 

808 ta2 = _t_2(a2) 

809 tdb = _t_2(db, points=i) 

810 R += _at2(tdb * (ta1 + ta2), 

811 _1_0 + ta1 * ta2) 

812 ta1, b1 = ta2, b2 

813 

814 if not p2.isequalTo(p1, eps=EPS): 

815 z, z2 = _bearingTo2(p1, p2, wrap=wrap) 

816 if z1 is not None: 

817 D += _w180(z - z1) # (z - z1 + 540) ... 

818 D += _w180(z2 - z) # (z2 - z + 540) % 360 - 180 

819 p1, z1 = p2, z2 

820 

821 R = abs(R * _2_0) 

822 if abs(D) < _90_0: # ispolar(points) 

823 R = abs(R - PI2) 

824 if radius: 

825 a = degrees(A.fover(len(A))) # mean lat 

826 R *= _mean_radius(radius, a)**2 

827 return float(R) 

828 

829 

830def _destination2(a, b, r, t): 

831 '''(INTERNAL) Destination lat- and longitude in C{radians}. 

832 

833 @arg a: Latitude (C{radians}). 

834 @arg b: Longitude (C{radians}). 

835 @arg r: Angular distance (C{radians}). 

836 @arg t: Bearing (compass C{radians}). 

837 

838 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}). 

839 ''' 

840 # see <https://www.EdWilliams.org/avform.htm#LL> 

841 sa, ca, sr, cr, st, ct = sincos2_(a, r, t) 

842 ca *= sr 

843 

844 a = asin1(ct * ca + cr * sa) 

845 d = atan2(st * ca, cr - sa * sin(a)) 

846 # note, in EdWilliams.org/avform.htm W is + and E is - 

847 return a, (b + d) # (mod(b + d + PI, PI2) - PI) 

848 

849 

850def _int3d2(s, end, wrap, _i_, Vector, hs): 

851 # see <https://www.EdWilliams.org/intersect.htm> (5) ff 

852 # and similar logic in .ellipsoidalBaseDI._intersect3 

853 a1, b1 = s.philam 

854 

855 if _isDegrees(end): # bearing, get pseudo-end point 

856 a2, b2 = _destination2(a1, b1, PI_4, radians(end)) 

857 else: # must be a point 

858 s.others(end, name=_end_ + _i_) 

859 hs.append(end.height) 

860 a2, b2 = end.philam 

861 if wrap: 

862 a2, b2 = _Wrap.philam(a2, b2) 

863 

864 db, b2 = unrollPI(b1, b2, wrap=wrap) 

865 if max(fabs(db), fabs(a2 - a1)) < EPS: 

866 raise _ValueError(_SPACE_(_line_ + _i_, _null_)) 

867 # note, in EdWilliams.org/avform.htm W is + and E is - 

868 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5, 

869 -(b1 + b2) * _0_5) 

870 cb21 *= sin(a1 - a2) # sa21 

871 sb21 *= sin(a1 + a2) # sa12 

872 x = Vector(fdot_(sb12, cb21, -cb12, sb21), 

873 fdot_(cb12, cb21, sb12, sb21), 

874 cos(a1) * cos(a2) * sin(db)) # ll=start 

875 return x.unit(), (db, (a2 - a1)) # negated d 

876 

877 

878def _intdot(ds, a1, b1, a, b, wrap): 

879 # compute dot product ds . (-b + b1, a - a1) 

880 db, _ = unrollPI(b1, b, wrap=wrap) 

881 return fdot(ds, db, a - a1) 

882 

883 

884def intersecant2(center, circle, point, other, **radius_exact_height_wrap): 

885 '''Compute the intersections of a circle and a (great circle) line given as 

886 two points or as a point and a bearing from North. 

887 

888 @arg center: Center of the circle (L{LatLon}). 

889 @arg circle: Radius of the circle (C{meter}, same units as the earth 

890 B{C{radius}}) or a point on the circle (L{LatLon}). 

891 @arg point: A point on the (great circle) line (L{LatLon}). 

892 @arg other: An other point on the (great circle) line (L{LatLon}) or 

893 the bearing at the B{C{point}} (compass C{degrees360}). 

894 @kwarg radius_exact_height_wrap: Optional keyword arguments, see method 

895 L{intersecant2<pygeodesy.sphericalBase.LatLonSphericalBase. 

896 intersecant2>} for further details. 

897 

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

899 an instance of the B{C{point}} class. Both points are the same 

900 instance if the (great circle) line is tangent to the circle. 

901 

902 @raise IntersectionError: The circle and line do not intersect. 

903 

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

905 not L{LatLon}. 

906 

907 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}}, 

908 B{C{exact}}, B{C{height}} or B{C{napieradius}}. 

909 ''' 

910 c = _T00.others(center=center) 

911 p = _T00.others(point=point) 

912 try: 

913 return _intersecant2(c, circle, p, other, **radius_exact_height_wrap) 

914 except (TypeError, ValueError) as x: 

915 raise _xError(x, center=center, circle=circle, point=point, other=other, 

916 **radius_exact_height_wrap) 

917 

918 

919def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3 

920 LatLon=LatLon, **LatLon_kwds): 

921 # (INTERNAL) Intersect two (spherical) lines, see L{intersection} 

922 # above, separated to allow callers to embellish any exceptions 

923 

924 s1, s2 = start1, start2 

925 if wrap: 

926 s2 = _Wrap.point(s2) 

927 hs = [s1.height, s2.height] 

928 

929 a1, b1 = s1.philam 

930 a2, b2 = s2.philam 

931 db, b2 = unrollPI(b1, b2, wrap=wrap) 

932 r12 = vincentys_(a2, a1, db) 

933 if fabs(r12) < EPS: # [nearly] coincident points 

934 a, b = favg(a1, a2), favg(b1, b2) 

935 

936 # see <https://www.EdWilliams.org/avform.htm#Intersection> 

937 elif _isDegrees(end1) and _isDegrees(end2): # both bearings 

938 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12) 

939 

940 x1, x2 = (sr12 * ca1), (sr12 * ca2) 

941 if isnear0(x1) or isnear0(x2): 

942 raise IntersectionError(_parallel_) 

943 # handle domain error for equivalent longitudes, 

944 # see also functions asin_safe and acos_safe at 

945 # <https://www.EdWilliams.org/avform.htm#Math> 

946 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1) 

947 t21, t23 = acos1((sa1 - sa2 * cr12) / x2), radiansPI2(end2) 

948 if sin(db) > 0: 

949 t21 = PI2 - t21 

950 else: 

951 t12 = PI2 - t12 

952 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3 

953 wrapPI(t21 - t23)) # angle 1-2-3) 

954 if isnear0(sx1) and isnear0(sx2): 

955 raise IntersectionError(_infinite_) 

956 sx3 = sx1 * sx2 

957# XXX if sx3 < 0: 

958# XXX raise ValueError(_ambiguous_) 

959 x3 = acos1(cr12 * sx3 - cx2 * cx1) 

960 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3)) 

961 

962 a, b = _destination2(a1, b1, r13, t13) 

963 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

965 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)): 

966 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple 

967 

968 else: # end point(s) or bearing(s) 

969 _N_vector_ = _MODS.nvectorBase._N_vector_ 

970 

971 x1, d1 = _int3d2(s1, end1, wrap, _1_, _N_vector_, hs) 

972 x2, d2 = _int3d2(s2, end2, wrap, _2_, _N_vector_, hs) 

973 x = x1.cross(x2) 

974 if x.length < EPS: # [nearly] colinear or parallel lines 

975 raise IntersectionError(_colinear_) 

976 a, b = x.philam 

977 # choose intersection similar to sphericalNvector 

978 if not (_intdot(d1, a1, b1, a, b, wrap) * 

979 _intdot(d2, a2, b2, a, b, wrap)) > 0: 

980 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple 

981 

982 h = fmean(hs) if height is None else Height(height) 

983 return _LL3Tuple(degrees90(a), degrees180(b), h, 

984 intersection, LatLon, LatLon_kwds) 

985 

986 

987def intersection(start1, end1, start2, end2, height=None, wrap=False, 

988 **LatLon_and_kwds): 

989 '''Compute the intersection point of two lines, each defined by 

990 two points or by a start point and a bearing from North. 

991 

992 @arg start1: Start point of the first line (L{LatLon}). 

993 @arg end1: End point of the first line (L{LatLon}) or the bearing 

994 at the first start point (compass C{degrees360}). 

995 @arg start2: Start point of the second line (L{LatLon}). 

996 @arg end2: End point of the second line (L{LatLon}) or the bearing 

997 at the second start point (compass C{degrees360}). 

998 @kwarg height: Optional height for the intersection point, 

999 overriding the mean height (C{meter}). 

1000 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{start2}} 

1001 and both B{C{end*}} points (C{bool}). 

1002 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=}L{LatLon} to use 

1003 for the intersection point and optionally additional 

1004 B{C{LatLon}} keyword arguments, ignored if C{B{LatLon} 

1005 is None}. 

1006 

1007 @return: The intersection point as a (B{C{LatLon}}) or if C{B{LatLon} 

1008 is None} a L{LatLon3Tuple}C{(lat, lon, height)}. An alternate 

1009 intersection point might be the L{antipode} to the returned result. 

1010 

1011 @raise IntersectionError: Ambiguous or infinite intersection or colinear, 

1012 parallel or otherwise non-intersecting lines. 

1013 

1014 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}} or B{C{end2}} 

1015 point not L{LatLon}. 

1016 

1017 @raise ValueError: Invalid B{C{height}} or C{null} line. 

1018 ''' 

1019 s1 = _T00.others(start1=start1) 

1020 s2 = _T00.others(start2=start2) 

1021 try: 

1022 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap, **LatLon_and_kwds) 

1023 except (TypeError, ValueError) as x: 

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

1025 

1026 

1027def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0, 

1028 height=None, wrap=False, # was=True 

1029 **LatLon_and_kwds): 

1030 '''Compute the intersection points of two circles each defined by a 

1031 center point and a radius. 

1032 

1033 @arg center1: Center of the first circle (L{LatLon}). 

1034 @arg rad1: Radius of the first circle (C{meter} or C{radians}, see 

1035 B{C{radius}}). 

1036 @arg center2: Center of the second circle (L{LatLon}). 

1037 @arg rad2: Radius of the second circle (C{meter} or C{radians}, see 

1038 B{C{radius}}). 

1039 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}}, 

1040 B{C{rad2}} and B{C{eps}} are given in C{radians}). 

1041 @kwarg eps: Required overlap (C{meter} or C{radians}, see B{C{radius}}). 

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

1043 conventionally) or C{None} for the I{"radical height"} 

1044 at the I{radical line} between both centers. 

1045 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{center2}} 

1046 (C{bool}). 

1047 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=}L{LatLon} to use for 

1048 the intersection points and optionally additional B{C{LatLon}} 

1049 keyword arguments, ignored if C{B{LatLon} is None}. 

1050 

1051 @return: 2-Tuple of the intersection points, each a B{C{LatLon}} 

1052 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, 

1053 lon, height)}. For abutting circles, both intersection 

1054 points are the same instance, aka the I{radical center}. 

1055 

1056 @raise IntersectionError: Concentric, antipodal, invalid or 

1057 non-intersecting circles. 

1058 

1059 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}. 

1060 

1061 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}}, 

1062 B{C{eps}} or B{C{height}}. 

1063 

1064 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}. 

1065 

1066 @see: This U{Answer<https://StackOverflow.com/questions/53324667/ 

1067 find-intersection-coordinates-of-two-circles-on-earth/53331953>}. 

1068 ''' 

1069 c1 = _T00.others(center1=center1) 

1070 c2 = _T00.others(center2=center2) 

1071 try: 

1072 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps, 

1073 height=height, wrap=wrap, 

1074 **LatLon_and_kwds) 

1075 except (TypeError, ValueError) as x: 

1076 raise _xError(x, center1=center1, rad1=rad1, 

1077 center2=center2, rad2=rad2, wrap=wrap) 

1078 

1079 

1080def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2 

1081 height=None, too_d=None, wrap=False, # was=True 

1082 LatLon=LatLon, **LatLon_kwds): 

1083 # (INTERNAL) Intersect two spherical circles, see L{intersections2} 

1084 # above, separated to allow callers to embellish any exceptions 

1085 

1086 def _dest3(bearing, h): 

1087 a, b = _destination2(a1, b1, r1, bearing) 

1088 return _LL3Tuple(degrees90(a), degrees180(b), h, 

1089 intersections2, LatLon, LatLon_kwds) 

1090 

1091 a1, b1 = c1.philam 

1092 a2, b2 = c2.philam 

1093 if wrap: 

1094 a2, b2 = _Wrap.philam(a2, b2) 

1095 

1096 r1, r2, f = _rads3(rad1, rad2, radius) 

1097 if f: # swapped radii, swap centers 

1098 a1, a2 = a2, a1 # PYCHOK swap! 

1099 b1, b2 = b2, b1 # PYCHOK swap! 

1100 

1101 db, b2 = unrollPI(b1, b2, wrap=wrap) 

1102 d = vincentys_(a2, a1, db) # radians 

1103 if d < max(r1 - r2, EPS): 

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

1105 

1106 r = eps if radius is None else (m2radians( 

1107 eps, radius=radius) if eps else _0_0) 

1108 if r < _0_0: 

1109 raise _ValueError(eps=r) 

1110 

1111 x = fsumf_(r1, r2, -d) # overlap 

1112 if x > max(r, EPS): 

1113 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2) 

1114 x = sd * sr1 

1115 if isnear0(x): 

1116 raise _ValueError(_invalid_) 

1117 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI 

1118 

1119 elif x < r: # PYCHOK no cover 

1120 t = (d * radius) if too_d is None else too_d 

1121 raise IntersectionError(_too_(_Fmt.distant(t))) 

1122 

1123 if height is None: # "radical height" 

1124 f = _radical2(d, r1, r2).ratio 

1125 h = Height(favg(c1.height, c2.height, f=f)) 

1126 else: 

1127 h = Height(height) 

1128 

1129 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap) 

1130 if x < EPS4: # externally ... 

1131 r = _dest3(b, h) 

1132 elif x > _PI_EPS4: # internally ... 

1133 r = _dest3(b + PI, h) 

1134 else: 

1135 return _dest3(b + x, h), _dest3(b - x, h) 

1136 return r, r # ... abutting circles 

1137 

1138 

1139@deprecated_function 

1140def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover 

1141 '''DEPRECATED, use function L{pygeodesy.ispolar}. 

1142 ''' 

1143 return ispolar(points, wrap=wrap) 

1144 

1145 

1146def _LL3Tuple(lat, lon, height, where, LatLon, LatLon_kwds): 

1147 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}. 

1148 ''' 

1149 n = where.__name__ 

1150 if LatLon is None: 

1151 r = LatLon3Tuple(lat, lon, height, name=n) 

1152 else: 

1153 kwds = _xkwds(LatLon_kwds, height=height, name=n) 

1154 r = LatLon(lat, lon, **kwds) 

1155 return r 

1156 

1157 

1158def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds): 

1159 '''Compute the I{geographic} mean of several points. 

1160 

1161 @arg points: Points to be averaged (L{LatLon}[]). 

1162 @kwarg height: Optional height at mean point, overriding the mean height 

1163 (C{meter}). 

1164 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}} (C{bool}). 

1165 @kwarg LatLon: Optional class to return the mean point (L{LatLon}) or C{None}. 

1166 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments, 

1167 ignored if C{B{LatLon} is None}. 

1168 

1169 @return: The geographic mean and height (B{C{LatLon}}) or if C{B{LatLon} 

1170 is None}, a L{LatLon3Tuple}C{(lat, lon, height)}. 

1171 

1172 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

1173 

1174 @raise ValueError: No B{C{points}} or invalid B{C{height}}. 

1175 ''' 

1176 def _N_vs(ps, w): 

1177 Ps = _T00.PointsIter(ps, wrap=w) 

1178 for p in Ps.iterate(closed=False): 

1179 yield p._N_vector 

1180 

1181 m = _MODS.nvectorBase 

1182 # geographic, vectorial mean 

1183 n = m.sumOf(_N_vs(points, wrap), h=height, Vector=m.NvectorBase) 

1184 lat, lon, h = n.latlonheight 

1185 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds) 

1186 

1187 

1188@deprecated_function 

1189def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover 

1190 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}. 

1191 

1192 @return: ... 2-tuple C{(closest, distance)} of the C{closest} 

1193 point (L{LatLon}) on the polygon and the C{distance} 

1194 between the C{closest} and the given B{C{point}}. The 

1195 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat, 

1196 lon)} if C{B{LatLon} is None} ... 

1197 ''' 

1198 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple 

1199 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None: 

1200 ll = LatLon2Tuple(ll.lat, ll.lon) 

1201 return ll, d 

1202 

1203 

1204def nearestOn3(point, points, closed=False, radius=R_M, wrap=False, adjust=True, 

1205 limit=9, **LatLon_and_kwds): 

1206 '''Locate the point on a path or polygon closest to a reference point. 

1207 

1208 Distances are I{approximated} using function L{equirectangular4 

1209 <pygeodesy.equirectangular4>}, subject to the supplied B{C{options}}. 

1210 

1211 @arg point: The reference point (L{LatLon}). 

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

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

1214 @kwarg radius: Mean earth radius (C{meter}). 

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

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

1217 @kwarg adjust: See function L{equirectangular4<pygeodesy.equirectangular4>} 

1218 (C{bool}). 

1219 @kwarg limit: See function L{equirectangular4<pygeodesy.equirectangular4>} 

1220 (C{degrees}), default C{9 degrees} is about C{1,000 Km} (for 

1221 (mean spherical earth radius L{R_KM}). 

1222 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=L{LatLon}} to return the 

1223 closest point and optionally additional C{B{LatLon}} keyword 

1224 arguments or specify C{B{LatLon}=None}. 

1225 

1226 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the 

1227 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat, 

1228 lon, height)} if C{B{LatLon} is None}. The C{distance} is 

1229 the L{equirectangular4<pygeodesy.equirectangular4>} distance 

1230 between the C{closest} and the given B{C{point}} converted to 

1231 C{meter}, same units as B{C{radius}}. The C{angle} from the 

1232 given B{C{point}} to the C{closest} is in compass C{degrees360}, 

1233 like function L{compassAngle<pygeodesy.compassAngle>}. The 

1234 C{height} is the (interpolated) height at the C{closest} point. 

1235 

1236 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}}, 

1237 see function L{equirectangular4<pygeodesy.equirectangular4>}. 

1238 

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

1240 

1241 @raise TypeError: Some B{C{points}} are not C{LatLon}. 

1242 

1243 @raise ValueError: Invalid B{C{radius}}. 

1244 

1245 @see: Functions L{equirectangular4<pygeodesy.equirectangular4>} and 

1246 L{nearestOn5<pygeodesy.nearestOn5>}. 

1247 ''' 

1248 t = _nearestOn5(point, points, closed=closed, wrap=wrap, 

1249 adjust=adjust, limit=limit) 

1250 d = degrees2m(t.distance, radius=radius) 

1251 h = t.height 

1252 n = nearestOn3.__name__ 

1253 

1254 LL, kwds = _xkwds_pop2(LatLon_and_kwds, LatLon=LatLon) 

1255 r = LatLon3Tuple(t.lat, t.lon, h, name=n) if LL is None else \ 

1256 LL(t.lat, t.lon, **_xkwds(kwds, height=h, name=n)) 

1257 return NearestOn3Tuple(r, d, t.angle, name=n) 

1258 

1259 

1260def perimeterOf(points, closed=False, radius=R_M, wrap=True): 

1261 '''Compute the perimeter of a (spherical) polygon or composite 

1262 (with great circle arcs joining the points). 

1263 

1264 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP} 

1265 or L{BooleanGH}). 

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

1267 @kwarg radius: Mean earth radius (C{meter}) or C{None}. 

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

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

1270 

1271 @return: Polygon perimeter (C{meter}, same units as B{C{radius}} 

1272 or C{radians} if C{B{radius} is None}). 

1273 

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

1275 

1276 @raise TypeError: Some B{C{points}} are not L{LatLon}. 

1277 

1278 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with 

1279 C{B{points}} a composite. 

1280 

1281 @note: Distances are based on function L{vincentys_<pygeodesy.vincentys_>}. 

1282 

1283 @see: Functions L{perimeterOf<pygeodesy.perimeterOf>}, 

1284 L{sphericalNvector.perimeterOf} and L{ellipsoidalKarney.perimeterOf}. 

1285 ''' 

1286 def _rads(ps, c, w): # angular edge lengths in radians 

1287 Ps = _T00.PointsIter(ps, loop=1, wrap=w) 

1288 a1, b1 = Ps[0].philam 

1289 for p in Ps.iterate(closed=c): 

1290 a2, b2 = p.philam 

1291 db, b2 = unrollPI(b1, b2, wrap=w and not (c and Ps.looped)) 

1292 yield vincentys_(a2, a1, db) 

1293 a1, b1 = a2, b2 

1294 

1295 if _MODS.booleans.isBoolean(points): 

1296 if not closed: 

1297 raise _ValueError(closed=closed, points=_composite_) 

1298 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap) 

1299 else: 

1300 r = fsum(_rads(points, closed, wrap)) 

1301 return _radians2m(r, radius) 

1302 

1303 

1304def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M, 

1305 excess=excessAbc_, 

1306 wrap=False): 

1307 '''Compute the angles, sides, and area of a (spherical) triangle. 

1308 

1309 @arg latA: First corner latitude (C{degrees}). 

1310 @arg lonA: First corner longitude (C{degrees}). 

1311 @arg latB: Second corner latitude (C{degrees}). 

1312 @arg lonB: Second corner longitude (C{degrees}). 

1313 @arg latC: Third corner latitude (C{degrees}). 

1314 @arg lonC: Third corner longitude (C{degrees}). 

1315 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, 

1316 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) 

1317 or C{None}. 

1318 @kwarg excess: I{Spherical excess} callable (L{excessAbc_}, 

1319 L{excessGirard_} or L{excessLHuilier_}). 

1320 @kwarg wrap: If C{True}, wrap and L{unroll180<pygeodesy.unroll180>} 

1321 longitudes (C{bool}). 

1322 

1323 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with 

1324 spherical angles C{A}, C{B} and C{C}, angular sides 

1325 C{a}, C{b} and C{c} all in C{degrees} and C{area} 

1326 in I{square} C{meter} or same units as B{C{radius}} 

1327 I{squared} or if C{B{radius}=0} or C{None}, a 

1328 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with 

1329 I{spherical deficit} C{D} and I{spherical excess} 

1330 C{E} as the C{unit area}, all in C{radians}. 

1331 ''' 

1332 t = triangle8_(Phid(latA=latA), Lamd(lonA=lonA), 

1333 Phid(latB=latB), Lamd(lonB=lonB), 

1334 Phid(latC=latC), Lamd(lonC=lonC), 

1335 excess=excess, wrap=wrap) 

1336 return _t7Tuple(t, radius) 

1337 

1338 

1339def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_, 

1340 wrap=False): 

1341 '''Compute the angles, sides, I{spherical deficit} and I{spherical 

1342 excess} of a (spherical) triangle. 

1343 

1344 @arg phiA: First corner latitude (C{radians}). 

1345 @arg lamA: First corner longitude (C{radians}). 

1346 @arg phiB: Second corner latitude (C{radians}). 

1347 @arg lamB: Second corner longitude (C{radians}). 

1348 @arg phiC: Third corner latitude (C{radians}). 

1349 @arg lamC: Third corner longitude (C{radians}). 

1350 @kwarg excess: I{Spherical excess} callable (L{excessAbc_}, 

1351 L{excessGirard_} or L{excessLHuilier_}). 

1352 @kwarg wrap: If C{True}, L{unrollPI<pygeodesy.unrollPI>} the 

1353 longitudinal deltas (C{bool}). 

1354 

1355 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with 

1356 spherical angles C{A}, C{B} and C{C}, angular sides 

1357 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and 

1358 I{spherical excess} C{E}, all in C{radians}. 

1359 ''' 

1360 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC): 

1361 d, _ = unrollPI(lamB, lamC, wrap=w) 

1362 a = vincentys_(phiC, phiB, d) 

1363 return a, (phiB, lamB, phiC, lamC, phiA, lamA) # rotate A, B, C 

1364 

1365 def _A_r(a, sa, ca, sb, cb, sc, cc): 

1366 s = sb * sc 

1367 A = acos1((ca - cb * cc) / s) if isnon0(s) else a 

1368 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2_'s 

1369 

1370 # notation: side C{a} is oposite to corner C{A}, etc. 

1371 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC) 

1372 b, r = _a_r(wrap, *r) 

1373 c, _ = _a_r(wrap, *r) 

1374 

1375 A, r = _A_r(a, *sincos2_(a, b, c)) 

1376 B, r = _A_r(b, *r) 

1377 C, _ = _A_r(c, *r) 

1378 

1379 D = fsumf_(PI2, -a, -b, -c) # deficit aka defect 

1380 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else ( 

1381 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else 

1382 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b)))) 

1383 

1384 return Triangle8Tuple(A, a, B, b, C, c, D, E) 

1385 

1386 

1387def _t7Tuple(t, radius): 

1388 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}. 

1389 ''' 

1390 if radius: # not in (None, _0_0) 

1391 r = radius if _isRadius(radius) else \ 

1392 _ellipsoidal_datum(radius).ellipsoid.Rmean 

1393 A, B, C = map1(degrees, t.A, t.B, t.C) 

1394 t = Triangle7Tuple(A, (r * t.a), 

1395 B, (r * t.b), 

1396 C, (r * t.c), t.E * r**2) 

1397 return t 

1398 

1399 

1400__all__ += _ALL_OTHER(Cartesian, LatLon, # classes 

1401 areaOf, # functions 

1402 intersecant2, intersection, intersections2, ispolar, 

1403 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1404 meanOf, 

1405 nearestOn2, nearestOn3, 

1406 perimeterOf, 

1407 sumOf, # XXX == vector3d.sumOf 

1408 triangle7, triangle8_) 

1409 

1410# **) MIT License 

1411# 

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

1413# 

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

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

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

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

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

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

1420# 

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

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

1423# 

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

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

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

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

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

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

1430# OTHER DEALINGS IN THE SOFTWARE.