Coverage for pygeodesy/sphericalTrigonometry.py: 95%

376 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-05 15:46 -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-2016} 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, isscalar, 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 _ValueError, IntersectionError, _xError, \ 

26 _xkwds, _xkwds_get, _xkwds_pop 

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

28from pygeodesy.fsums import Fsum, fsum, fsum_ 

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_, _LatLon_, _near_, _not_, _null_, \ 

35 _points_, _SPACE_, _too_ 

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

37# from pygeodesy.named import notImplemented # from .points 

38from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \ 

39 NearestOn3Tuple, Triangle7Tuple, \ 

40 Triangle8Tuple 

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

42 notImplemented, Fmt as _Fmt # XXX shadowed 

43from pygeodesy.props import deprecated_function, deprecated_method 

44from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, _intersecant2, \ 

45 LatLonSphericalBase, _rads3, _trilaterate5 

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

47from pygeodesy.units import Bearing_, Height, Lam_, Phi_, Radius, \ 

48 Radius_, Scalar 

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

50 m2radians, radiansPI2, sincos2_, tan_2, unrollPI, \ 

51 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__ = '23.04.05' 

58 

59_parallel_ = 'parallel' 

60_path_ = 'path' 

61 

62_PI_EPS4 = PI - EPS4 

63if _PI_EPS4 >= PI: 

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

65 

66 

67class Cartesian(CartesianSphericalBase): 

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

69 spherical, geodetic L{LatLon}. 

70 ''' 

71 

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

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

74 

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

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

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

78 

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

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

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

82 

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

84 ''' 

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

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

87 

88 

89class LatLon(LatLonSphericalBase): 

90 '''New point on spherical model earth model. 

91 

92 @example: 

93 

94 >>> p = LatLon(52.205, 0.119) # height=0 

95 ''' 

96 

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

98 '''Compute the (angular) distance (signed) from the start to 

99 the closest point on the great circle path defined by a 

100 start and an end point. 

101 

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

103 great circle path, the along-track distance is the distance 

104 from the start point to the point where the perpendicular 

105 crosses the path. 

106 

107 @arg start: Start point of great circle path (L{LatLon}). 

108 @arg end: End point of great circle path (L{LatLon}). 

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

110 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

111 

112 @return: Distance along the great circle path (C{meter}, 

113 same units as B{C{radius}}) or C{radians} if 

114 C{B{radius} is None}, positive if after the B{C{start}} 

115 toward the B{C{end}} point of the path or negative 

116 if before the B{C{start}} point. 

117 

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

119 

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

121 

122 @example: 

123 

124 >>> p = LatLon(53.2611, -0.7972) 

125 

126 >>> s = LatLon(53.3206, -1.7297) 

127 >>> e = LatLon(53.1887, 0.1334) 

128 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58 

129 ''' 

130 r, x, b = self._trackDistanceTo3(start, end, radius, wrap) 

131 cx = cos(x) 

132 return _0_0 if isnear0(cx) else \ 

133 _r2m(copysign0(acos1(cos(r) / cx), cos(b)), radius) 

134 

135 @deprecated_method 

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

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

138 ''' 

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

140 

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

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

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

144 

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

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

147 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

148 

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

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

151 ''' 

152 self.others(other) 

153 

154 a1, b1 = self.philam 

155 a2, b2 = other.philam 

156 

157 a = radians(lat) 

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

159 

160 sa, ca, sa1, ca1, \ 

161 sa2, ca2, sdb, cdb = sincos2_(a, a1, a2, db) 

162 sa1 *= ca2 * ca 

163 

164 x = sa1 * sdb 

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

166 z = ca1 * sdb * ca2 * sa 

167 

168 h = hypot(x, y) 

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

170 return None # great circle doesn't reach latitude 

171 

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

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

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

175 

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

177 '''Compute the (angular) distance (signed) from this point to 

178 the great circle defined by a start and an end point. 

179 

180 @arg start: Start point of great circle path (L{LatLon}). 

181 @arg end: End point of great circle path (L{LatLon}). 

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

183 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

184 

185 @return: Distance to great circle (negative if to the left 

186 or positive if to the right of the path) (C{meter}, 

187 same units as B{C{radius}} or C{radians} if 

188 B{C{radius}} is C{None}). 

189 

190 @raise TypeError: The B{C{start}} or B{C{end}} point is not L{LatLon}. 

191 

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

193 

194 @example: 

195 

196 >>> p = LatLon(53.2611, -0.7972) 

197 

198 >>> s = LatLon(53.3206, -1.7297) 

199 >>> e = LatLon(53.1887, 0.1334) 

200 >>> d = p.crossTrackDistanceTo(s, e) # -307.5 

201 ''' 

