Coverage for pygeodesy/sphericalTrigonometry.py: 95%

375 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-12 11:45 -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_, _invalid_, \ 

34 _LatLon_, _line_, _near_, _not_, _null_, _points_, \ 

35 _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.11' 

58 

59_parallel_ = 'parallel' 

60 

61_PI_EPS4 = PI - EPS4 

62if _PI_EPS4 >= PI: 

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

64 

65 

66class Cartesian(CartesianSphericalBase): 

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

68 spherical, geodetic L{LatLon}. 

69 ''' 

70 

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

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

73 

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

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

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

77 

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

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

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

81 

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

83 ''' 

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

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

86 

87 

88class LatLon(LatLonSphericalBase): 

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

90 

91 @example: 

92 

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

94 ''' 

95 

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

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

98 the closest point on the great circle line defined by a 

99 start and an end point. 

100 

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

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

103 from the start point to the point where the perpendicular 

104 crosses the line. 

105 

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

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

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

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

110 

111 @return: Distance along the great circle line (C{meter}, 

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

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

114 toward the B{C{end}} point of the line or negative 

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

116 

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

118 

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

120 

121 @example: 

122 

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

124 

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

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

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

128 ''' 

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

130 cx = cos(x) 

131 return _0_0 if isnear0(cx) else \ 

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

133 

134 @deprecated_method 

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

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

137 ''' 

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

139 

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

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

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

143 

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

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

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

147 

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

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

150 ''' 

151 self.others(other) 

152 

153 a1, b1 = self.philam 

154 a2, b2 = other.philam 

155 

156 a = radians(lat) 

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

158 

159 sa, ca, sa1, ca1, \ 

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

161 sa1 *= ca2 * ca 

162 

163 x = sa1 * sdb 

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

165 z = ca1 * sdb * ca2 * sa 

166 

167 h = hypot(x, y) 

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

169 return None # great circle doesn't reach latitude 

170 

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

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

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

174 

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

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

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

178 

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

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

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

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

183 

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

185 or positive if to the right of the line) (C{meter}, 

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

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

188 

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

190 

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

192 

193 @example: 

194 

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

196 

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

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

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

200 ''' 

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

202 return _r2m(x, radius) 

203 

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

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

206 travelled the given distance on the given initial bearing. 

207 

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

209 B{C{radius}}). 

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

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

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

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

214 

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

216 

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

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

219 

220 @example: 

221 

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

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

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

225 ''' 

226 a, b = self.philam 

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

228 

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

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

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

232 

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

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

235 

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

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

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

239 

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

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

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

243 

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

245 

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

247 

248 @example: 

249 

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

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

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

253 ''' 

254 self.others(other) 

255 

256 a1, b1 = self.philam 

257 a2, b2 = other.philam 

258 

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

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

261 

262# @Property_RO 

263# def Ecef(self): 

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

265# ''' 

266# return _MODS.ecef.EcefKarney 

267 

268 def greatCircle(self, bearing): 

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

270 on the given initial bearing from this point. 

271 

272 Direction of vector is such that initial bearing vector 

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

274 

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

276 

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

278 

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

280 

281 @example: 

282 

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

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

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

286 ''' 

287 a, b = self.philam 

288 

289 t = Bearing_(bearing) 

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

291 

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

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

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

295 

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

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

298 to an other point. 

299 

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

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

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

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

304 C{sphericalNvector.LatLon.initialBearingTo}. 

305 

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

307 

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

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

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

311 

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

313 

314 @example: 

315 

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

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

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

319 ''' 

320 self.others(other) 

321 

322 a1, b1 = self.philam 

323 a2, b2 = other.philam 

324 

325 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

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

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

328 

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

330 

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

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

333 and an other point. 

334 

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

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

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

338 @kwarg height: Optional height, overriding the intermediate 

339 height (C{meter}). 

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

341 

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

343 

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

345 

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

347 

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

349 

350 @example: 

351 

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

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

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

355 ''' 

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

357 if isnear0(f): # PYCHOK no cover 

358 r = self 

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

360 r = self.others(other) 

361 else: 

362 a1, b1 = self.philam 

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

364 

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

366 r = vincentys_(a2, a1, db) 

367 sr = sin(r) 

368 if isnon0(sr): 

369 sa1, ca1, sa2, ca2, \ 

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

371 

372 t = f * r 

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

374 b = sin( t) # / sr superflous 

375 

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

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

378 z = a * sa1 + b * sa2 

379 

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

381 b = atan2(y, x) 

382 

383 else: # PYCHOK no cover 

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

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

386 

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

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

389 return r 

390 

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

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

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

394 

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

396 initial bearing at this point (compass 

397 C{degrees360}). 

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

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

400 the initial bearing at the other start point 

401 (compass C{degrees360}). 

402 @kwarg height: Optional height for intersection point, 

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

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

405 

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

407 intersection point might be the L{antipode} to 

408 the returned result. 

409 

410 @raise IntersectionError: Ambiguous or infinite intersection 

411 or colinear, parallel or otherwise 

412 non-intersecting lines. 

413 

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

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

416 

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

418 

419 @example: 

420 

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

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

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

424 ''' 

425 try: 

426 s2 = self.others(other) 

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

428 LatLon=self.classof) 

