Coverage for pygeodesy/sphericalTrigonometry.py: 96%

374 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-03-31 10:52 -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.03.30' 

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 encloses this point. 

474 

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

476 

477 @return: C{True} if the polygon encloses this point, 

478 C{False} otherwise. 

479 

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

481 

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

483 

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

485 

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

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

488 enclose a pole or wrap around the earth longitudinally. 

489 

490 @example: 

491 

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

493 >>> p = LatLon(45,1, 1.1) 

494 >>> inside = p.isEnclosedBy(b) # True 

495 ''' 

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

497 n0 = self._N_vector 

498 

499 v2 = Ps[0]._N_vector 

500 v1 = Ps[1]._N_vector 

501 gc1 = v2.cross(v1) 

502 # check whether this point on same side of all 

503 # polygon edges (to the left or right depending 

504 # on anti-/clockwise polygon direction) 

505 t0 = gc1.angleTo(n0) > PI_2 

506 s0 = None 

507 # get great-circle vector for each edge 

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

509 v2 = p._N_vector 

510 gc = v1.cross(v2) 

511 v1 = v2 

512 

513 t = gc.angleTo(n0) > PI_2 

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

515 return False # outside 

516 

517 # check for convex polygon: angle between 

518 # gc vectors, signed by direction of n0 

519 # (otherwise the test above is not reliable) 

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

521 if s != s0: 

522 if s0 is None: 

523 s0 = s 

524 else: 

525 t = _Fmt.SQUARE(points=i) 

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

527 gc1 = gc 

528 

529 return True # inside 

530 

531 @deprecated_method 

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

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

534 return self.isenclosedBy(points) 

535 

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

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

538 

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

540 @kwarg height: Optional height for midpoint, overriding 

541 the mean height (C{meter}). 

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

543 may be negative or greater than 1.0. 

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

545 

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

547 

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

549 

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

551 

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

553 

554 @example: 

555 

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

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

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

559 ''' 

560 if fraction is _0_5: 

561 self.others(other) 

562 

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

564 a1, b1 = self.philam 

565 a2, b2 = other.philam 

566 

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

568 

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

570 

571 x = ca2 * cdb + ca1 

572 y = ca2 * sdb 

573 

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

575 b = atan2(y, x) + b1 

576 

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

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

579 else: 

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

581 return r 

582 

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

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

585 

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

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

588 

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

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

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

592 @kwarg options: Optional keyword arguments for function 

593 L{pygeodesy.equirectangular_}. 

594 

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

596 

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

598 see function L{pygeodesy.equirectangular_}. 

599 

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

601 is not (yet) supported. 

602 

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

604 

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

606 

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

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

609 ''' 

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

611 within = _xkwds_pop(options, within=True) 

612 if not within: 

613 notImplemented(self, within=within) 

614 

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

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

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

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

619# return point1 

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

621# if isnear0(d): 

622# return point1 # or point2 

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

624# return point2 

625# f = a / d 

626# if within: 

627# if f > EPS1: 

628# return point2 

629# elif f < EPS: 

630# return point1 

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

632 

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

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

635 **options)[0] 

636 

637 @deprecated_method 

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

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

640 

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

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

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

644 units of B{C{radius}}. 

645 ''' 

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

647 return r.closest, r.distance 

648 

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

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

651 

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

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

654 

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

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

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

658 @kwarg options: Optional keyword arguments for function 

659 L{pygeodesy.equirectangular_}. 

660 

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

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

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

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

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

666 function L{pygeodesy.compassAngle}. 

667 

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

669 see function L{pygeodesy.equirectangular_}. 

670 

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

672 

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

674 

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

676 

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

678 and L{pygeodesy.nearestOn5}. 

679 ''' 

680 if _LatLon_ in options: 

681 raise _ValueError(**options) 

682 

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

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

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

686 

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

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

689 coordinates. 

690 

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

692 and other keyword arguments, ignored 

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

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

695 this L{Cartesian} class or specify 

696 C{B{Cartesian}=None}. 

697 

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

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

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

701 

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

703 ''' 

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

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

706 

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

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