202 _, x, _ = self._trackDistanceTo3(start, end, radius, wrap) 

203 return _r2m(x, radius) 

204 

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

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

207 travelled the given distance on the given initial bearing. 

208 

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

210 B{C{radius}}). 

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

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

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

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

215 

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

217 

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

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

220 

221 @example: 

222 

223 >>> p1 = LatLon(51.4778, -0.0015) 

224 >>> p2 = p1.destination(7794, 300.7) 

225 >>> p2.toStr() # '51.5135°N, 000.0983°W' 

226 ''' 

227 a, b = self.philam 

228 r, t = _angular(distance, radius), Bearing_(bearing) 

229 

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

231 h = self.height if height is None else Height(height) 

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

233 

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

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

236 

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

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

239 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

240 

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

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

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

244 

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

246 

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

248 

249 @example: 

250 

251 >>> p1 = LatLon(52.205, 0.119) 

252 >>> p2 = LatLon(48.857, 2.351); 

253 >>> d = p1.distanceTo(p2) # 404300 

254 ''' 

255 self.others(other) 

256 

257 a1, b1 = self.philam 

258 a2, b2 = other.philam 

259 

260 db, _ = unrollPI(b1, b2, wrap=wrap) 

261 return _r2m(vincentys_(a2, a1, db), radius) 

262 

263# @Property_RO 

264# def Ecef(self): 

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

266# ''' 

267# return _MODS.ecef.EcefKarney 

268 

269 def greatCircle(self, bearing): 

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

271 on the given initial bearing from this point. 

272 

273 Direction of vector is such that initial bearing vector 

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

275 

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

277 

278 @return: Vector representing great circle (L{Vector3d}). 

279 

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

281 

282 @example: 

283 

284 >>> p = LatLon(53.3206, -1.7297) 

285 >>> g = p.greatCircle(96.0) 

286 >>> g.toStr() # (-0.794, 0.129, 0.594) 

287 ''' 

288 a, b = self.philam 

289 

290 t = Bearing_(bearing) 

291 sa, ca, sb, cb, st, ct = sincos2_(a, b, t) 

292 

293 return Vector3d(sb * ct - cb * sa * st, 

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

295 ca * st) # XXX .unit()? 

296 

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

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

299 to an other point. 

300 

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

302 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

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

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

305 C{sphericalNvector.LatLon.initialBearingTo}. 

306 

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

308 

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

310 provided both B{C{raiser}} is C{True} and 

311 L{pygeodesy.crosserrors} is C{True}. 

312 

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

314 

315 @example: 

316 

317 >>> p1 = LatLon(52.205, 0.119) 

318 >>> p2 = LatLon(48.857, 2.351) 

319 >>> b = p1.initialBearingTo(p2) # 156.2 

320 ''' 

321 self.others(other) 

322 

323 a1, b1 = self.philam 

324 a2, b2 = other.philam 

325 

326 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

327 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(b2 - b1)) < EPS: 

328 raise CrossError(_points_, self, txt=_coincident_) 

329 

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

331 

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

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

334 and an other point. 

335 

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

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

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

339 @kwarg height: Optional height, overriding the intermediate 

340 height (C{meter}). 

341 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

342 

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

344 

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

346 

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

348 

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

350 

351 @example: 

352 

353 >>> p1 = LatLon(52.205, 0.119) 

354 >>> p2 = LatLon(48.857, 2.351) 

355 >>> p = p1.intermediateTo(p2, 0.25) # 51.3721°N, 000.7073°E 