429 except (TypeError, ValueError) as x: 

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

431 

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

433 height=None, wrap=True): 

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

435 by a center point and radius. 

436 

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

438 see B{C{radius}}). 

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

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

441 see B{C{radius}}). 

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

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

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

445 B{C{radius}}). 

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

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

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

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

450 

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

452 instance. For abutting circles, both intersection 

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

454 

455 @raise IntersectionError: Concentric, antipodal, invalid or 

456 non-intersecting circles. 

457 

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

459 

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

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

462 ''' 

463 try: 

464 c2 = self.others(other) 

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

466 height=height, wrap=wrap, 

467 LatLon=self.classof) 

468 except (TypeError, ValueError) as x: 

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

470 

471 def isenclosedBy(self, points): 

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

473 

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

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

476 

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

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 I{longitudinally}. 

489 ''' 

490 if _MODS.booleans.isBoolean(points): 

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

492 

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

494 n0 = self._N_vector 

495 

496 v2 = Ps[0]._N_vector 

497 v1 = Ps[1]._N_vector 

498 gc1 = v2.cross(v1) 

499 # check whether this point on same side of all 

500 # polygon edges (to the left or right depending 

501 # on anti-/clockwise polygon direction) 

502 t0 = gc1.angleTo(n0) > PI_2 

503 s0 = None 

504 # get great-circle vector for each edge 

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

506 v2 = p._N_vector 

507 gc = v1.cross(v2) 

508 v1 = v2 

509 

510 t = gc.angleTo(n0) > PI_2 

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

512 return False # outside 

513 

514 # check for convex polygon: angle between 

515 # gc vectors, signed by direction of n0 

516 # (otherwise the test above is not reliable) 

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

518 if s != s0: 

519 if s0 is None: 

520 s0 = s 

521 else: 

522 t = _Fmt.SQUARE(points=i) 

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

524 gc1 = gc 

525 

526 return True # inside 

527 

528 @deprecated_method 

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

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

531 return self.isenclosedBy(points) 

532 

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

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

535 

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

537 @kwarg height: Optional height for midpoint, overriding 

538 the mean height (C{meter}). 

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

540 may be negative or greater than 1.0. 

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

542 

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

544 

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

546 

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

548 

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

550 

551 @example: 

552 

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

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

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

556 ''' 

557 if fraction is _0_5: 

558 self.others(other) 

559 

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

561 a1, b1 = self.philam 

562 a2, b2 = other.philam 

563 

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

565 

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

567 

568 x = ca2 * cdb + ca1 

569 y = ca2 * sdb 

570 

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

572 b = atan2(y, x) + b1 

573 

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

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

576 else: 

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

578 return r 

579 

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

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

582 

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

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

585 

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

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

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

589 @kwarg options: Optional keyword arguments for function 

590 L{pygeodesy.equirectangular_}. 

591 

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

593 

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

595 see function L{pygeodesy.equirectangular_}. 

596 

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

598 is not (yet) supported. 

599 

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

601 

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

603 

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

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

606 ''' 

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

608 within = _xkwds_pop(options, within=True) 

609 if not within: 

610 notImplemented(self, within=within) 

611 

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

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

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

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

616# return point1 

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

618# if isnear0(d): 

619# return point1 # or point2 

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

621# return point2 

622# f = a / d 

623# if within: 

624# if f > EPS1: 