709 ''' 

710 self.others(start=start) 

711 self.others(end=end) 

712 

713 r = Radius_(radius) 

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

715 

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

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

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

719 return r, x, (e - b) 

720 

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

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

723 

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

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

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

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

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

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

730 

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

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

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

734 

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

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

737 ''' 

738 self.others(otherB=otherB) 

739 self.others(otherC=otherC) 

740 

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

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

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

744 

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

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

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

748 of three corresponding circles. 

749 

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

751 as B{C{radius}}). 

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

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

754 B{C{radius}}). 

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

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

757 B{C{radius}}). 

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

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

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

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

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

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

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

765 

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

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

768 the corresponding trilaterated points C{minPoint} and 

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

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

771 

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

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

774 

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

776 respectively largest I{radial} overlap found. 

777 

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

779 nearest respectively farthest intersection margin. 

780 

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

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

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

784 largest B{C{distance#}}. 

785 

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

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

788 no intersection or all (near-)concentric for 

789 C{B{area}=False}. 

790 

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

792 

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

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

795 ''' 

796 return _trilaterate5(self, distance1, 

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

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

799 area=area, radius=radius, eps=eps, wrap=wrap) 

800 

801 

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

803 

804 

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

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

807 (with the pointsjoined by great circle arcs). 

808 

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

810 or L{BooleanGH}). 

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

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

813 or C{None}. 

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

815 

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

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

818 

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

820 

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

822 

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

824 

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

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

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

828 

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

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

831 L{ellipsoidalKarney.areaOf}. 

832 

833 @example: 

834 

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

836 >>> areaOf(b) # 8666058750.718977 

837 

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

839 >>> areaOf(c) # 6.18e9 

840 ''' 

841 if _MODS.booleans.isBoolean(points): 

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

843 

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

845 p1 = p2 = Ps[0] 

846 a1, b1 = p1.philam 

847 ta1, z1 = tan_2(a1), None 

848 

849 A = Fsum() # mean phi 

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

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

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

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

854 D = Fsum() 

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

856 a2, b2 = p2.philam 

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

858 ta2 = tan_2(a2) 

859 A += a2 

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

861 _1_0 + ta1 * ta2) 

862 ta1, b1 = ta2, b2 

863 

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

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

866 if z1 is not None: 

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

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

869 p1, z1 = p2, z2 

870 

871 R = fabs(E * _2_0) 

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

873 R = fabs(R - PI2) 

874 if radius: 

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

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

877 return float(R) 

878 

879 

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

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

882 

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

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

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

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

887 

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

889 ''' 

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

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

892 ca *= sr 

893 

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

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

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

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

898 

899 

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

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

902 # and similar logic in .ellipsoidalBaseDI._intersect3 

903 a1, b1 = start.philam 

904 

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

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

907 else: # must be a point 

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

909 hs.append(end.height) 

910 a2, b2 = end.philam 

911 

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

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

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

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

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

917 -(b1 + b2) * _0_5) 

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

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

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

921 cb12 * cb21 + sb12 * sb21, 

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

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

924 

925 

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

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

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

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

930 

931 

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

933 height=None, wrap=True): 

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

935 

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

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

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

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

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

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

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

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

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

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

946 circle} methods. 

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

948 conventionally) or C{None}. 

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

950 

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

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

953 point C{is} this very instance. 

954 

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

956 

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

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

959 

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

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

962 ''' 

963 c = _T00.others(center=center) 

964 p = _T00.others(point=point) 

965 try: 

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

967 height=height, wrap=wrap) 

968 except (TypeError, ValueError) as x: 

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

970 

971 

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

973 LatLon=None, **LatLon_kwds): 

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

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

976 

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

978 

979 a1, b1 = start1.philam 

980 a2, b2 = start2.philam 

981 

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

983 r12 = vincentys_(a2, a1, db) 

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

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

986 

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

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

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

990 

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

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

993 raise IntersectionError(_parallel_) 

994 # handle domain error for equivalent longitudes, 

995 # see also functions asin_safe and acos_safe at 

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

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

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

999 if sin(db) > 0: 

1000 t12, t21 = t1, PI2 - t2 

1001 else: 

1002 t12, t21 = PI2 - t1, t2 

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

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

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

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

