1

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

3

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

5

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

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

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

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

10

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

12 spherical trigonometry, transcoded from JavaScript originals by

13 I{(C) Chris Veness 2011-2024} published under the same MIT Licence**, see

14 U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}.

15'''

16# make sure int/int division yields float quotient, see .basics

17 from __future__ import division as _; del _ # PYCHOK semicolon

18

19 from pygeodesy.basics import copysign0, map1, signOf

20 from 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

23 from pygeodesy.datums import _ellipsoidal_datum, _mean_radius

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

25 _TypeError, _ValueError, IntersectionError, \

26 _xError, _xkwds, _xkwds_get, _xkwds_pop2

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

28 from pygeodesy.fsums import Fsum, fsum, fsumf_

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

30 excessGirard_, excessLHuilier_, opposing_, _radical2, \

31 vincentys_

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

33 _concentric_, _convex_, _end_, _infinite_, \

34 _invalid_, _line_, _near_, _null_, _parallel_, \

35 _point_, _SPACE_, _too_

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

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

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

39 Triangle7Tuple, Triangle8Tuple

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

41 Fmt as _Fmt # XXX shadowed

42 from pygeodesy.props import deprecated_function, deprecated_method

43 from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \

44 _intersecant2, LatLonSphericalBase, \

45 _rads3, _radians2m, _trilaterate5

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

47 from pygeodesy.units import Bearing_, Height, _isDegrees, _isRadius, Lamd, \

48 Phid, Radius_, Scalar

49 from pygeodesy.utily import acos1, asin1, atan1d, atan2, atan2d, degrees90, \

50 degrees180, degrees2m, m2radians, radiansPI2, \

51 sincos2_, tan_2, unrollPI, _unrollon, _unrollon3, \

52 wrap180, wrapPI, _Wrap

53 from pygeodesy.vector3d import sumOf, Vector3d

54

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

56

57__all__ = _ALL_LAZY.sphericalTrigonometry

58__version__ = '24.11.24'

59

60_PI_EPS4 = PI - EPS4

61 if _PI_EPS4 >= PI:

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

63

64

65 class 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 C{B{LatLon} is 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

87 class 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 a great

187 circle from a start to 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 a bearing from North.

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 C{B{radius} is 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 from this point on the bearing from North.

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 sa *= st

276 return Vector(fdot_(sb, ct, -cb, sa),

277 -fdot_(cb, ct, sb, sa),

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

279

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

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

282 to an other point.

283

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

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

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

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

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

289 C{sphericalNvector.LatLon.initialBearingTo}.

290

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

292

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

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

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

296

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

298 '''

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

300 # XXX behavior like sphericalNvector.LatLon.initialBearingTo

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

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

303

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

305

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

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

308 and an other point.

309

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

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

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

313 @kwarg height: Optional height, overriding the intermediate

314 height (C{meter}).

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

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

317

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

319

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

321

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

323

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

325 '''

326 p = self

327 f = Scalar(fraction=fraction)

328 if not isnear0(f):

329 p = p.others(other)

330 if wrap:

331 p = _Wrap.point(p)

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

333 a1, b1 = self.philam

334 a2, b2 = p.philam

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

336 r = vincentys_(a2, a1, db)

337 sr = sin(r)

338 if isnon0(sr):

339 sa1, ca1, sa2, ca2, \

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

341

342 t = f * r

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

344 b = sin( t) # / sr superflous

345

346 x = fdot_(a, ca1 * cb1, b, ca2 * cb2)

347 y = fdot_(a, ca1 * sb1, b, ca2 * sb2)

348 z = fdot_(a, sa1, b, sa2)

349

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

351 b = atan2d(y, x)

352

353 else: # PYCHOK no cover

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

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

356

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

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

359 return p

360

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

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

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

364

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

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

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

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

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

370 C{degrees360}).

371 @kwarg height: Optional height for intersection point,

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

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

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

375

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

377 intersection point might be the L{antipode} to

378 the returned result.

379

380 @raise IntersectionError: Ambiguous or infinite intersection

381 or colinear, parallel or otherwise

382 non-intersecting lines.

383

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

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

386

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

388 '''