356 ''' 

357 f = Scalar(fraction=fraction) # not high=_1_0 

358 if isnear0(f): # PYCHOK no cover 

359 r = self 

360 elif isnear1(f) and not wrap: # PYCHOK no cover 

361 r = self.others(other) 

362 else: 

363 a1, b1 = self.philam 

364 a2, b2 = self.others(other).philam 

365 

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

367 r = vincentys_(a2, a1, db) 

368 sr = sin(r) 

369 if isnon0(sr): 

370 sa1, ca1, sa2, ca2, \ 

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

372 

373 t = f * r 

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

375 b = sin( t) # / sr superflous 

376 

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

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

379 z = a * sa1 + b * sa2 

380 

381 a = atan2(z, hypot(x, y)) 

382 b = atan2(y, x) 

383 

384 else: # PYCHOK no cover 

385 a = favg(a1, a2, f=f) # coincident 

386 b = favg(b1, b2, f=f) 

387 

388 h = self._havg(other, f=f) if height is None else Height(height) 

389 r = self.classof(degrees90(a), degrees180(b), height=h) 

390 return r 

391 

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

393 '''Compute the intersection point of two paths, each defined 

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

395 

396 @arg end1: End point of this path (L{LatLon}) or the 

397 initial bearing at this point (compass 

398 C{degrees360}). 

399 @arg other: Start point of the other path (L{LatLon}). 

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

401 the initial bearing at the other start point 

402 (compass C{degrees360}). 

403 @kwarg height: Optional height for intersection point, 

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

405 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

406 

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

408 intersection point might be the L{antipode} to 

409 the returned result. 

410 

411 @raise IntersectionError: Ambiguous or infinite intersection 

412 or colinear, parallel or otherwise 

413 non-intersecting paths. 

414 

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

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

417 

418 @raise ValueError: Invalid B{C{height}} or C{null} path. 

419 

420 @example: 

421 

422 >>> p = LatLon(51.8853, 0.2545) 

423 >>> s = LatLon(49.0034, 2.5735) 

424 >>> i = p.intersection(108.547, s, 32.435) # '50.9078°N, 004.5084°E' 

425 ''' 

426 try: 

427 s2 = self.others(other) 

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

429 LatLon=self.classof) 

430 except (TypeError, ValueError) as x: 

431 raise _xError(x, start1=self, end1=end1, other=other, end2=end2) 

432 

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

434 height=None, wrap=True): 

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

436 by a center point and radius. 

437 

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

439 see B{C{radius}}). 

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

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

442 see B{C{radius}}). 

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

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

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

446 B{C{radius}}). 

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

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

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

450 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

451 

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

453 instance. For abutting circles, both intersection 

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

455 

456 @raise IntersectionError: Concentric, antipodal, invalid or 

457 non-intersecting circles. 

458 

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

460 

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

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

463 ''' 

464 try: 

465 c2 = self.others(other) 

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

467 height=height, wrap=wrap, 

468 LatLon=self.classof) 

469 except (TypeError, ValueError) as x: 

470 raise _xError(x, center=self, rad1=rad1, other=other, rad2=rad2) 

471 

472 def isenclosedBy(self, points): 

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

474 

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

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

477 

478 @return: C{True} if this point is inside the polygon or composite, 

479 C{False} otherwise. 

480 

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

482 

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

484 

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

486 

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

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

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

490 ''' 

491 if _MODS.booleans.isBoolean(points): 

492 return points._encloses(self.lat, self.lon) 

493 

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

495 n0 = self._N_vector 

496 

497 v2 = Ps[0]._N_vector 

498 v1 = Ps[1]._N_vector 

499 gc1 = v2.cross(v1) 

500 # check whether this point on same side of all 

501 # polygon edges (to the left or right depending 

502 # on anti-/clockwise polygon direction) 

503 t0 = gc1.angleTo(n0) > PI_2 

504 s0 = None 

505 # get great-circle vector for each edge 

506 for i, p in Ps.enumerate(closed=True): 

507 v2 = p._N_vector 

508 gc = v1.cross(v2) 

509 v1 = v2 

510 

511 t = gc.angleTo(n0) > PI_2 

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

513 return False # outside 

514 

515 # check for convex polygon: angle between 

516 # gc vectors, signed by direction of n0 

517 # (otherwise the test above is not reliable) 

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

519 if s != s0: 

520 if s0 is None: 

521 s0 = s 

522 else: 

523 t = _Fmt.SQUARE(points=i) 

524 raise _ValueError(t, p, txt=_not_(_convex_)) 

525 gc1 = gc 

526 

527 return True # inside 

528 

529 @deprecated_method 

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

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

532 return self.isenclosedBy(points) 

533 

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

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

536 

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

538 @kwarg height: Optional height for midpoint, overriding 

539 the mean height (C{meter}). 

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

541 may be negative or greater than 1.0. 

542 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

543 

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

545 

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

547 

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

549 

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

551 

552 @example: 

553 

554 >>> p1 = LatLon(52.205, 0.119) 

555 >>> p2 = LatLon(48.857, 2.351) 

556 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E' 

557 ''' 