1007 raise IntersectionError(_infinite_) 

1008 sx3 = sx1 * sx2 

1009# XXX if sx3 < 0: 

1010# XXX raise ValueError(_ambiguous_) 

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

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

1013 

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

1015 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

1019 

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

1021 _N_vector_ = _MODS.nvectorBase._N_vector_ 

1022 

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

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

1025 x = x1.cross(x2) 

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

1027 raise IntersectionError(_colinear_) 

1028 a, b = x.philam 

1029 # choose intersection similar to sphericalNvector 

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

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

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

1033 

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

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

1036 intersection, LatLon, LatLon_kwds) 

1037 

1038 

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

1040 LatLon=LatLon, **LatLon_kwds): 

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

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

1043 

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

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

1046 the initial bearing at the first start point 

1047 (compass C{degrees360}). 

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

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

1050 the initial bearing at the second start point 

1051 (compass C{degrees360}). 

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

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

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

1055 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1059 

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

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

1062 height)}. An alternate intersection point might 

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

1064 

1065 @raise IntersectionError: Ambiguous or infinite intersection 

1066 or colinear, parallel or otherwise 

1067 non-intersecting paths. 

1068 

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

1070 

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

1072 

1073 @example: 

1074 

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

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

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

1078 ''' 

1079 s1 = _T00.others(start1=start1) 

1080 s2 = _T00.others(start2=start2) 

1081 try: 

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

1083 LatLon=LatLon, **LatLon_kwds) 

1084 except (TypeError, ValueError) as x: 

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

1086 

1087 

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

1089 height=None, wrap=True, 

1090 LatLon=LatLon, **LatLon_kwds): 

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

1092 by a center point and a radius. 

1093 

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

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

1096 see B{C{radius}}). 

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

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

1099 see B{C{radius}}). 

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

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

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

1103 B{C{radius}}). 

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

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

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

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

1108 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1112 

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

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

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

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

1117 

1118 @raise IntersectionError: Concentric, antipodal, invalid or 

1119 non-intersecting circles. 

1120 

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

1122 

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

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

1125 

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

1127 

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

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

1130 ''' 

1131 c1 = _T00.others(center1=center1) 

1132 c2 = _T00.others(center2=center2) 

1133 try: 

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

1135 height=height, wrap=wrap, 

1136 LatLon=LatLon, **LatLon_kwds) 

1137 except (TypeError, ValueError) as x: 

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

1139 

1140 

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

1142 height=None, too_d=None, wrap=True, 

1143 LatLon=LatLon, **LatLon_kwds): 

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

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

1146 

1147 def _dest3(bearing, h): 

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

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

1150 intersections2, LatLon, LatLon_kwds) 

1151 

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

1153 if f: # swapped 

1154 c1, c2 = c2, c1 # PYCHOK swap 

1155 

1156 a1, b1 = c1.philam 

1157 a2, b2 = c2.philam 

1158 

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

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

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

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

1163 

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

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

1166 if r < _0_0: 

1167 raise _ValueError(eps=r) 

1168 

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

1170 if x > max(r, EPS): 

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

1172 x = sd * sr1 

1173 if isnear0(x): 

1174 raise _ValueError(_invalid_) 

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

1176 

1177 elif x < r: # PYCHOK no cover 

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

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

1180 

1181 if height is None: # "radical height" 

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

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

1184 else: 

1185 h = Height(height) 

1186 

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

1188 if x < EPS4: # externally ... 

1189 r = _dest3(b, h) 

1190 elif x > _PI_EPS4: # internally ... 

1191 r = _dest3(b + PI, h) 

1192 else: 

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

1194 return r, r # ... abutting circles 

1195 

1196 

1197@deprecated_function 

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

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

1200 ''' 

1201 return ispolar(points, wrap=wrap) 

1202 

1203 

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

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

1206 ''' 

1207 n = func.__name__ 

1208 if LatLon is None: 

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

1210 else: 

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

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

1213 return r 

1214 

1215 

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

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

1218 

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

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

1221 height (C{meter}). 

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

1223 or C{None}. 

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

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

1226 

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

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

1229 is C{None}. 