389 try:

390 s2 = self.others(other)

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

392 LatLon=self.classof)

393 except (TypeError, ValueError) as x:

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

395 other=other, end2=end2, wrap=wrap)

396

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

398 height=None, wrap=True):

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

400 by a center point and a radius.

401

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

403 see B{C{radius}}).

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

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

406 see B{C{radius}}).

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

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

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

410 B{C{radius}}).

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

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

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

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

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

416

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

418 instance. For abutting circles, both intersection

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

420

421 @raise IntersectionError: Concentric, antipodal, invalid or

422 non-intersecting circles.

423

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

425

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

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

428 '''

429 try:

430 c2 = self.others(other)

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

432 height=height, wrap=wrap,

433 LatLon=self.classof)

434 except (TypeError, ValueError) as x:

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

436 other=other, rad2=rad2, wrap=wrap)

437

438 @deprecated_method

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

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

441 return self.isenclosedBy(points)

442

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

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

445

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

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

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

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

450

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

452 composite, C{False} otherwise.

453

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

455

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

457

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

459

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

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

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

463 '''

464 if _MODS.booleans.isBoolean(points):

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

466

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

468 n0 = self._N_vector

469

470 v2 = Ps[0]._N_vector

471 p1 = Ps[1]

472 v1 = p1._N_vector

473 # check whether this point on same side of all

474 # polygon edges (to the left or right depending

475 # on the anti-/clockwise polygon direction)

476 gc1 = v2.cross(v1)

477 t0 = gc1.angleTo(n0) > PI_2

478 s0 = None

479 # get great-circle vector for each edge

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

481 if wrap and not Ps.looped:

482 p2 = _unrollon(p1, p2)

483 p1 = p2

484 v2 = p2._N_vector

485 gc = v1.cross(v2)

486 t = gc.angleTo(n0) > PI_2

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

488 return False # outside

489

490 # check for convex polygon: angle between

491 # gc vectors, signed by direction of n0

492 # (otherwise the test above is not reliable)

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

494 if s != s0:

495 if s0 is None:

496 s0 = s

497 else:

498 t = _Fmt.SQUARE(points=i)

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

500 gc1, v1 = gc, v2

501

502 return True # inside

503

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

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

506

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

508 @kwarg height: Optional height for midpoint, overriding

509 the mean height (C{meter}).

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

511 may be negative or greater than 1.0.

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

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

514

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

516

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

518

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

520

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

522 '''

523 if fraction is _0_5:

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

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

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

527

528 x = ca2 * cdb + ca1

529 y = ca2 * sdb

530

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

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

533

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

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

536 else:

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

538 return r

539

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

541 '''Locate the point between two other points closest to this point.

542

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

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

545

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

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

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

549 @kwarg wrap_adjust_limit: Optional keyword arguments for functions

550 L{sphericalTrigonometry.nearestOn3} and

551 L{pygeodesy.equirectangular4},

552

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

554

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

556 see function L{pygeodesy.equirectangular4}.

557

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

559 is not (yet) supported.

560

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

562

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

564

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

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

567 '''

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

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

570 if not w:

571 self._notImplemented(within=w)

572

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

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

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

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

577# return point1

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

579# if isnear0(d):

580# return point1 # or point2

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

582# return point2

583# f = a / d

584# if within:

585# if f > EPS1:

586# return point2

587# elif f < EPS:

588# return point1

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

590

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

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

593 **kwds)[0]

594

595 @deprecated_method

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

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

598

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

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

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

602 units of B{C{radius}}.

603 '''

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

605 return r.closest, r.distance

606

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

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

609

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

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

612

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

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

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

616 @kwarg wrap_adjust_limit: Optional keyword arguments for function

617 L{sphericalTrigonometry.nearestOn3} and

618 L{pygeodesy.equirectangular4},

619

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

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

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

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

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

625 function L{pygeodesy.compassAngle}.

