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-2016} 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 _ValueError, IntersectionError, _xError, \

26 _xkwds, _xkwds_get, _xkwds_pop

27 from pygeodesy.fmath import favg, 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_, _invalid_,\

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

35 _SPACE_, _too_

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

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

40 Triangle7Tuple, Triangle8Tuple

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

42 Fmt as _Fmt, notImplemented # XXX shadowed

43 from pygeodesy.props import deprecated_function, deprecated_method

44 from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \

45 _intersecant2, LatLonSphericalBase, \

46 _rads3, _radians2m, _trilaterate5

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

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

49 Phi_, Radius_, Scalar

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

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

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

53 from pygeodesy.vector3d import sumOf, Vector3d

54

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

56

57__all__ = _ALL_LAZY.sphericalTrigonometry

58__version__ = '23.12.03'

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

87 class LatLon(LatLonSphericalBase):

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

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 trilaterate5(self, distance1, point2, distance2, point3, distance3,

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

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

690 intersection} of three corresponding circles.

691

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

693 as B{C{radius}}).

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

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

696 B{C{radius}}).

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

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

699 B{C{radius}}).

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

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

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

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

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

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

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

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

708

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

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

711 the corresponding trilaterated points C{minPoint} and

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

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

714

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

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

717

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

719 respectively largest I{radial} overlap found.

720

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

722 nearest respectively farthest intersection margin.

723

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

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

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

727 largest B{C{distance#}}.

728

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

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

731 no intersection or all (near-)concentric for

732 C{B{area}=False}.

733

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

735

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

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

738 or B{C{radius}}.

739 '''

740 return _trilaterate5(self, distance1,

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

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

743 area=area, radius=radius, eps=eps, wrap=wrap)

744

745

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

747

748

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

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

751 (with the pointsjoined by great circle arcs).

752

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

754 or L{BooleanGH}).

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

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

757 or C{None}.

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

759 (C{bool}).

760

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

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

763

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

765

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

767

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

769

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

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

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

773

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

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

776 L{ellipsoidalKarney.areaOf}.

777 '''

778 if _MODS.booleans.isBoolean(points):

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

780

781 _at2, _t_2 = atan2, tan_2

782 _un, _w180 = unrollPI, wrap180

783

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

785 p1 = p2 = Ps[0]

786 a1, b1 = p1.philam

787 ta1, z1 = _t_2(a1), None

788

789 A = Fsum() # mean phi

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

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

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

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

794 D = Fsum()

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

796 a2, b2 = p2.philam

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

798 A += a2

799 ta2 = _t_2(a2)

800 tdb = _t_2(db, points=i)

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

802 _1_0 + ta1 * ta2)

803 ta1, b1 = ta2, b2

804

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

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

807 if z1 is not None:

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

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

810 p1, z1 = p2, z2

811

812 R = abs(R * _2_0)

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

814 R = abs(R - PI2)

815 if radius:

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

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

818 return float(R)

819

820

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

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

823

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

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

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

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

828

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

830 '''

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

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

833 ca *= sr

834

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

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

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

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

839

840

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

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

843 # and similar logic in .ellipsoidalBaseDI._intersect3

844 a1, b1 = s.philam

845

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

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

848 else: # must be a point

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

850 hs.append(end.height)

851 a2, b2 = end.philam

852 if wrap:

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

854

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

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

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

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

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

860 -(b1 + b2) * _0_5)

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

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

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

864 cb12 * cb21 + sb12 * sb21,

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

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

867

868

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

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

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

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

873

874

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

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

877 two points or as a point and bearing.

878

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

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

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

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

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

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

885 @kwarg radius_exact_height_wrap: Optional keyword arguments, see

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

887

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

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

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

891

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

893

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

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

896

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

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

899 '''

900 c = _T00.others(center=center)

901 p = _T00.others(point=point)

902 try:

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

904 except (TypeError, ValueError) as x:

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

906 **radius_exact_height_wrap)

907

908

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