558 if fraction is _0_5: 

559 self.others(other) 

560 

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

562 a1, b1 = self.philam 

563 a2, b2 = other.philam 

564 

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

566 

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

568 

569 x = ca2 * cdb + ca1 

570 y = ca2 * sdb 

571 

572 a = atan2(sa1 + sa2, hypot(x, y)) 

573 b = atan2(y, x) + b1 

574 

575 h = self._havg(other) if height is None else Height(height) 

576 r = self.classof(degrees90(a), degrees180(b), height=h) 

577 else: 

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

579 return r 

580 

581 def nearestOn(self, point1, point2, radius=R_M, **options): 

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

583 

584 Distances are approximated by function L{pygeodesy.equirectangular_}, 

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

586 

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

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

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

590 @kwarg options: Optional keyword arguments for function 

591 L{pygeodesy.equirectangular_}. 

592 

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

594 

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

596 see function L{pygeodesy.equirectangular_}. 

597 

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

599 is not (yet) supported. 

600 

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

602 

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

604 

605 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5} 

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

607 ''' 

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

609 within = _xkwds_pop(options, within=True) 

610 if not within: 

611 notImplemented(self, within=within) 

612 

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

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

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

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

617# return point1 

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

619# if isnear0(d): 

620# return point1 # or point2 

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

622# return point2 

623# f = a / d 

624# if within: 

625# if f > EPS1: 

626# return point2 

627# elif f < EPS: 

628# return point1 

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

630 

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

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

633 **options)[0] 

634 

635 @deprecated_method 

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

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

638 

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

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

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

642 units of B{C{radius}}. 

643 ''' 

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

645 return r.closest, r.distance 

646 

647 def nearestOn3(self, points, closed=False, radius=R_M, **options): 

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

649 

650 Distances are approximated by function L{pygeodesy.equirectangular_}, 

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

652 

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

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

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

656 @kwarg options: Optional keyword arguments for function 

657 L{pygeodesy.equirectangular_}. 

658 

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

660 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_} 

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

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

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

664 function L{pygeodesy.compassAngle}. 

665 

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

667 see function L{pygeodesy.equirectangular_}. 

668 

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

670 

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

672 

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

674 

675 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_} 

676 and L{pygeodesy.nearestOn5}. 

677 ''' 

678 if _LatLon_ in options: 

679 raise _ValueError(**options) 

680 

681 lat, lon, d, c, h = _nearestOn5(self, points, closed=closed, **options) 

682 return NearestOn3Tuple(self.classof(lat, lon, height=h), 

683 degrees2m(d, radius=radius), c) 

684 

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

686 '''Convert this point to C{Karney}-based cartesian (ECEF) 

687 coordinates. 

688 

689 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}} 

690 and other keyword arguments, ignored 

691 if C{B{Cartesian} is None}. Use 

692 C{B{Cartesian}=...} to override 

693 this L{Cartesian} class or specify 

694 C{B{Cartesian}=None}. 

695 

696 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} 

697 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, 

698 C, M, datum)} with C{C} and C{M} if available. 

699 

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

701 ''' 

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

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

704 

705 def _trackDistanceTo3(self, start, end, radius, wrap): 

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

707 ''' 

708 self.others(start=start) 

709 self.others(end=end) 

710 

711 r = Radius_(radius) 

712 r = start.distanceTo(self, r, wrap=wrap) / r 

713 

714 b = radians(start.initialBearingTo(self, wrap=wrap)) 

715 e = radians(start.initialBearingTo(end, wrap=wrap)) 

716 x = asin(sin(r) * sin(b - e)) 

717 return r, x, (e - b) 

718 

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

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

721 

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

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

724 @kwarg radius: Mean earth radius, ellipsoid or datum 

725 (C{meter}, L{Ellipsoid}, L{Ellipsoid2}, 

726 L{Datum} or L{a_f2Tuple}) or C{None}. 

727 @kwarg wrap: Wrap/unroll angular distances (C{bool}). 

728 

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

730 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A, 

731 a, B, b, C, c, D, E)}. 

732 

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

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

