Coverage for pygeodesy/sphericalTrigonometry.py: 93%

393 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-02-05 16:22 -0500

1 

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

3 

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

5 

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

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

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

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

10 

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

12spherical trigonometry, transcoded from JavaScript originals by 

13I{(C) Chris Veness 2011-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, map1, signOf 

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

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

22 _1_0, _2_0, _90_0 

23from pygeodesy.datums import _ellipsoidal_datum, _mean_radius 

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

25 _TypeError, _ValueError, IntersectionError, \ 

26 _xError, _xkwds, _xkwds_get, _xkwds_pop 

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

28from pygeodesy.fsums import Fsum, fsum, fsumf_ 

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

30 excessGirard_, excessLHuilier_, opposing_, _radical2, \ 

31 vincentys_ 

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

33 _concentric_, _convex_, _end_, _infinite_, _invalid_,\ 

34 _line_, _near_, _not_, _null_, _parallel_, _point_, \ 

35 _SPACE_, _too_ 

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

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

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

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

40 Triangle7Tuple, Triangle8Tuple 

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

42 Fmt as _Fmt, notImplemented # XXX shadowed 

43from pygeodesy.props import deprecated_function, deprecated_method 

44from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \ 

45 _intersecant2, LatLonSphericalBase, \ 

46 _rads3, _radians2m, _trilaterate5 

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

48from pygeodesy.units import Bearing_, Height, _isDegrees, _isRadius, Lam_, \ 

49 Phi_, Radius_, Scalar 

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

51 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \ 

52 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI 

53from pygeodesy.vector3d import sumOf, Vector3d 

54 

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

56 

57__all__ = _ALL_LAZY.sphericalTrigonometry 

58__version__ = '23.12.18' 

59 

60_PI_EPS4 = PI - EPS4 

61if _PI_EPS4 >= PI: 

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

63 

64 

65class Cartesian(CartesianSphericalBase): 

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

67 spherical, geodetic L{LatLon}. 

68 ''' 

69 

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

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

72 

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

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

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

76 

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

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

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

80 

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

82 ''' 

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

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

85 

86 

87class LatLon(LatLonSphericalBase): 

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

89 ''' 

90 

91 def _ab1_ab2_db5(self, other, wrap): 

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

93 ''' 

94 a1, b1 = self.philam 

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

96 if wrap: 

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

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

99 else: # unrollPI shortcut 

100 db = b2 - b1 

101 return a1, b1, a2, b2, db 

102 

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

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

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

106 end point. 

107 

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

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

110 from the start point to the point where the perpendicular 

111 crosses the line. 

112 

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

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

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

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

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

118 

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

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

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

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

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

124 B{C{start}} point. 

125 

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

127 

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

129 ''' 

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

131 cx = cos(x) 

132 return _0_0 if isnear0(cx) else \ 

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

134 

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

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

137 ''' 

138 s = self.others(start=start) 

139 e = self.others(end=end) 

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

141 

142 r = Radius_(radius) 

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

144 

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

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

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

148 return r, x, -b 

149 

150 @deprecated_method 

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

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

153 ''' 

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

155 

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

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

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

159 

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

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

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

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

164 

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

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

167 ''' 

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

169 sa, ca, sa1, ca1, \ 

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

171 sa1 *= ca2 * ca 

172 

173 x = sa1 * sdb 

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

175 z = ca1 * sdb * ca2 * sa 

176 

177 h = hypot(x, y) 

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

179 return None # great circle doesn't reach latitude 

180 

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

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

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

184 

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

186 '''Compute the (signed) distance from this point to 

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

188 

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

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

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

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

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

194 

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

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

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

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

199 

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

201 

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

203 ''' 

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

205 return _radians2m(x, radius) 

206 

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

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

209 travelled the given distance on the given initial bearing. 

210 

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

212 B{C{radius}}). 

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

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

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

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

217 

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

219 

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

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

222 ''' 

223 a, b = self.philam 

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

225 

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

227 h = self._heigHt(height) 

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

229 

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

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

232 

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

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

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

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

237 

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

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

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

241 

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

243 

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

245 ''' 

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

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

