Coverage for pygeodesy/sphericalTrigonometry.py: 93%

386 statements  

« prev     ^ index     » next       coverage.py v7.6.1, created at 2024-10-22 18:16 -0400

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, 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.10.12' 

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 

186 the great circle defined by a start and 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 the given initial bearing. 

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 on the given initial bearing from this point. 

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 return Vector(sb * ct - cb * sa * st, 

275 -cb * ct - sb * sa * st, 

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

277 

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

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

280 to an other point. 

281 

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

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

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

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

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

287 C{sphericalNvector.LatLon.initialBearingTo}. 

288 

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

290 

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

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

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

294 

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

296 ''' 

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

298 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

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

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

301 

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

303 

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

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

306 and an other point. 

307 

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

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

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

311 @kwarg height: Optional height, overriding the intermediate 

312 height (C{meter}). 

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

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

315 

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

317 

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

319 

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

321 

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

323 ''' 

324 p = self 

325 f = Scalar(fraction=fraction) 

326 if not isnear0(f): 

327 p = p.others(other) 

328 if wrap: 

329 p = _Wrap.point(p) 

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

331 a1, b1 = self.philam 

332 a2, b2 = p.philam 

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

334 r = vincentys_(a2, a1, db) 

335 sr = sin(r) 

336 if isnon0(sr): 

337 sa1, ca1, sa2, ca2, \ 

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

339 

340 t = f * r 

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

342 b = sin( t) # / sr superflous 

343 

344 x = a * ca1 * cb1 + b * ca2 * cb2 

345 y = a * ca1 * sb1 + b * ca2 * sb2 

346 z = a * sa1 + b * sa2 

347 

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

349 b = atan2d(y, x) 

350 

351 else: # PYCHOK no cover 

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

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

354 

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

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

357 return p 

358 

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

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

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

362 

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

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

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

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

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

368 C{degrees360}). 

369 @kwarg height: Optional height for intersection point, 

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

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

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

373 

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

375 intersection point might be the L{antipode} to 

376 the returned result. 

377 

378 @raise IntersectionError: Ambiguous or infinite intersection 

379 or colinear, parallel or otherwise 

380 non-intersecting lines. 

381 

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

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

384 

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

386 ''' 

387 try: 

388 s2 = self.others(other) 

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

390 LatLon=self.classof) 

391 except (TypeError, ValueError) as x: 

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

393 other=other, end2=end2, wrap=wrap) 

394 

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

396 height=None, wrap=True): 

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

398 by a center point and radius. 

399 

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

401 see B{C{radius}}). 

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

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

404 see B{C{radius}}). 

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

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

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

408 B{C{radius}}). 

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

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

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

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

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

414 

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

416 instance. For abutting circles, both intersection 

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

418 

419 @raise IntersectionError: Concentric, antipodal, invalid or 

420 non-intersecting circles. 

421 

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

423 

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

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

426 ''' 

427 try: 

428 c2 = self.others(other) 

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

430 height=height, wrap=wrap, 

431 LatLon=self.classof) 

432 except (TypeError, ValueError) as x: 

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

434 other=other, rad2=rad2, wrap=wrap) 

435 

436 @deprecated_method 

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

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

439 return self.isenclosedBy(points) 

440 

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

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

443 

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

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

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

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

448 

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

450 composite, C{False} otherwise. 

451 

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

453 

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

455 

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

457 

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

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

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

461 ''' 

462 if _MODS.booleans.isBoolean(points): 

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

464 

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

466 n0 = self._N_vector 

467 

468 v2 = Ps[0]._N_vector 

469 p1 = Ps[1] 

470 v1 = p1._N_vector 

471 # check whether this point on same side of all 

472 # polygon edges (to the left or right depending 

473 # on the anti-/clockwise polygon direction) 

474 gc1 = v2.cross(v1) 

475 t0 = gc1.angleTo(n0) > PI_2 

476 s0 = None 

477 # get great-circle vector for each edge 

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

479 if wrap and not Ps.looped: 

480 p2 = _unrollon(p1, p2) 

481 p1 = p2 

482 v2 = p2._N_vector 

483 gc = v1.cross(v2) 

484 t = gc.angleTo(n0) > PI_2 

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

486 return False # outside 

487 

488 # check for convex polygon: angle between 

489 # gc vectors, signed by direction of n0 

490 # (otherwise the test above is not reliable) 

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

492 if s != s0: 

493 if s0 is None: 

494 s0 = s 

495 else: 

496 t = _Fmt.SQUARE(points=i) 

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

498 gc1, v1 = gc, v2 

499 

500 return True # inside 

501 

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

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

504 

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

506 @kwarg height: Optional height for midpoint, overriding 

507 the mean height (C{meter}). 

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

509 may be negative or greater than 1.0. 

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

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

512 

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

514 

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

516 

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

518 

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

520 ''' 