625# return point2 

626# elif f < EPS: 

627# return point1 

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

629 

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

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

632 **options)[0] 

633 

634 @deprecated_method 

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

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

637 

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

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

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

641 units of B{C{radius}}. 

642 ''' 

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

644 return r.closest, r.distance 

645 

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

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

648 

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

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

651 

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

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

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

655 @kwarg options: Optional keyword arguments for function 

656 L{pygeodesy.equirectangular_}. 

657 

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

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

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

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

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

663 function L{pygeodesy.compassAngle}. 

664 

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

666 see function L{pygeodesy.equirectangular_}. 

667 

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

669 

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

671 

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

673 

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

675 and L{pygeodesy.nearestOn5}. 

676 ''' 

677 if _LatLon_ in options: 

678 raise _ValueError(**options) 

679 

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

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

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

683 

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

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

686 coordinates. 

687 

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

689 and other keyword arguments, ignored 

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

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

692 this L{Cartesian} class or specify 

693 C{B{Cartesian}=None}. 

694 

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

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

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

698 

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

700 ''' 

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

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

703 

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

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

706 ''' 

707 self.others(start=start) 

708 self.others(end=end) 

709 

710 r = Radius_(radius) 

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

712 

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

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

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

716 return r, x, (e - b) 

717 

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

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

720 

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

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

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

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

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

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

727 

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

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

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

731 

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

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

734 ''' 

735 self.others(otherB=otherB) 

736 self.others(otherC=otherC) 

737 

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

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

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

741 

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

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

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

745 of three corresponding circles. 

746 

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

748 as B{C{radius}}). 

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

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

751 B{C{radius}}). 

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

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

754 B{C{radius}}). 

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

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

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

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

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

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

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

762 

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

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

765 the corresponding trilaterated points C{minPoint} and 

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

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

768 

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

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

771 

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

773 respectively largest I{radial} overlap found. 

774 

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

776 nearest respectively farthest intersection margin. 

777 

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

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

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

781 largest B{C{distance#}}. 

782 

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

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

785 no intersection or all (near-)concentric for 

786 C{B{area}=False}. 

787 

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

789 

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

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

792 ''' 

793 return _trilaterate5(self, distance1, 

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

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

796 area=area, radius=radius, eps=eps, wrap=wrap) 

797 

798 

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

800 

801 

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

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

804 (with the pointsjoined by great circle arcs). 

805 

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

807 or L{BooleanGH}). 

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

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

810 or C{None}. 

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

812 

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

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

815 

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

817 

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

819 

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

821 

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

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

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

825 

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

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

828 L{ellipsoidalKarney.areaOf}. 

829 

830 @example: 

831 

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

833 >>> areaOf(b) # 8666058750.718977 

834 

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

836 >>> areaOf(c) # 6.18e9 

837 ''' 

838 if _MODS.booleans.isBoolean(points): 

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

840 

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

842 p1 = p2 = Ps[0] 

843 a1, b1 = p1.philam 

844 ta1, z1 = tan_2(a1), None 

845 

846 A = Fsum() # mean phi 

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

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

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

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

851 D = Fsum() 

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

853 a2, b2 = p2.philam 

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

855 ta2 = tan_2(a2) 

856 A += a2 

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

858 _1_0 + ta1 * ta2) 

859 ta1, b1 = ta2, b2 

860 

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

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

863 if z1 is not None: 

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

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

866 p1, z1 = p2, z2 

867 

868 R = fabs(E * _2_0) 

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

870 R = fabs(R - PI2) 

871 if radius: 

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

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

874 return float(R) 

875 

876 

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

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

879 

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

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

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

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

884 

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

886 ''' 

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

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

889 ca *= sr 

890 

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

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

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

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

895 

896 

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

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

899 # and similar logic in .ellipsoidalBaseDI._intersect3 

900 a1, b1 = start.philam 

901 

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

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

904 else: # must be a point 

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

906 hs.append(end.height) 

907 a2, b2 = end.philam 

908 

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

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

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

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

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

914 -(b1 + b2) * _0_5) 

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

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

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

918 cb12 * cb21 + sb12 * sb21, 

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

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

921 

922 

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

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

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

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

927 

928 

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

930 height=None, wrap=True): 

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

932 

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

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

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

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

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

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

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

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

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

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

943 circle} methods. 

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