910 LatLon=None, **LatLon_kwds):

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

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

913

914 s1, s2 = start1, start2

915 if wrap:

916 s2 = _Wrap.point(s2)

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

918

919 a1, b1 = s1.philam

920 a2, b2 = s2.philam

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

922 r12 = vincentys_(a2, a1, db)

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

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

925

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

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

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

929

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

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

932 raise IntersectionError(_parallel_)

933 # handle domain error for equivalent longitudes,

934 # see also functions asin_safe and acos_safe at

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

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

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

938 if sin(db) > 0:

939 t21 = PI2 - t21

940 else:

941 t12 = PI2 - t12

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

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

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

945 raise IntersectionError(_infinite_)

946 sx3 = sx1 * sx2

947# XXX if sx3 < 0:

948# XXX raise ValueError(_ambiguous_)

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

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

951

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

953 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

957

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

959 _N_vector_ = _MODS.nvectorBase._N_vector_

960

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

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

963 x = x1.cross(x2)

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

965 raise IntersectionError(_colinear_)

966 a, b = x.philam

967 # choose intersection similar to sphericalNvector

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

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

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

971

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

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

974 intersection, LatLon, LatLon_kwds)

975

976

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

978 LatLon=LatLon, **LatLon_kwds):

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

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

981

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

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

984 the initial bearing at the first start point

985 (compass C{degrees360}).

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

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

988 the initial bearing at the second start point

989 (compass C{degrees360}).

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

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

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

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

994 @kwarg LatLon: Optional class to return the intersection

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

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

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

998

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

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

1001 height)}. An alternate intersection point might

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

1003

1004 @raise IntersectionError: Ambiguous or infinite intersection

1005 or colinear, parallel or otherwise

1006 non-intersecting lines.

1007

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

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

1010

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

1012 '''

1013 s1 = _T00.others(start1=start1)

1014 s2 = _T00.others(start2=start2)

1015 try:

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

1017 LatLon=LatLon, **LatLon_kwds)

1018 except (TypeError, ValueError) as x:

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

1020

1021

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

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

1024 LatLon=LatLon, **LatLon_kwds):

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

1026 by a center point and a radius.

1027

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

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

1030 see B{C{radius}}).

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

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

1033 see B{C{radius}}).

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

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

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

1037 B{C{radius}}).

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

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

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

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

1042 (C{bool}).

1043 @kwarg LatLon: Optional class to return the intersection

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

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

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

1047

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

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

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

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

1052

1053 @raise IntersectionError: Concentric, antipodal, invalid or

1054 non-intersecting circles.

1055

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

1057

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

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

1060

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

1062

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

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

1065 '''

1066 c1 = _T00.others(center1=center1)

1067 c2 = _T00.others(center2=center2)

1068 try:

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

1070 height=height, wrap=wrap,

1071 LatLon=LatLon, **LatLon_kwds)

1072 except (TypeError, ValueError) as x:

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

1074 center2=center2, rad2=rad2, wrap=wrap)

1075

1076

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

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

1079 LatLon=LatLon, **LatLon_kwds):

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

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

1082

1083 def _dest3(bearing, h):

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

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

1086 intersections2, LatLon, LatLon_kwds)

1087

1088 a1, b1 = c1.philam

1089 a2, b2 = c2.philam

1090 if wrap:

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

1092

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

1094 if f: # swapped radii, swap centers

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

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

1097

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

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

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

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

1102

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

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

1105 if r < _0_0:

1106 raise _ValueError(eps=r)

1107

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

1109 if x > max(r, EPS):

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

1111 x = sd * sr1

1112 if isnear0(x):

1113 raise _ValueError(_invalid_)

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

1115

1116 elif x < r: # PYCHOK no cover

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

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

1119

1120 if height is None: # "radical height"

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

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

1123 else:

1124 h = Height(height)

1125

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

1127 if x < EPS4: # externally ...

1128 r = _dest3(b, h)

1129 elif x > _PI_EPS4: # internally ...