521 if fraction is _0_5: 

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

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

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

525 

526 x = ca2 * cdb + ca1 

527 y = ca2 * sdb 

528 

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

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

531 

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

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

534 else: 

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

536 return r 

537 

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

539 '''Locate the point between two points closest to this point. 

540 

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

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

543 

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

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

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

547 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

548 L{sphericalTrigonometry.nearestOn3} and 

549 L{pygeodesy.equirectangular4}, 

550 

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

552 

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

554 see function L{pygeodesy.equirectangular4}. 

555 

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

557 is not (yet) supported. 

558 

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

560 

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

562 

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

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

565 ''' 

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

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

568 if not w: 

569 self._notImplemented(within=w) 

570 

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

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

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

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

575# return point1 

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

577# if isnear0(d): 

578# return point1 # or point2 

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

580# return point2 

581# f = a / d 

582# if within: 

583# if f > EPS1: 

584# return point2 

585# elif f < EPS: 

586# return point1 

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

588 

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

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

591 **kwds)[0] 

592 

593 @deprecated_method 

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

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

596 

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

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

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

600 units of B{C{radius}}. 

601 ''' 

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

603 return r.closest, r.distance 

604 

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

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

607 

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

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

610 

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

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

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

614 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

615 L{sphericalTrigonometry.nearestOn3} and 

616 L{pygeodesy.equirectangular4}, 

617 

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

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

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

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

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

623 function L{pygeodesy.compassAngle}. 

624 

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

626 see function L{pygeodesy.equirectangular4}. 

627 

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

629 

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

631 

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

633 

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

635 and L{pygeodesy.nearestOn5}. 

636 ''' 

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

638 LatLon=self.classof, **wrap_adjust_limit) 

639 

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

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

642 

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

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

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

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

647 

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

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

650 and C{M} if available. 

651 

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

653 ''' 

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

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

656 

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

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

659 

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

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

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

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

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

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

666 

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

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

669 

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

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

672 ''' 

673 B = self.others(otherB=otherB) 

674 C = self.others(otherC=otherC) 

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

676 

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

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

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

680 

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

682 '''Locate a point given this, an other point and the (initial) bearing 

683 at this and at the other point. 

684 

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

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

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

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

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

690 

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

692 

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

694 ''' 

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

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

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

698 

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

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

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

702 of three corresponding circles. 

703 

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

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

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

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

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

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

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

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

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

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

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

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

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

717 

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

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

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

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

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

723 

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

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

726 

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

728 largest I{radial} overlap found. 

729 

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

731 respectively farthest intersection margin. 

732 

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

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

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

736 

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

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

739 no intersection or all (near-)concentric if 

740 C{B{area}=False}. 

741 

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

743 

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

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

746 or B{C{radius}}. 

747 ''' 

748 return _trilaterate5(self, distance1, 

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

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

751 area=area, radius=radius, eps=eps, wrap=wrap) 

752 

753 

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

755 

756 

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

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

759 (with the pointsjoined by great circle arcs). 

760 

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

762 or L{BooleanGH}). 

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

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

765 or C{None}. 

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

767 (C{bool}). 

768 

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

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

771 

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

773 

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

775 

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

777 

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

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

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

781 

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

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

784 L{ellipsoidalKarney.areaOf}. 

785 ''' 

786 if _MODS.booleans.isBoolean(points): 

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

788 

789 _at2, _t_2 = atan2, tan_2 

790 _un, _w180 = unrollPI, wrap180 

791 

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

793 p1 = p2 = Ps[0] 

794 a1, b1 = p1.philam 

795 ta1, z1 = _t_2(a1), None 

796 

797 A = Fsum() # mean phi 

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

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

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

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

802 D = Fsum() 

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

804 a2, b2 = p2.philam 

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

806 A += a2 

807 ta2 = _t_2(a2) 

808 tdb = _t_2(db, points=i) 

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

810 _1_0 + ta1 * ta2) 

811 ta1, b1 = ta2, b2 

812 

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

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

815 if z1 is not None: 

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

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

818 p1, z1 = p2, z2 

819 

820 R = abs(R * _2_0) 

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

822 R = abs(R - PI2) 

823 if radius: 

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

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

826 return float(R) 

827 

828 

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

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

831 

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

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

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

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

836 

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