626

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

628 see function L{pygeodesy.equirectangular4}.

629

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

631

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

633

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

635

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

637 and L{pygeodesy.nearestOn5}.

638 '''

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

640 LatLon=self.classof, **wrap_adjust_limit)

641

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

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

644

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

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

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

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

649

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

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

652 and C{M} if available.

653

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

655 '''

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

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

658

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

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

661

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

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

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

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

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

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

668

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

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

671

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

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

674 '''

675 B = self.others(otherB=otherB)

676 C = self.others(otherC=otherC)

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

678

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

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

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

682

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

684 '''Locate a point given this, an other point and a bearing from

685 North at both points.

686

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

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

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

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

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

692

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

694

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

696 '''

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

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

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

700

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

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

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

704 of three corresponding circles.

705

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

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

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

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

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

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

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

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

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

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

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

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

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

719

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

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

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

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

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

725

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

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

728

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

730 largest I{radial} overlap found.

731

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

733 respectively farthest intersection margin.

734

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

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

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

738

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

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

741 no intersection or all (near-)concentric if

742 C{B{area}=False}.

743

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

745

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

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

748 or B{C{radius}}.

749 '''

750 return _trilaterate5(self, distance1,

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

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

753 area=area, radius=radius, eps=eps, wrap=wrap)

754

755

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

757

758

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

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

761 points joined by great circle arcs).

762

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

764 or L{BooleanGH}).

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

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

767 or C{None}.

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

769 (C{bool}).

770

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

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

773

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

775

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

777

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

779

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

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

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

783

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

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

786 L{ellipsoidalKarney.areaOf}.

787 '''

788 if _MODS.booleans.isBoolean(points):

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

790

791 _at2, _t_2 = atan2, tan_2

792 _un, _w180 = unrollPI, wrap180

793

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

795 p1 = p2 = Ps[0]

796 a1, b1 = p1.philam

797 ta1, z1 = _t_2(a1), None

798

799 A = Fsum() # mean phi

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

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

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

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

804 D = Fsum()

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

806 a2, b2 = p2.philam

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

808 A += a2

809 ta2 = _t_2(a2)

810 tdb = _t_2(db, points=i)

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

812 _1_0 + ta1 * ta2)

813 ta1, b1 = ta2, b2

814

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

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

817 if z1 is not None:

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

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

820 p1, z1 = p2, z2

821

822 R = abs(R * _2_0)

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

824 R = abs(R - PI2)

825 if radius:

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

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

828 return float(R)

829

830

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

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

833

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

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

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

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

838

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

840 '''

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

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

843 ca *= sr

844

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

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

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

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

849

850

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

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

853 # and similar logic in .ellipsoidalBaseDI._intersect3

854 a1, b1 = s.philam

855

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

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

858 else: # must be a point

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

860 hs.append(end.height)

861 a2, b2 = end.philam

862 if wrap:

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

864

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

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

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

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

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

870 -(b1 + b2) * _0_5)

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

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