735 ''' 

736 self.others(otherB=otherB) 

737 self.others(otherC=otherC) 

738 

739 r = self.philam + otherB.philam + otherC.philam 

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

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

742 

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

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

745 '''Trilaterate three points by area overlap or perimeter intersection 

746 of three corresponding circles. 

747 

748 @arg distance1: Distance to this point (C{meter}, same units 

749 as B{C{radius}}). 

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

751 @arg distance2: Distance to point2 (C{meter}, same units as 

752 B{C{radius}}). 

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

754 @arg distance3: Distance to point3 (C{meter}, same units as 

755 B{C{radius}}). 

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

757 perimeter intersection of the circles (C{bool}). 

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

759 or the I{intersection margin} for C{B{area}=False} 

760 (C{meter}, same units as B{C{radius}}). 

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

762 @kwarg wrap: Wrap/unroll angular distances (C{bool}). 

763 

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

765 with C{min} and C{max} in C{meter}, same units as B{C{eps}}, 

766 the corresponding trilaterated points C{minPoint} and 

767 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number 

768 of trilatered points found for the given B{C{eps}}. 

769 

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

771 max}, C{minPoint I{is} maxPoint} and C{n = 1}. 

772 

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

774 respectively largest I{radial} overlap found. 

775 

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

777 nearest respectively farthest intersection margin. 

778 

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

780 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}} 

781 with the smallest B{C{distance#}} C{min} and C{max} the 

782 largest B{C{distance#}}. 

783 

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

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

786 no intersection or all (near-)concentric for 

787 C{B{area}=False}. 

788 

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

790 

791 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}}, 

792 B{C{distance2}}, B{C{distance3}} or B{C{radius}}. 

793 ''' 

794 return _trilaterate5(self, distance1, 

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

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

797 area=area, radius=radius, eps=eps, wrap=wrap) 

798 

799 

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

801 

802 

803def areaOf(points, radius=R_M, wrap=True): 

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

805 (with the pointsjoined by great circle arcs). 

806 

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

808 or L{BooleanGH}). 

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

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

811 or C{None}. 

812 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

813 

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

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

816 

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

818 

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

820 

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

822 

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

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

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

826 

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

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

829 L{ellipsoidalKarney.areaOf}. 

830 

831 @example: 

832 

833 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1) 

834 >>> areaOf(b) # 8666058750.718977 

835 

836 >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1) 

837 >>> areaOf(c) # 6.18e9 

838 ''' 

839 if _MODS.booleans.isBoolean(points): 

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

841 

842 Ps = _T00.PointsIter(points, loop=1) 

843 p1 = p2 = Ps[0] 

844 a1, b1 = p1.philam 

845 ta1, z1 = tan_2(a1), None 

846 

847 A = Fsum() # mean phi 

848 E = Fsum() # see L{pygeodesy.excessKarney_} 

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

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

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

852 D = Fsum() 

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

854 a2, b2 = p2.philam 

855 db, b2 = unrollPI(b1, b2, wrap=wrap if i else False) 

856 ta2 = tan_2(a2) 

857 A += a2 

858 E += atan2(tan_2(db, points=i) * (ta1 + ta2), 

859 _1_0 + ta1 * ta2) 

860 ta1, b1 = ta2, b2 

861 

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

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

864 if z1 is not None: 

865 D += wrap180(z - z1) # (z - z1 + 540) ... 

866 D += wrap180(z2 - z) # (z2 - z + 540) % 360 - 180 

867 p1, z1 = p2, z2 

868 

869 R = fabs(E * _2_0) 

870 if fabs(D) < _90_0: # ispolar(points) 

871 R = fabs(R - PI2) 

872 if radius: 

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

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

875 return float(R) 

876 

877 

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

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

880 

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

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

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

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

885 

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

887 ''' 

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

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

890 ca *= sr 

891 

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

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

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

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

896 

897 

898def _int3d2(start, end, wrap, _i_, Vector, hs): 

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

900 # and similar logic in .ellipsoidalBaseDI._intersect3 

901 a1, b1 = start.philam 

902 

903 if isscalar(end): # bearing, get pseudo-end point 

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

905 else: # must be a point 

906 start.others(end, name=_end_ + _i_) 

907 hs.append(end.height) 

908 a2, b2 = end.philam 

909 

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

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

912 raise _ValueError(_SPACE_(_path_ + _i_, _null_)) 

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

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

915 -(b1 + b2) * _0_5) 

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

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

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

919 cb12 * cb21 + sb12 * sb21, 

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

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

922 

923 

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

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

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

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

928 

929 

930def intersecant2(center, circle, point, bearing, radius=R_M, exact=False, 

931 height=None, wrap=True): 

932 '''Compute the intersections of a circle and a line. 

933 

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

935 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}}) 