838 ''' 

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

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

841 ca *= sr 

842 

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

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

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

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

847 

848 

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

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

851 # and similar logic in .ellipsoidalBaseDI._intersect3 

852 a1, b1 = s.philam 

853 

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

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

856 else: # must be a point 

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

858 hs.append(end.height) 

859 a2, b2 = end.philam 

860 if wrap: 

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

862 

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

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

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

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

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

868 -(b1 + b2) * _0_5) 

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

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

871 x = Vector(sb12 * cb21 - cb12 * sb21, 

872 cb12 * cb21 + sb12 * sb21, 

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

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

875 

876 

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

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

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

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

881 

882 

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

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

885 two points or as a point and bearing. 

886 

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

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

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

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

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

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

893 @kwarg radius_exact_height_wrap: Optional keyword arguments, see method 

894 L{intersecant2<pygeodesy.sphericalBase.LatLonSphericalBase. 

895 intersecant2>} for further details. 

896 

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

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

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

900 

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

902 

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

904 not L{LatLon}. 

905 

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

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

908 ''' 

909 c = _T00.others(center=center) 

910 p = _T00.others(point=point) 

911 try: 

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

913 except (TypeError, ValueError) as x: 

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

915 **radius_exact_height_wrap) 

916 

917 

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

919 LatLon=None, **LatLon_kwds): 

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

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

922 

923 s1, s2 = start1, start2 

924 if wrap: 

925 s2 = _Wrap.point(s2) 

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

927 

928 a1, b1 = s1.philam 

929 a2, b2 = s2.philam 

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

931 r12 = vincentys_(a2, a1, db) 

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

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

934 

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

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

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

938 

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

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

941 raise IntersectionError(_parallel_) 

942 # handle domain error for equivalent longitudes, 

943 # see also functions asin_safe and acos_safe at 

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

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

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

947 if sin(db) > 0: 

948 t21 = PI2 - t21 

949 else: 

950 t12 = PI2 - t12 

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

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

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

954 raise IntersectionError(_infinite_) 

955 sx3 = sx1 * sx2 

956# XXX if sx3 < 0: 

957# XXX raise ValueError(_ambiguous_) 

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

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

960 

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

962 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

966 

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

968 _N_vector_ = _MODS.nvectorBase._N_vector_ 

969 

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

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

972 x = x1.cross(x2) 

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

974 raise IntersectionError(_colinear_) 

975 a, b = x.philam 

976 # choose intersection similar to sphericalNvector 

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

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

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

980 

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

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

983 intersection, LatLon, LatLon_kwds) 

984 

985 

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

987 LatLon=LatLon, **LatLon_kwds): 

988 '''Compute the intersection point of two lines, each defined 

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

990 

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

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

993 the initial bearing at the first start point 

994 (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 

997 the initial bearing at the second start point 

998 (compass C{degrees360}). 

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

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

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

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

1003 @kwarg LatLon: Optional class to return the intersection 

1004 point (L{LatLon}) or C{None}. 

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

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

1007 

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

1009 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon, 

1010 height)}. An alternate intersection point might 

1011 be the L{antipode} to the returned result. 

1012 

1013 @raise IntersectionError: Ambiguous or infinite intersection 

1014 or colinear, parallel or otherwise 

1015 non-intersecting lines. 

1016 

1017 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}} 

1018 or B{C{end2}} point not L{LatLon}. 

1019 

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

1021 ''' 

1022 s1 = _T00.others(start1=start1) 

1023 s2 = _T00.others(start2=start2) 

1024 try: 

1025 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap, 

1026 LatLon=LatLon, **LatLon_kwds) 

1027 except (TypeError, ValueError) as x: 

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

1029 

1030 

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

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

1033 LatLon=LatLon, **LatLon_kwds): 

1034 '''Compute the intersection points of two circles each defined 

1035 by a center point and a radius. 

1036 

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

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

1039 see B{C{radius}}). 

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

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

1042 see B{C{radius}}). 

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

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

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

1046 B{C{radius}}). 

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

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

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

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

1051 (C{bool}). 

1052 @kwarg LatLon: Optional class to return the intersection 

1053 points (L{LatLon}) or C{None}. 

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

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

1056 

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

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

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

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

1061 

1062 @raise IntersectionError: Concentric, antipodal, invalid or 

1063 non-intersecting circles. 

1064 

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

1066 

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

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

1069 

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

1071 

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

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

1074 ''' 

1075 c1 = _T00.others(center1=center1) 

1076 c2 = _T00.others(center2=center2) 

1077 try: 

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

1079 height=height, wrap=wrap, 

1080 LatLon=LatLon, **LatLon_kwds) 

1081 except (TypeError, ValueError) as x: 

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

1083 center2=center2, rad2=rad2, wrap=wrap) 

1084 

1085 

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

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

1088 LatLon=LatLon, **LatLon_kwds): 

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

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

1091 

1092 def _dest3(bearing, h): 

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

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

1095 intersections2, LatLon, LatLon_kwds) 

1096 

1097 a1, b1 = c1.philam 

1098 a2, b2 = c2.philam 

1099 if wrap: 

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

1101 

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

1103 if f: # swapped radii, swap centers 

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

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