873 x = Vector(fdot_(sb12, cb21, -cb12, sb21),

874 fdot_(cb12, cb21, sb12, sb21),

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

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

877

878

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

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

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

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

883

884

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

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

887 two points or as a point and a bearing from North.

888

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

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

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

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

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

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

895 @kwarg radius_exact_height_wrap: Optional keyword arguments, see method

896 L{intersecant2<pygeodesy.sphericalBase.LatLonSphericalBase.

897 intersecant2>} for further details.

898

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

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

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

902

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

904

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

906 not L{LatLon}.

907

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

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

910 '''

911 c = _T00.others(center=center)

912 p = _T00.others(point=point)

913 try:

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

915 except (TypeError, ValueError) as x:

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

917 **radius_exact_height_wrap)

918

919

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

921 LatLon=LatLon, **LatLon_kwds):

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

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

924

925 s1, s2 = start1, start2

926 if wrap:

927 s2 = _Wrap.point(s2)

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

929

930 a1, b1 = s1.philam

931 a2, b2 = s2.philam

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

933 r12 = vincentys_(a2, a1, db)

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

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

936

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

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

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

940

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

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

943 raise IntersectionError(_parallel_)

944 # handle domain error for equivalent longitudes,

945 # see also functions asin_safe and acos_safe at

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

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

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

949 if sin(db) > 0:

950 t21 = PI2 - t21

951 else:

952 t12 = PI2 - t12

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

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

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

956 raise IntersectionError(_infinite_)

957 sx3 = sx1 * sx2

958# XXX if sx3 < 0:

959# XXX raise ValueError(_ambiguous_)

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

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

962

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

964 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

968

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

970 _N_vector_ = _MODS.nvectorBase._N_vector_

971

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

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

974 x = x1.cross(x2)

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

976 raise IntersectionError(_colinear_)

977 a, b = x.philam

978 # choose intersection similar to sphericalNvector

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

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

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

982

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

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

985 intersection, LatLon, LatLon_kwds)

986

987

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

989 **LatLon_and_kwds):

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

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

992

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

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

995 at the first start point (compass C{degrees360}).

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

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

998 at the second start point (compass C{degrees360}).

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

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

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

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

1003 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=}L{LatLon} to use

1004 for the intersection point and optionally additional

1005 B{C{LatLon}} keyword arguments, ignored if C{B{LatLon}

1006 is None}.

1007

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

1009 is None} a L{LatLon3Tuple}C{(lat, lon, height)}. An alternate

1010 intersection point might be the L{antipode} to the returned result.

1011

1012 @raise IntersectionError: Ambiguous or infinite intersection or colinear,

1013 parallel or otherwise non-intersecting lines.

1014

1015 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}} or B{C{end2}}

1016 point not L{LatLon}.

1017

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

1019 '''

1020 s1 = _T00.others(start1=start1)

1021 s2 = _T00.others(start2=start2)

1022 try:

1023 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap, **LatLon_and_kwds)

1024 except (TypeError, ValueError) as x:

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

1026

1027

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

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

1030 **LatLon_and_kwds):

1031 '''Compute the intersection points of two circles each defined by a

1032 center point and a radius.

1033

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

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

1036 B{C{radius}}).

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

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

1039 B{C{radius}}).

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

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

1042 @kwarg eps: Required overlap (C{meter} or C{radians}, see B{C{radius}}).

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

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

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

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

1047 (C{bool}).

1048 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=}L{LatLon} to use for

1049 the intersection points and optionally additional B{C{LatLon}}

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

1051

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

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

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

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

1056

1057 @raise IntersectionError: Concentric, antipodal, invalid or

1058 non-intersecting circles.

1059

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

1061

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

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

1064

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

1066

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

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

1069 '''

1070 c1 = _T00.others(center1=center1)

1071 c2 = _T00.others(center2=center2)

1072 try:

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

1074 height=height, wrap=wrap,

1075 **LatLon_and_kwds)

1076 except (TypeError, ValueError) as x:

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

1078 center2=center2, rad2=rad2, wrap=wrap)

1079

1080

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

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

1083 LatLon=LatLon, **LatLon_kwds):

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

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

1086

1087 def _dest3(bearing, h):

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

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

1090 intersections2, LatLon, LatLon_kwds)

1091

1092 a1, b1 = c1.philam

1093 a2, b2 = c2.philam

1094 if wrap:

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

1096

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

1098 if f: # swapped radii, swap centers

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

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

1101

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

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

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

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

1106

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

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

1109 if r < _0_0:

1110 raise _ValueError(eps=r)

1111

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

1113 if x > max(r, EPS):

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

1115 x = sd * sr1

1116 if isnear0(x):

1117 raise _ValueError(_invalid_)

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

1119

1120 elif x < r: # PYCHOK no cover

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

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

1123

1124 if height is None: # "radical height"

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

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

1127 else:

1128 h = Height(height)

1129

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

1131 if x < EPS4: # externally ...

1132 r = _dest3(b, h)