936 or a point on the circle (L{LatLon}). 

937 @arg point: A point in- or outside the circle (L{LatLon}). 

938 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or 

939 a second point on the line (L{LatLon}). 

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

941 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth, 

942 destination and distance, if C{False} use the basic 

943 rhumb methods (C{bool}) or if C{None} use the I{great 

944 circle} methods. 

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

946 conventionally) or C{None}. 

947 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

948 

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

950 each an instance of this class. For a tangent line, each 

951 point C{is} this very instance. 

952 

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

954 

955 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or 

956 B{C{circle}} or B{C{bearing}} invalid. 

957 

958 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}}, 

959 B{C{exact}} or B{C{height}}. 

960 ''' 

961 c = _T00.others(center=center) 

962 p = _T00.others(point=point) 

963 try: 

964 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact, 

965 height=height, wrap=wrap) 

966 except (TypeError, ValueError) as x: 

967 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing, exact=exact) 

968 

969 

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

971 LatLon=None, **LatLon_kwds): 

972 # (INTERNAL) Intersect two (spherical) path, see L{intersection} 

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

974 

975 hs = [start1.height, start2.height] 

976 

977 a1, b1 = start1.philam 

978 a2, b2 = start2.philam 

979 

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

981 r12 = vincentys_(a2, a1, db) 

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

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

984 

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

986 elif isscalar(end1) and isscalar(end2): # both bearings 

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

988 

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

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

991 raise IntersectionError(_parallel_) 

992 # handle domain error for equivalent longitudes, 

993 # see also functions asin_safe and acos_safe at 

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

995 t1, t2 = acos1((sa2 - sa1 * cr12) / x1), \ 

996 acos1((sa1 - sa2 * cr12) / x2) 

997 if sin(db) > 0: 

998 t12, t21 = t1, PI2 - t2 

999 else: 

1000 t12, t21 = PI2 - t1, t2 

1001 t13, t23 = radiansPI2(end1), radiansPI2(end2) 

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

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

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

1005 raise IntersectionError(_infinite_) 

1006 sx3 = sx1 * sx2 

1007# XXX if sx3 < 0: 

1008# XXX raise ValueError(_ambiguous_) 

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

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

1011 

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

1013 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

1017 

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

1019 _N_vector_ = _MODS.nvectorBase._N_vector_ 

1020 

1021 x1, d1 = _int3d2(start1, end1, wrap, _1_, _N_vector_, hs) 

1022 x2, d2 = _int3d2(start2, end2, wrap, _2_, _N_vector_, hs) 

1023 x = x1.cross(x2) 

1024 if x.length < EPS: # [nearly] colinear or parallel paths 

1025 raise IntersectionError(_colinear_) 

1026 a, b = x.philam 

1027 # choose intersection similar to sphericalNvector 

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

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

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

1031 

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

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

1034 intersection, LatLon, LatLon_kwds) 

1035 

1036 

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

1038 LatLon=LatLon, **LatLon_kwds): 

1039 '''Compute the intersection point of two paths, each defined 

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

1041 

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

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

1044 the initial bearing at the first start point 

1045 (compass C{degrees360}). 

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

1047 @arg end2: End point of the second path (L{LatLon}) or 

1048 the initial bearing at the second start point 

1049 (compass C{degrees360}). 

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

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

1052 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

1053 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1057 

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

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

1060 height)}. An alternate intersection point might 

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

1062 

1063 @raise IntersectionError: Ambiguous or infinite intersection 

1064 or colinear, parallel or otherwise 

1065 non-intersecting paths. 

1066 

1067 @raise TypeError: A B{C{start}} or B{C{end}} point not L{LatLon}. 

1068 

1069 @raise ValueError: Invalid B{C{height}} or C{null} path. 

1070 

1071 @example: 

1072 

1073 >>> p = LatLon(51.8853, 0.2545) 

1074 >>> s = LatLon(49.0034, 2.5735) 

1075 >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E' 