1106 

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

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

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

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

1111 

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

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

1114 if r < _0_0: 

1115 raise _ValueError(eps=r) 

1116 

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

1118 if x > max(r, EPS): 

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

1120 x = sd * sr1 

1121 if isnear0(x): 

1122 raise _ValueError(_invalid_) 

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

1124 

1125 elif x < r: # PYCHOK no cover 

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

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

1128 

1129 if height is None: # "radical height" 

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

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

1132 else: 

1133 h = Height(height) 

1134 

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

1136 if x < EPS4: # externally ... 

1137 r = _dest3(b, h) 

1138 elif x > _PI_EPS4: # internally ... 

1139 r = _dest3(b + PI, h) 

1140 else: 

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

1142 return r, r # ... abutting circles 

1143 

1144 

1145@deprecated_function 

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

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

1148 ''' 

1149 return ispolar(points, wrap=wrap) 

1150 

1151 

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

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

1154 ''' 

1155 n = where.__name__ 

1156 if LatLon is None: 

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

1158 else: 

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

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

1161 return r 

1162 

1163 

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

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

1166 

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

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

1169 (C{meter}). 

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

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

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

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

1174 

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

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

1177 

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

1179 

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

1181 ''' 

1182 def _N_vs(ps, w): 

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

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

1185 yield p._N_vector 

1186 

1187 m = _MODS.nvectorBase 

1188 # geographic, vectorial mean 

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

1190 lat, lon, h = n.latlonheight 

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

1192 

1193 

1194@deprecated_function 

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

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

1197 

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

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

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

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

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

1203 ''' 

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

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

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

1207 return ll, d 

1208 

1209 

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

1211 limit=9, **LatLon_and_kwds): 

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

1213 

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

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

1216 

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

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

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

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

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

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

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

1224 (C{bool}). 

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

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

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

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

1229 closest point and optional, additional C{B{LatLon}} keyword 

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

1231 

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

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

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

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

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

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

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

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

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

1241 

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

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

1244 

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

1246 

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

1248 

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

1250 

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

1252 L{nearestOn5<pygeodesy.nearestOn5>}. 

1253 ''' 

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

1255 adjust=adjust, limit=limit) 

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

1257 h = t.height 

1258 n = nearestOn3.__name__ 

1259 

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

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

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

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

1264 

1265 

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

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

1268 (with great circle arcs joining the points). 

1269 

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

1271 or L{BooleanGH}). 

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

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

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

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

1276 

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

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

1279 

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

1281 

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

1283 

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

1285 C{B{points}} a composite. 

1286 

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

1288 

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

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

1291 ''' 

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

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

1294 a1, b1 = Ps[0].philam 

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

1296 a2, b2 = p.philam 

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

1298 yield vincentys_(a2, a1, db) 

1299 a1, b1 = a2, b2 

1300 

1301 if _MODS.booleans.isBoolean(points): 

1302 if not closed: 

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

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

1305 else: 

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

1307 return _radians2m(r, radius) 

1308 

1309 

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

1311 excess=excessAbc_, 

1312 wrap=False): 

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

1314 

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

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

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

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

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

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

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

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

1323 or C{None}. 

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

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

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

1327 longitudes (C{bool}). 

1328 

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

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

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

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

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

1334 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in 

1335 C{radians} with the I{spherical excess} C{E} as the 

1336 C{unit area} in C{radians}. 

1337 ''' 

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

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

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

1341 excess=excess, wrap=wrap) 

1342 return _t7Tuple(t, radius) 

1343 

1344 

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

1346 wrap=False): 

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

1348 excess} of a (spherical) triangle. 

1349 

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

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

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

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

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

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

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

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

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

1359 longitudinal deltas (C{bool}). 

1360 

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

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

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

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

1365 ''' 

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

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

1368 a = vincentys_(phiC, phiB, d) 

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

1370 

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

1372 s = sb * sc 

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

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

1375 

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

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

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

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

1380 

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

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

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

1384 

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

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

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

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

1389 

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

1391 

1392 

1393def _t7Tuple(t, radius): 

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

1395 ''' 

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

1397 r = radius if _isRadius(radius) else \ 

1398 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1401 B, (r * t.b), 

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

1403 return t 

1404 

1405 

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

1407 areaOf, # functions 

1408 intersecant2, intersection, intersections2, ispolar, 

1409 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1410 meanOf, 

1411 nearestOn2, nearestOn3, 

1412 perimeterOf, 

1413 sumOf, # XXX == vector3d.sumOf 

1414 triangle7, triangle8_) 

1415 

1416# **) MIT License 

1417# 

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

1419# 

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

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

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

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

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

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

1426# 

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

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

1429# 

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

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

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

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

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

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

1436# OTHER DEALINGS IN THE SOFTWARE.