248 

249# @Property_RO 

250# def Ecef(self): 

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

252# ''' 

253# return _MODS.ecef.EcefKarney 

254 

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

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

257 on the given initial bearing from this point. 

258 

259 Direction of vector is such that initial bearing vector 

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

261 

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

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

264 overriding the default L{Vector3d}. 

265 @kwarg Vector_kwds: Optional, additional keyword argunents 

266 for B{C{Vector}}. 

267 

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

269 

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

271 ''' 

272 a, b = self.philam 

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

274 

275 return Vector(sb * ct - cb * sa * st, 

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

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

278 

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

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

281 to an other point. 

282 

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

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

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

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

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

288 C{sphericalNvector.LatLon.initialBearingTo}. 

289 

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

291 

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

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

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

295 

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

297 ''' 

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

299 # XXX behavior like sphericalNvector.LatLon.initialBearingTo 

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

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

302 

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

304 

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

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

307 and an other point. 

308 

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

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

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

312 @kwarg height: Optional height, overriding the intermediate 

313 height (C{meter}). 

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

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

316 

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

318 

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

320 

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

322 

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

324 ''' 

325 p = self 

326 f = Scalar(fraction=fraction) 

327 if not isnear0(f): 

328 p = p.others(other) 

329 if wrap: 

330 p = _Wrap.point(p) 

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

332 a1, b1 = self.philam 

333 a2, b2 = p.philam 

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

335 r = vincentys_(a2, a1, db) 

336 sr = sin(r) 

337 if isnon0(sr): 

338 sa1, ca1, sa2, ca2, \ 

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

340 

341 t = f * r 

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

343 b = sin( t) # / sr superflous 

344 

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

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

347 z = a * sa1 + b * sa2 

348 

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

350 b = atan2d(y, x) 

351 

352 else: # PYCHOK no cover 

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

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

355 

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

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

358 return p 

359 

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

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

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

363 

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

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

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

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

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

369 C{degrees360}). 

370 @kwarg height: Optional height for intersection point, 

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

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

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

374 

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

376 intersection point might be the L{antipode} to 

377 the returned result. 

378 

379 @raise IntersectionError: Ambiguous or infinite intersection 

380 or colinear, parallel or otherwise 

381 non-intersecting lines. 

382 

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

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

385 

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

387 ''' 

388 try: 

389 s2 = self.others(other) 

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

391 LatLon=self.classof) 

392 except (TypeError, ValueError) as x: 

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

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

395 

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

397 height=None, wrap=True): 

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

399 by a center point and radius. 

400 

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

402 see B{C{radius}}). 

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

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

405 see B{C{radius}}). 

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

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

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

409 B{C{radius}}). 

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

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

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

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

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

415 

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

417 instance. For abutting circles, both intersection 

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

419 

420 @raise IntersectionError: Concentric, antipodal, invalid or 

421 non-intersecting circles. 

422 

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

424 

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

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

427 ''' 

428 try: 

429 c2 = self.others(other) 

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

431 height=height, wrap=wrap, 

432 LatLon=self.classof) 

433 except (TypeError, ValueError) as x: 

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

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

436 

437 @deprecated_method 

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

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

440 return self.isenclosedBy(points) 

441 

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

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

444 

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

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

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

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

449 

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

451 composite, C{False} otherwise. 

452 

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

454 

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

456 

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

458 

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

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

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

462 ''' 

463 if _MODS.booleans.isBoolean(points): 

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

465 

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

467 n0 = self._N_vector 

468 

469 v2 = Ps[0]._N_vector 

470 p1 = Ps[1] 

471 v1 = p1._N_vector 

472 # check whether this point on same side of all 

473 # polygon edges (to the left or right depending 

474 # on the anti-/clockwise polygon direction) 

475 gc1 = v2.cross(v1) 

476 t0 = gc1.angleTo(n0) > PI_2 

477 s0 = None 

478 # get great-circle vector for each edge 

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