1130 r = _dest3(b + PI, h)

1131 else:

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

1133 return r, r # ... abutting circles

1134

1135

1136@deprecated_function

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

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

1139 '''

1140 return ispolar(points, wrap=wrap)

1141

1142

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

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

1145 '''

1146 n = where.__name__

1147 if LatLon is None:

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

1149 else:

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

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

1152 return r

1153

1154

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

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

1157

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

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

1160 height (C{meter}).

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

1162 (C{bool}).

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

1164 or C{None}.

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

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

1167

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

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

1170 is C{None}.

1171

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

1173

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

1175 '''

1176 def _N_vs(ps, w):

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

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

1179 yield p._N_vector

1180

1181 m = _MODS.nvectorBase

1182 # geographic, vectorial mean

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

1184 lat, lon, h = n.latlonheight

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

1186

1187

1188@deprecated_function

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

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

1191

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

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

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

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

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

1197 '''

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

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

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

1201 return ll, d

1202

1203

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

1205 limit=9, **LatLon_and_kwds):

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

1207

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

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

1210

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

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

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

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

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

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

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

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

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

1220 spherical earth radius L{R_KM}).

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

1222 or C{None}.

1223 @kwarg options: Optional keyword arguments for function

1224 L{pygeodesy.equirectangular_}.

1225

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

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

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

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

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

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

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

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

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

1235

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

1237 see function L{pygeodesy.equirectangular_}.

1238

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

1240

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

1242

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

1244

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

1246 '''

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

1248 adjust=adjust, limit=limit)

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

1250 h = t.height

1251 n = nearestOn3.__name__

1252

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

1254 LL = _xkwds_pop(kwds, LatLon=LatLon)

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

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

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

1258

1259

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

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

1262 (with great circle arcs joining the points).

1263

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

1265 or L{BooleanGH}).

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

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

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

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

1270

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

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

1273

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

1275

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

1277

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

1279 C{B{points}} a composite.

1280

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

1282

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

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

1285 '''

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

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

1288 a1, b1 = Ps[0].philam

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

1290 a2, b2 = p.philam

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

1292 yield vincentys_(a2, a1, db)

1293 a1, b1 = a2, b2

1294

1295 if _MODS.booleans.isBoolean(points):

1296 if not closed:

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

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

1299 else:

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

1301 return _radians2m(r, radius)

1302

1303

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

1305 excess=excessAbc_,

1306 wrap=False):

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

1308

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

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

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

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

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

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

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

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

1317 or C{None}.

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

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

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

1321 longitudes (C{bool}).

1322

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

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

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

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

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

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

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

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

1331 '''

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

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

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

1335 excess=excess, wrap=wrap)

1336 return _t7Tuple(t, radius)

1337

1338

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

1340 wrap=False):

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

1342 excess} of a (spherical) triangle.

1343

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

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

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

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

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

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

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

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

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

1353 longitudinal deltas (C{bool}).

1354

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

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

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

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

1359 '''

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

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

1362 a = vincentys_(phiC, phiB, d)

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

1364

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

1366 s = sb * sc

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

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

1369

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

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

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

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

1374

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

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

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

1378

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

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

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

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

1383

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

1385

1386

1387 def _t7Tuple(t, radius):

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

1389 '''

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

1391 r = radius if _isRadius(radius) else \

1392 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1395 B, (r * t.b),

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

1397 return t

1398

1399

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

1401 areaOf, # functions

1402 intersecant2, intersection, intersections2, ispolar,

1403 isPoleEnclosedBy, # DEPRECATED, use ispolar

1404 meanOf,

1405 nearestOn2, nearestOn3,

1406 perimeterOf,

1407 sumOf, # XXX == vector3d.sumOf

1408 triangle7, triangle8_)

1409

1410# **) MIT License

1411#

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

1413#

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

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

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

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

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

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

1420#

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

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

1423#

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

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

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

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

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

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

1430# OTHER DEALINGS IN THE SOFTWARE.