1133 elif x > _PI_EPS4: # internally ...

1134 r = _dest3(b + PI, h)

1135 else:

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

1137 return r, r # ... abutting circles

1138

1139

1140@deprecated_function

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

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

1143 '''

1144 return ispolar(points, wrap=wrap)

1145

1146

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

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

1149 '''

1150 n = where.__name__

1151 if LatLon is None:

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

1153 else:

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

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

1156 return r

1157

1158

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

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

1161

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

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

1164 (C{meter}).

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

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

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

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

1169

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

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

1172

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

1174

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

1176 '''

1177 def _N_vs(ps, w):

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

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

1180 yield p._N_vector

1181

1182 m = _MODS.nvectorBase

1183 # geographic, vectorial mean

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

1185 lat, lon, h = n.latlonheight

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

1187

1188

1189@deprecated_function

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

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

1192

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

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

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

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

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

1198 '''

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

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

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

1202 return ll, d

1203

1204

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

1206 limit=9, **LatLon_and_kwds):

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

1208

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

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

1211

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

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

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

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

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

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

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

1219 (C{bool}).

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

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

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

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

1224 closest point and optionally additional C{B{LatLon}} keyword

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

1226

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

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

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

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

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

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

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

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

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

1236

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

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

1239

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

1241

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

1243

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

1245

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

1247 L{nearestOn5<pygeodesy.nearestOn5>}.

1248 '''

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

1250 adjust=adjust, limit=limit)

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

1252 h = t.height

1253 n = nearestOn3.__name__

1254

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

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

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

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

1259

1260

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

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

1263 (with great circle arcs joining the points).

1264

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

1266 or L{BooleanGH}).

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

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

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

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

1271

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

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

1274

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

1276

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

1278

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

1280 C{B{points}} a composite.

1281

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

1283

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

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

1286 '''

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

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

1289 a1, b1 = Ps[0].philam

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

1291 a2, b2 = p.philam

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

1293 yield vincentys_(a2, a1, db)

1294 a1, b1 = a2, b2

1295

1296 if _MODS.booleans.isBoolean(points):

1297 if not closed:

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

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

1300 else:

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

1302 return _radians2m(r, radius)

1303

1304

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

1306 excess=excessAbc_,

1307 wrap=False):

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

1309

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

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

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

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

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

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

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

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

1318 or C{None}.

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

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

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

1322 longitudes (C{bool}).

1323

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

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

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

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

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

1329 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with

1330 I{spherical deficit} C{D} and I{spherical excess}

1331 C{E} as the C{unit area}, all in C{radians}.

1332 '''

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

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

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

1336 excess=excess, wrap=wrap)

1337 return _t7Tuple(t, radius)

1338

1339

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

1341 wrap=False):

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

1343 excess} of a (spherical) triangle.

1344

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

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

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

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

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

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

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

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

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

1354 longitudinal deltas (C{bool}).

1355

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

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

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

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

1360 '''

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

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

1363 a = vincentys_(phiC, phiB, d)

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

1365

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

1367 s = sb * sc

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

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

1370

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

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

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

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

1375

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

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

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

1379

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

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

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

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

1384

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

1386

1387

1388 def _t7Tuple(t, radius):

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

1390 '''

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

1392 r = radius if _isRadius(radius) else \

1393 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1396 B, (r * t.b),

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

1398 return t

1399

1400

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

1402 areaOf, # functions

1403 intersecant2, intersection, intersections2, ispolar,

1404 isPoleEnclosedBy, # DEPRECATED, use ispolar

1405 meanOf,

1406 nearestOn2, nearestOn3,

1407 perimeterOf,

1408 sumOf, # XXX == vector3d.sumOf

1409 triangle7, triangle8_)

1410

1411# **) MIT License

1412#

1413# Copyright (C) 2016-2025 -- mrJean1 at Gmail -- All Rights Reserved.

1414#

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

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

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

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

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

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

1421#

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

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

1424#

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

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

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

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

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

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

1431# OTHER DEALINGS IN THE SOFTWARE.