480 if wrap and not Ps.looped: 

481 p2 = _unrollon(p1, p2) 

482 p1 = p2 

483 v2 = p2._N_vector 

484 gc = v1.cross(v2) 

485 t = gc.angleTo(n0) > PI_2 

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

487 return False # outside 

488 

489 # check for convex polygon: angle between 

490 # gc vectors, signed by direction of n0 

491 # (otherwise the test above is not reliable) 

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

493 if s != s0: 

494 if s0 is None: 

495 s0 = s 

496 else: 

497 t = _Fmt.SQUARE(points=i) 

498 raise _ValueError(t, p2, wrap=wrap, txt=_not_(_convex_)) 

499 gc1, v1 = gc, v2 

500 

501 return True # inside 

502 

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

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

505 

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

507 @kwarg height: Optional height for midpoint, overriding 

508 the mean height (C{meter}). 

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

510 may be negative or greater than 1.0. 

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

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

513 

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

515 

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

517 

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

519 

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

521 ''' 

522 if fraction is _0_5: 

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

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

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

526 

527 x = ca2 * cdb + ca1 

528 y = ca2 * sdb 

529 

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

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

532 

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

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

535 else: 

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

537 return r 

538 

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

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

541 

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

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

544 

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

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

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

548 @kwarg wrap_adjust_limit: Optional keyword arguments for functions 

549 L{sphericalTrigonometry.nearestOn3} and 

550 L{pygeodesy.equirectangular_}, 

551 

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

553 

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

555 see function L{pygeodesy.equirectangular_}. 

556 

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

558 is not (yet) supported. 

559 

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

561 

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

563 

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

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

566 ''' 

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

568 w = _xkwds_pop(wrap_adjust_limit, within=True) 

569 if not w: 

570 notImplemented(self, within=w) 

571 

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

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

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

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

576# return point1 

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

578# if isnear0(d): 

579# return point1 # or point2 

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

581# return point2 

582# f = a / d 

583# if within: 

584# if f > EPS1: 

585# return point2 

586# elif f < EPS: 

587# return point1 

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

589 

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

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

592 **wrap_adjust_limit)[0] 

593 

594 @deprecated_method 

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

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

597 

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

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

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

601 units of B{C{radius}}. 

602 ''' 

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

604 return r.closest, r.distance 

605 

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

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

608 

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

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

611 

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

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

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

615 @kwarg wrap_adjust_limit: Optional keyword arguments for function 

616 L{sphericalTrigonometry.nearestOn3} and 

617 L{pygeodesy.equirectangular_}, 

618 

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

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

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

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

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

624 function L{pygeodesy.compassAngle}. 

625 

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

627 see function L{pygeodesy.equirectangular_}. 

628 

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

630 

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

632 

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

634 

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

636 and L{pygeodesy.nearestOn5}. 

637 ''' 

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

639 LatLon=self.classof, **wrap_adjust_limit) 

640 

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

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

643 coordinates. 

644 

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

646 and other keyword arguments, ignored 

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

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

649 this L{Cartesian} class or specify 

650 C{B{Cartesian}=None}. 

651 

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

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

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

655 

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

657 ''' 

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

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

660 

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

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

663 

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

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

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

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

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

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

670 B{C{otherB}} and B{C{otherC}} points (C{bool}). 

671 

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

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

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

675 

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

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

678 ''' 

679 B = self.others(otherB=otherB) 

680 C = self.others(otherC=otherC) 

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

682 

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

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

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

686 

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

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

689 at this and at the other point. 

690 

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

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

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

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

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

696 

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

698 

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

700 ''' 

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

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

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

704 

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

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

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

708 intersection} of three corresponding circles. 

709 

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

711 as B{C{radius}}). 

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

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

714 B{C{radius}}). 

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

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

717 B{C{radius}}). 

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

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

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

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

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

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

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

725 B{C{point2}} and B{C{point3}} (C{bool}). 

726 

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

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

729 the corresponding trilaterated points C{minPoint} and 

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

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

732 

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

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

735 

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

737 respectively largest I{radial} overlap found. 

738 

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

740 nearest respectively farthest intersection margin. 

741 

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

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

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

745 largest B{C{distance#}}. 

746 

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

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

749 no intersection or all (near-)concentric for 

750 C{B{area}=False}. 

751 

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

753 

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

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

756 or B{C{radius}}. 

757 ''' 