1076 ''' 

1077 s1 = _T00.others(start1=start1) 

1078 s2 = _T00.others(start2=start2) 

1079 try: 

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

1081 LatLon=LatLon, **LatLon_kwds) 

1082 except (TypeError, ValueError) as x: 

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

1084 

1085 

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

1087 height=None, wrap=True, 

1088 LatLon=LatLon, **LatLon_kwds): 

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

1090 by a center point and a radius. 

1091 

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

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

1094 see B{C{radius}}). 

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

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

1097 see B{C{radius}}). 

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

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

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

1101 B{C{radius}}). 

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

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

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

1105 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

1106 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1110 

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

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

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

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

1115 

1116 @raise IntersectionError: Concentric, antipodal, invalid or 

1117 non-intersecting circles. 

1118 

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

1120 

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

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

1123 

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

1125 

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

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

1128 ''' 

1129 c1 = _T00.others(center1=center1) 

1130 c2 = _T00.others(center2=center2) 

1131 try: 

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

1133 height=height, wrap=wrap, 

1134 LatLon=LatLon, **LatLon_kwds) 

1135 except (TypeError, ValueError) as x: 

1136 raise _xError(x, center1=center1, rad1=rad1, center2=center2, rad2=rad2) 

1137 

1138 

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

1140 height=None, too_d=None, wrap=True, 

1141 LatLon=LatLon, **LatLon_kwds): 

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

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

1144 

1145 def _dest3(bearing, h): 

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

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

1148 intersections2, LatLon, LatLon_kwds) 

1149 

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

1151 if f: # swapped 

1152 c1, c2 = c2, c1 # PYCHOK swap 

1153 

1154 a1, b1 = c1.philam 

1155 a2, b2 = c2.philam 

1156 

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

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

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

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

1161 

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

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

1164 if r < _0_0: 

1165 raise _ValueError(eps=r) 

1166 

1167 x = fsum_(r1, r2, -d) # overlap 

1168 if x > max(r, EPS): 

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

1170 x = sd * sr1 

1171 if isnear0(x): 

1172 raise _ValueError(_invalid_) 

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

1174 

1175 elif x < r: # PYCHOK no cover 

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

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

1178 

1179 if height is None: # "radical height" 

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

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

1182 else: 

1183 h = Height(height) 

1184 

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

1186 if x < EPS4: # externally ... 

1187 r = _dest3(b, h) 

1188 elif x > _PI_EPS4: # internally ... 

1189 r = _dest3(b + PI, h) 

1190 else: 

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

1192 return r, r # ... abutting circles 

1193 

1194 

1195@deprecated_function 

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

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

1198 ''' 

1199 return ispolar(points, wrap=wrap) 

1200 

1201 

1202def _LL3Tuple(lat, lon, height, func, LatLon, LatLon_kwds): 

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

1204 ''' 

1205 n = func.__name__ 

1206 if LatLon is None: 

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

1208 else: 

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

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

1211 return r 

1212 

1213 

1214def meanOf(points, height=None, LatLon=LatLon, **LatLon_kwds): 

1215 '''Compute the geographic mean of several points. 

1216 

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

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

1219 height (C{meter}). 

1220 @kwarg LatLon: Optional class to return the mean point (L{LatLon}) 

1221 or C{None}. 

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

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

1224 

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

1226 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}} 

1227 is C{None}. 

1228 

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

1230 

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

1232 ''' 

1233 # geographic mean 

1234 m = sumOf(p._N_vector for p in _T00.PointsIter(points).iterate(closed=False)) 

1235 lat, lon = m._N_vector.latlon 

1236 

1237 if height is None: 

1238 h = fmean(p.height for p in _T00.PointsIter(points).iterate(closed=False)) 

1239 else: # PYCHOK no cover 

1240 h = Height(height) 

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

1242 

1243 

1244@deprecated_function 

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

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

1247 

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

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

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

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

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

1253 ''' 

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

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

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

1257 return ll, d 

1258 

1259 

1260def nearestOn3(point, points, closed=False, radius=R_M, 

1261 LatLon=LatLon, **options): 

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

1263 

1264 Distances are I{approximated} using function L{pygeodesy.equirectangular_}, 

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

1266 

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

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

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

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

1271 @kwarg LatLon: Optional class to return the closest point (L{LatLon}) 

1272 or C{None}. 

1273 @kwarg options: Optional keyword arguments for function 

1274 L{pygeodesy.equirectangular_}. 

1275 

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

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

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

1279 is the L{pygeodesy.equirectangular_} distance between the 

1280 C{closest} and the given B{C{point}} converted to C{meter}, 

1281 same units as B{C{radius}}. The C{angle} from the given 

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

1283 like function L{pygeodesy.compassAngle}. The C{height} is 

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

1285 

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

1287 see function L{pygeodesy.equirectangular_}. 

1288 

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

1290 

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

1292 

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

1294 

1295 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}. 

1296 ''' 