945 conventionally) or C{None}. 

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

947 

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

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

950 point C{is} this very instance. 

951 

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

953 

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

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

956 

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

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

959 ''' 

960 c = _T00.others(center=center) 

961 p = _T00.others(point=point) 

962 try: 

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

964 height=height, wrap=wrap) 

965 except (TypeError, ValueError) as x: 

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

967 

968 

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

970 LatLon=None, **LatLon_kwds): 

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

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

973 

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

975 

976 a1, b1 = start1.philam 

977 a2, b2 = start2.philam 

978 

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

980 r12 = vincentys_(a2, a1, db) 

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

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

983 

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

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

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

987 

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

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

990 raise IntersectionError(_parallel_) 

991 # handle domain error for equivalent longitudes, 

992 # see also functions asin_safe and acos_safe at 

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

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

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

996 if sin(db) > 0: 

997 t12, t21 = t1, PI2 - t2 

998 else: 

999 t12, t21 = PI2 - t1, t2 

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

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

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

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

1004 raise IntersectionError(_infinite_) 

1005 sx3 = sx1 * sx2 

1006# XXX if sx3 < 0: 

1007# XXX raise ValueError(_ambiguous_) 

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

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

1010 

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

1012 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

1016 

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

1018 _N_vector_ = _MODS.nvectorBase._N_vector_ 

1019 

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

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

1022 x = x1.cross(x2) 

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

1024 raise IntersectionError(_colinear_) 

1025 a, b = x.philam 

1026 # choose intersection similar to sphericalNvector 

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

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

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

1030 

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

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

1033 intersection, LatLon, LatLon_kwds) 

1034 

1035 

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

1037 LatLon=LatLon, **LatLon_kwds): 

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

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

1040 

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

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

1043 the initial bearing at the first start point 

1044 (compass C{degrees360}). 

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

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

1047 the initial bearing at the second start point 

1048 (compass C{degrees360}). 

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

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

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

1052 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1056 

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

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

1059 height)}. An alternate intersection point might 

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

1061 

1062 @raise IntersectionError: Ambiguous or infinite intersection 

1063 or colinear, parallel or otherwise 

1064 non-intersecting lines. 

1065 

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

1067 

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

1069 

1070 @example: 

1071 

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

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

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

1075 ''' 

1076 s1 = _T00.others(start1=start1) 

1077 s2 = _T00.others(start2=start2) 

1078 try: 

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

1080 LatLon=LatLon, **LatLon_kwds) 

1081 except (TypeError, ValueError) as x: 

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

1083 

1084 

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

1086 height=None, wrap=True, 

1087 LatLon=LatLon, **LatLon_kwds): 

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

1089 by a center point and a radius. 

1090 

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

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

1093 see B{C{radius}}). 

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

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

1096 see B{C{radius}}). 

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

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

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

1100 B{C{radius}}). 

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

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

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

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

1105 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1109 

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

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

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

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

1114 

1115 @raise IntersectionError: Concentric, antipodal, invalid or 

1116 non-intersecting circles. 

1117 

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

1119 

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

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

1122 

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

1124 

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

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

1127 ''' 

1128 c1 = _T00.others(center1=center1) 

1129 c2 = _T00.others(center2=center2) 

1130 try: 

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

1132 height=height, wrap=wrap, 

1133 LatLon=LatLon, **LatLon_kwds) 

1134 except (TypeError, ValueError) as x: 

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

1136 

1137 

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

1139 height=None, too_d=None, wrap=True, 

1140 LatLon=LatLon, **LatLon_kwds): 

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

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

1143 

1144 def _dest3(bearing, h): 

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

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

1147 intersections2, LatLon, LatLon_kwds) 

1148 

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

1150 if f: # swapped 

1151 c1, c2 = c2, c1 # PYCHOK swap 

1152 

1153 a1, b1 = c1.philam 

1154 a2, b2 = c2.philam 

1155 

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

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

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

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

1160 

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

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

1163 if r < _0_0: 

1164 raise _ValueError(eps=r) 

1165 

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

1167 if x > max(r, EPS): 

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

1169 x = sd * sr1 

1170 if isnear0(x): 

1171 raise _ValueError(_invalid_) 

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

1173 

1174 elif x < r: # PYCHOK no cover 

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

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

1177 

1178 if height is None: # "radical height" 

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

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

1181 else: 

1182 h = Height(height) 

1183 

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

1185 if x < EPS4: # externally ... 

1186 r = _dest3(b, h) 

1187 elif x > _PI_EPS4: # internally ... 

1188 r = _dest3(b + PI, h) 

1189 else: 

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

1191 return r, r # ... abutting circles 

1192 

1193 

1194@deprecated_function 

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

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

1197 ''' 