758 return _trilaterate5(self, distance1, 

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

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

761 area=area, radius=radius, eps=eps, wrap=wrap) 

762 

763 

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

765 

766 

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

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

769 (with the pointsjoined by great circle arcs). 

770 

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

772 or L{BooleanGH}). 

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

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

775 or C{None}. 

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

777 (C{bool}). 

778 

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

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

781 

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

783 

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

785 

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

787 

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

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

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

791 

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

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

794 L{ellipsoidalKarney.areaOf}. 

795 ''' 

796 if _MODS.booleans.isBoolean(points): 

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

798 

799 _at2, _t_2 = atan2, tan_2 

800 _un, _w180 = unrollPI, wrap180 

801 

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

803 p1 = p2 = Ps[0] 

804 a1, b1 = p1.philam 

805 ta1, z1 = _t_2(a1), None 

806 

807 A = Fsum() # mean phi 

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

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

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

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

812 D = Fsum() 

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

814 a2, b2 = p2.philam 

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

816 A += a2 

817 ta2 = _t_2(a2) 

818 tdb = _t_2(db, points=i) 

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

820 _1_0 + ta1 * ta2) 

821 ta1, b1 = ta2, b2 

822 

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

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

825 if z1 is not None: 

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

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

828 p1, z1 = p2, z2 

829 

830 R = abs(R * _2_0) 

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

832 R = abs(R - PI2) 

833 if radius: 

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

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

836 return float(R) 

837 

838 

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

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

841 

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

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

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

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

846 

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

848 ''' 

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

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

851 ca *= sr 

852 

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

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

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

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

857 

858 

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

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

861 # and similar logic in .ellipsoidalBaseDI._intersect3 

862 a1, b1 = s.philam 

863 

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

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

866 else: # must be a point 

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

868 hs.append(end.height) 

869 a2, b2 = end.philam 

870 if wrap: 

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

872 

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

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

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

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

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

878 -(b1 + b2) * _0_5) 

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

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

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

882 cb12 * cb21 + sb12 * sb21, 

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

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

885 

886 

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

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

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

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

891 

892 

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

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

895 two points or as a point and bearing. 

896 

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

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

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

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

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

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

903 @kwarg radius_exact_height_wrap: Optional keyword arguments, see 

904 method L{LatLon.intersecant2} for further details. 

905 

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

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

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

909 

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

911 

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

913 B{C{circle}} or B{C{other}} invalid. 

914 

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

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

917 ''' 

918 c = _T00.others(center=center) 

919 p = _T00.others(point=point) 

920 try: 

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

922 except (TypeError, ValueError) as x: 

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

924 **radius_exact_height_wrap) 

925 

926 

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

928 LatLon=None, **LatLon_kwds): 

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

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

931 

932 s1, s2 = start1, start2 

933 if wrap: 

934 s2 = _Wrap.point(s2) 

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

936 

937 a1, b1 = s1.philam 

938 a2, b2 = s2.philam 

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

940 r12 = vincentys_(a2, a1, db) 

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

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

943 

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

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

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

947 

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

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

950 raise IntersectionError(_parallel_) 

951 # handle domain error for equivalent longitudes, 

952 # see also functions asin_safe and acos_safe at 

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

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

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

956 if sin(db) > 0: 

957 t21 = PI2 - t21 

958 else: 

959 t12 = PI2 - t12 

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

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

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

963 raise IntersectionError(_infinite_) 

964 sx3 = sx1 * sx2 

965# XXX if sx3 < 0: 

966# XXX raise ValueError(_ambiguous_) 

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

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

969 

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

971 # like .ellipsoidalBaseDI,_intersect3, if this intersection 

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

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

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

975 

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

977 _N_vector_ = _MODS.nvectorBase._N_vector_ 

978 

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

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

981 x = x1.cross(x2) 

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

983 raise IntersectionError(_colinear_) 

984 a, b = x.philam 

985 # choose intersection similar to sphericalNvector 

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

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

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

989 

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

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

992 intersection, LatLon, LatLon_kwds) 

993 

994 

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

996 LatLon=LatLon, **LatLon_kwds): 

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

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

999 

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

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

1002 the initial bearing at the first start point 

1003 (compass C{degrees360}). 

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

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

1006 the initial bearing at the second start point 

1007 (compass C{degrees360}). 

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

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

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

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

1012 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1016 

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

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

1019 height)}. An alternate intersection point might 

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

1021 

1022 @raise IntersectionError: Ambiguous or infinite intersection 

1023 or colinear, parallel or otherwise 

1024 non-intersecting lines. 

1025 

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

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

1028 

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

1030 ''' 