1297 lat, lon, d, c, h = _nearestOn5(point, points, closed=closed, 

1298 LatLon=None, **options) 

1299 d = degrees2m(d, radius=radius) 

1300 n = nearestOn3.__name__ 

1301 r = LatLon3Tuple(lat, lon, h, name=n) if LatLon is None else \ 

1302 LatLon(lat, lon, height=h, name=n) 

1303 return NearestOn3Tuple(r, d, c, name=n) 

1304 

1305 

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

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

1308 (with great circle arcs joining the points). 

1309 

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

1311 or L{BooleanGH}). 

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

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

1314 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

1315 

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

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

1318 

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

1320 

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

1322 

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

1324 C{B{points}} a composite. 

1325 

1326 @note: Distances are based on function L{pygeodesy.vincentys_}. 

1327 

1328 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf} 

1329 and L{ellipsoidalKarney.perimeterOf}. 

1330 ''' 

1331 def _rads(points, closed, wrap): # angular edge lengths in radians 

1332 Ps = _T00.PointsIter(points, loop=1) 

1333 a1, b1 = Ps[0].philam 

1334 for i, p in Ps.enumerate(closed=closed): 

1335 a2, b2 = p.philam 

1336 db, b2 = unrollPI(b1, b2, wrap=wrap if i else False) 

1337 yield vincentys_(a2, a1, db) 

1338 a1, b1 = a2, b2 

1339 

1340 if _MODS.booleans.isBoolean(points): 

1341 if not closed: 

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

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

1344 else: 

1345 r = fsum(_rads(points, closed, wrap), floats=True) 

1346 return _r2m(r, radius) 

1347 

1348 

1349def _r2m(r, radius): 

1350 '''(INTERNAL) Angular distance in C{radians} to C{meter}. 

1351 ''' 

1352 if radius is not None: # not in (None, _0_0) 

1353 r *= R_M if radius is R_M else Radius(radius) 

1354 return r 

1355 

1356 

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

1358 excess=excessAbc_, 

1359 wrap=False): 

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

1361 

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

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

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

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

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

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

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

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

1370 or C{None}. 

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

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

1373 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}). 

1374 

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

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

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

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

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

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

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

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

1383 ''' 

1384 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA), 

1385 Phi_(latB=latB), Lam_(lonB=lonB), 

1386 Phi_(latC=latC), Lam_(lonC=lonC), 

1387 excess=excess, wrap=wrap) 

1388 return _t7Tuple(t, radius) 

1389 

1390 

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

1392 wrap=False): 

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

1394 excess} of a (spherical) triangle. 

1395 

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

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

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

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

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

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

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

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

1404 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}). 

1405 

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

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

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

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

1410 ''' 

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

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

1413 a = vincentys_(phiC, phiB, d) 

1414 return a, (phiB, lamB, phiC, lamC, phiA, lamA) 

1415 

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

1417 s = sb * sc 

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

1419 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2's 

1420 

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

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

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

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

1425 

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

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

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

1429 

1430 D = fsum_(PI2, -a, -b, -c, floats=True) # deficit aka defect 

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

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

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

1434 

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

1436 

1437 

1438def _t7Tuple(t, radius): 

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

1440 ''' 

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

1442 r = radius if isscalar(radius) else \ 

1443 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1446 B, (r * t.b), 

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

1448 return t 

1449 

1450 

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

1452 areaOf, # functions 

1453 intersecant2, intersection, intersections2, ispolar, 

1454 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1455 meanOf, 

1456 nearestOn2, nearestOn3, 

1457 perimeterOf, 

1458 sumOf, # == vector3d.sumOf 

1459 triangle7, triangle8_) 

1460 

1461# **) MIT License 

1462# 

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

1464# 

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

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

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

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

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

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

1471# 

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

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

1474# 

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

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

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

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

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

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

1481# OTHER DEALINGS IN THE SOFTWARE.