1230 

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

1232 

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

1234 ''' 

1235 # geographic mean 

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

1237 lat, lon = m._N_vector.latlon 

1238 

1239 if height is None: 

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

1241 else: # PYCHOK no cover 

1242 h = Height(height) 

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

1244 

1245 

1246@deprecated_function 

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

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

1249 

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

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

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

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

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

1255 ''' 

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

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

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

1259 return ll, d 

1260 

1261 

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

1263 LatLon=LatLon, **options): 

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

1265 

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

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

1268 

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

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

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

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

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

1274 or C{None}. 

1275 @kwarg options: Optional keyword arguments for function 

1276 L{pygeodesy.equirectangular_}. 

1277 

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

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

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

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

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

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

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

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

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

1287 

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

1289 see function L{pygeodesy.equirectangular_}. 

1290 

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

1292 

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

1294 

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

1296 

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

1298 ''' 

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

1300 LatLon=None, **options) 

1301 d = degrees2m(d, radius=radius) 

1302 n = nearestOn3.__name__ 

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

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

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

1306 

1307 

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

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

1310 (with great circle arcs joining the points). 

1311 

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

1313 or L{BooleanGH}). 

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

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

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

1317 

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

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

1320 

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

1322 

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

1324 

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

1326 C{B{points}} a composite. 

1327 

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

1329 

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

1331 and L{ellipsoidalKarney.perimeterOf}. 

1332 ''' 

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

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

1335 a1, b1 = Ps[0].philam 

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

1337 a2, b2 = p.philam 

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

1339 yield vincentys_(a2, a1, db) 

1340 a1, b1 = a2, b2 

1341 

1342 if _MODS.booleans.isBoolean(points): 

1343 if not closed: 

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

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

1346 else: 

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

1348 return _r2m(r, radius) 

1349 

1350 

1351def _r2m(r, radius): 

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

1353 ''' 

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

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

1356 return r 

1357 

1358 

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

1360 excess=excessAbc, 

1361 wrap=False): 

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

1363 

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

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

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

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

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

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

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

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

1372 or C{None}. 

1373 @kwarg excess: I{Spherical excess} callable (L{excessAbc}, 

1374 L{excessGirard} or L{excessLHuilier}). 

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

1376 

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

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

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

1380 in I{square} C{meter} or units of B{C{radius}} 

1381 I{squared} or if B{C{radius}} is C{None}, a 

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

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

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

1385 ''' 

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

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

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

1389 excess=excess, wrap=wrap) 

1390 return _t7Tuple(t, radius) 

1391 

1392 

1393def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc, 

1394 wrap=False): 

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

1396 excess} of a (spherical) triangle. 

1397 

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

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

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

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

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

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

1404 @kwarg excess: I{Spherical excess} callable (L{excessAbc}, 

1405 L{excessGirard} or L{excessLHuilier}). 

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

1407 

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

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

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

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

1412 ''' 

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

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

1415 a = vincentys_(phiC, phiB, d) 

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

1417 

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

1419 s = sb * sc 

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

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

1422 

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

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

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

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

1427 

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

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

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

1431 

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

1433 E = excessGirard(A, B, C) if excess in (excessGirard, True) else ( 

1434 excessLHuilier(a, b, c) if excess in (excessLHuilier, False) else 

1435 excessAbc(*max((A, b, c), (B, c, a), (C, a, b)))) 

1436 

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

1438 

1439 

1440def _t7Tuple(t, radius): 

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

1442 ''' 

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

1444 r = radius if isscalar(radius) else \ 

1445 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1448 B, (r * t.b), 

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

1450 return t 

1451 

1452 

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

1454 areaOf, # functions 

1455 intersecant2, intersection, intersections2, ispolar, 

1456 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1457 meanOf, 

1458 nearestOn2, nearestOn3, 

1459 perimeterOf, 

1460 sumOf, # == vector3d.sumOf 

1461 triangle7, triangle8_) 

1462 

1463# **) MIT License 

1464# 

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

1466# 

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

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

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

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

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

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

1473# 

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

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

1476# 

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

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

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

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

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

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

1483# OTHER DEALINGS IN THE SOFTWARE.