1031 s1 = _T00.others(start1=start1) 

1032 s2 = _T00.others(start2=start2) 

1033 try: 

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

1035 LatLon=LatLon, **LatLon_kwds) 

1036 except (TypeError, ValueError) as x: 

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

1038 

1039 

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

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

1042 LatLon=LatLon, **LatLon_kwds): 

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

1044 by a center point and a radius. 

1045 

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

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

1048 see B{C{radius}}). 

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

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

1051 see B{C{radius}}). 

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

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

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

1055 B{C{radius}}). 

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

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

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

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

1060 (C{bool}). 

1061 @kwarg LatLon: Optional class to return the intersection 

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

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

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

1065 

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

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

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

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

1070 

1071 @raise IntersectionError: Concentric, antipodal, invalid or 

1072 non-intersecting circles. 

1073 

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

1075 

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

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

1078 

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

1080 

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

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

1083 ''' 

1084 c1 = _T00.others(center1=center1) 

1085 c2 = _T00.others(center2=center2) 

1086 try: 

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

1088 height=height, wrap=wrap, 

1089 LatLon=LatLon, **LatLon_kwds) 

1090 except (TypeError, ValueError) as x: 

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

1092 center2=center2, rad2=rad2, wrap=wrap) 

1093 

1094 

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

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

1097 LatLon=LatLon, **LatLon_kwds): 

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

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

1100 

1101 def _dest3(bearing, h): 

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

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

1104 intersections2, LatLon, LatLon_kwds) 

1105 

1106 a1, b1 = c1.philam 

1107 a2, b2 = c2.philam 

1108 if wrap: 

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

1110 

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

1112 if f: # swapped radii, swap centers 

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

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

1115 

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

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

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

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

1120 

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

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

1123 if r < _0_0: 

1124 raise _ValueError(eps=r) 

1125 

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

1127 if x > max(r, EPS): 

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

1129 x = sd * sr1 

1130 if isnear0(x): 

1131 raise _ValueError(_invalid_) 

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

1133 

1134 elif x < r: # PYCHOK no cover 

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

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

1137 

1138 if height is None: # "radical height" 

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

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

1141 else: 

1142 h = Height(height) 

1143 

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

1145 if x < EPS4: # externally ... 

1146 r = _dest3(b, h) 

1147 elif x > _PI_EPS4: # internally ... 

1148 r = _dest3(b + PI, h) 

1149 else: 

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

1151 return r, r # ... abutting circles 

1152 

1153 

1154@deprecated_function 

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

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

1157 ''' 

1158 return ispolar(points, wrap=wrap) 

1159 

1160 

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

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

1163 ''' 

1164 n = where.__name__ 

1165 if LatLon is None: 

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

1167 else: 

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

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

1170 return r 

1171 

1172 

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

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

1175 

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

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

1178 height (C{meter}). 

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

1180 (C{bool}). 

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

1182 or C{None}. 

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

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

1185 

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

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

1188 is C{None}. 

1189 

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

1191 

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

1193 ''' 

1194 def _N_vs(ps, w): 

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

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

1197 yield p._N_vector 

1198 