1198 return ispolar(points, wrap=wrap) 

1199 

1200 

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

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

1203 ''' 

1204 n = func.__name__ 

1205 if LatLon is None: 

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

1207 else: 

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

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

1210 return r 

1211 

1212 

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

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

1215 

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

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

1218 height (C{meter}). 

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

1220 or C{None}. 

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

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

1223 

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

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

1226 is C{None}. 

1227 

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

1229 

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

1231 ''' 

1232 # geographic mean 

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

1234 lat, lon = m._N_vector.latlon 

1235 

1236 if height is None: 

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

1238 else: # PYCHOK no cover 

1239 h = Height(height) 

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

1241 

1242 

1243@deprecated_function 

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

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

1246 

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

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

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

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

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

1252 ''' 

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

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

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

1256 return ll, d 

1257 

1258 

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

1260 LatLon=LatLon, **options): 

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

1262 

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

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

1265 

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

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

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

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

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

1271 or C{None}. 

1272 @kwarg options: Optional keyword arguments for function 

1273 L{pygeodesy.equirectangular_}. 

1274 

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

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

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

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

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

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

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

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

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

1284 

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

1286 see function L{pygeodesy.equirectangular_}. 

1287 

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

1289 

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

1291 

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

1293 

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

1295 ''' 

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

1297 LatLon=None, **options) 

1298 d = degrees2m(d, radius=radius) 

1299 n = nearestOn3.__name__ 

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

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

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

1303 

1304 

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

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

1307 (with great circle arcs joining the points). 

1308 

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

1310 or L{BooleanGH}). 

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

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

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

1314 

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

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

1317 

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

1319 

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

1321 

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

1323 C{B{points}} a composite. 

1324 

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

1326 

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

1328 and L{ellipsoidalKarney.perimeterOf}. 

1329 ''' 

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

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

1332 a1, b1 = Ps[0].philam 

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

1334 a2, b2 = p.philam 

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

1336 yield vincentys_(a2, a1, db) 

1337 a1, b1 = a2, b2 

1338 

1339 if _MODS.booleans.isBoolean(points): 

1340 if not closed: 

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

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

1343 else: 

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

1345 return _r2m(r, radius) 

1346 

1347 

1348def _r2m(r, radius): 

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

1350 ''' 

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

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

1353 return r 

1354 

1355 

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

1357 excess=excessAbc_, 

1358 wrap=False): 

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

1360 

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

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

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

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

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

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

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

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

1369 or C{None}. 

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

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

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

1373 

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

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

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

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

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

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

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

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

1382 ''' 

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

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

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

1386 excess=excess, wrap=wrap) 

1387 return _t7Tuple(t, radius) 

1388 

1389 

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

1391 wrap=False): 

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

1393 excess} of a (spherical) triangle. 

1394 

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

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

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

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

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

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

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

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

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

1404 

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

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

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

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

1409 ''' 

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

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

1412 a = vincentys_(phiC, phiB, d) 

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

1414 

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

1416 s = sb * sc 

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

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

1419 

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

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

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

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

1424 

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

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

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

1428 

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

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

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

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

1433 

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

1435 

1436 

1437def _t7Tuple(t, radius): 

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

1439 ''' 

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

1441 r = radius if isscalar(radius) else \ 

1442 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1445 B, (r * t.b), 

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

1447 return t 

1448 

1449 

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

1451 areaOf, # functions 

1452 intersecant2, intersection, intersections2, ispolar, 

1453 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1454 meanOf, 

1455 nearestOn2, nearestOn3, 

1456 perimeterOf, 

1457 sumOf, # == vector3d.sumOf 

1458 triangle7, triangle8_) 

1459 

1460# **) MIT License 

1461# 

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

1463# 

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

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

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

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

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

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

1470# 

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

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

1473# 

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

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

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

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

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

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

1480# OTHER DEALINGS IN THE SOFTWARE.