1199 m = _MODS.nvectorBase 

1200 # geographic, vectorial mean 

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

1202 lat, lon, h = n.latlonheight 

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

1204 

1205 

1206@deprecated_function 

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

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

1209 

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

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

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

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

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

1215 ''' 

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

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

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

1219 return ll, d 

1220 

1221 

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

1223 limit=9, **LatLon_and_kwds): 

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

1225 

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

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

1228 

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

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

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

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

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

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

1235 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}). 

1236 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}), 

1237 default C{9 degrees} is about C{1,000 Kmeter} (for mean 

1238 spherical earth radius L{R_KM}). 

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

1240 or C{None}. 

1241 @kwarg options: Optional keyword arguments for function 

1242 L{pygeodesy.equirectangular_}. 

1243 

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

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

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

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

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

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

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

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

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

1253 

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

1255 see function L{pygeodesy.equirectangular_}. 

1256 

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

1258 

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

1260 

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

1262 

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

1264 ''' 

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

1266 adjust=adjust, limit=limit) 

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

1268 h = t.height 

1269 n = nearestOn3.__name__ 

1270 

1271 kwds = _xkwds(LatLon_and_kwds, height=h, name=n) 

1272 LL = _xkwds_pop(kwds, LatLon=LatLon) 

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

1274 LL(t.lat, t.lon, **kwds) 

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

1276 

1277 

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

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

1280 (with great circle arcs joining the points). 

1281 

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

1283 or L{BooleanGH}). 

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

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

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

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

1288 

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

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

1291 

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

1293 

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

1295 

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

1297 C{B{points}} a composite. 

1298 

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

1300 

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

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

1303 ''' 

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

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

1306 a1, b1 = Ps[0].philam 

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

1308 a2, b2 = p.philam 

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

1310 yield vincentys_(a2, a1, db) 

1311 a1, b1 = a2, b2 

1312 

1313 if _MODS.booleans.isBoolean(points): 

1314 if not closed: 

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

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

1317 else: 

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

1319 return _radians2m(r, radius) 

1320 

1321 

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

1323 excess=excessAbc_, 

1324 wrap=False): 

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

1326 

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

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

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

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

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

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

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

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

1335 or C{None}. 

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

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

1338 @kwarg wrap: If C{True}, wrap and L{pygeodesy.unroll180} 

1339 longitudes (C{bool}). 

1340 

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

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

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

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

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

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

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

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

1349 ''' 

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

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

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

1353 excess=excess, wrap=wrap) 

1354 return _t7Tuple(t, radius) 

1355 

1356 

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

1358 wrap=False): 

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

1360 excess} of a (spherical) triangle. 

1361 

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

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

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

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

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

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

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

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

1370 @kwarg wrap: If C{True}, L{pygeodesy.unrollPI} the 

1371 longitudinal deltas (C{bool}). 

1372 

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

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

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

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

1377 ''' 

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

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

1380 a = vincentys_(phiC, phiB, d) 

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

1382 

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

1384 s = sb * sc 

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

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

1387 

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

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

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

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

1392 

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

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

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

1396 

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

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

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

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

1401 

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

1403 

1404 

1405def _t7Tuple(t, radius): 

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

1407 ''' 

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

1409 r = radius if _isRadius(radius) else \ 

1410 _ellipsoidal_datum(radius).ellipsoid.Rmean 

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

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

1413 B, (r * t.b), 

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

1415 return t 

1416 

1417 

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

1419 areaOf, # functions 

1420 intersecant2, intersection, intersections2, ispolar, 

1421 isPoleEnclosedBy, # DEPRECATED, use ispolar 

1422 meanOf, 

1423 nearestOn2, nearestOn3, 

1424 perimeterOf, 

1425 sumOf, # XXX == vector3d.sumOf 

1426 triangle7, triangle8_) 

1427 

1428# **) MIT License 

1429# 

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

1431# 

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

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

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

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

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

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

1438# 

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

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

1441# 

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

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

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

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

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

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

1448# OTHER DEALINGS IN THE SOFTWARE.