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 _TypeError, _ValueError, IntersectionError, \

26 _xError, _xkwds, _xkwds_get, _xkwds_pop2

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_, \

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, atan2d, degrees90, degrees180, \

50 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \

51 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI

52 from pygeodesy.vector3d import sumOf, Vector3d

53

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

55

56__all__ = _ALL_LAZY.sphericalTrigonometry

57__version__ = '24.08.13'

58

59_PI_EPS4 = PI - EPS4

60 if _PI_EPS4 >= PI:

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

62

63

64 class Cartesian(CartesianSphericalBase):

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

66 spherical, geodetic L{LatLon}.

67 '''

68

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

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

71

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

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

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

75

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

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

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

79

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

81 '''

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

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

84

85

86 class LatLon(LatLonSphericalBase):

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

88 '''

89

90 def _ab1_ab2_db5(self, other, wrap):

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

92 '''

93 a1, b1 = self.philam

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

95 if wrap:

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

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

98 else: # unrollPI shortcut

99 db = b2 - b1

100 return a1, b1, a2, b2, db

101

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

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

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

105 end point.

106

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

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

109 from the start point to the point where the perpendicular

110 crosses the line.

111

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

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

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

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

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

117

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

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

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

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

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

123 B{C{start}} point.

124

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

126

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

128 '''

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

130 cx = cos(x)

131 return _0_0 if isnear0(cx) else \

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

133

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

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

136 '''

137 s = self.others(start=start)

138 e = self.others(end=end)

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

140

141 r = Radius_(radius)

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

143

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

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

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

147 return r, x, -b

148

149 @deprecated_method

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

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

152 '''

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

154

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

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

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

158

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

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

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

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

163

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

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

166 '''

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

168 sa, ca, sa1, ca1, \

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

170 sa1 *= ca2 * ca

171

172 x = sa1 * sdb

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

174 z = ca1 * sdb * ca2 * sa

175

176 h = hypot(x, y)

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

178 return None # great circle doesn't reach latitude

179

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

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

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

183

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

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

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

187

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

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

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

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

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

193

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

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

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

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

198

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

200

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

202 '''

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

204 return _radians2m(x, radius)

205

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

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

208 travelled the given distance on the given initial bearing.

209

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

211 B{C{radius}}).

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

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

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

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

216

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

218

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

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

221 '''

222 a, b = self.philam

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

224

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

226 h = self._heigHt(height)

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

228

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

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

231

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

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

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

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

236

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

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

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

240

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

242

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

244 '''

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

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

247

248# @Property_RO

249# def Ecef(self):

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

251# '''

252# return _MODS.ecef.EcefKarney

253

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

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

256 on the given initial bearing from this point.

257

258 Direction of vector is such that initial bearing vector

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

260

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

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

263 overriding the default L{Vector3d}.

264 @kwarg Vector_kwds: Optional, additional keyword argunents

265 for B{C{Vector}}.

266

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

268

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

270 '''

271 a, b = self.philam

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

273

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

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

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

277

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

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

280 to an other point.

281

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

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

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

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

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

287 C{sphericalNvector.LatLon.initialBearingTo}.

288

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

290

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

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

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

294

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

296 '''

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

298 # XXX behavior like sphericalNvector.LatLon.initialBearingTo

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

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

301

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

303

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

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

306 and an other point.

307

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

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

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

311 @kwarg height: Optional height, overriding the intermediate

312 height (C{meter}).

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

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

315

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

317

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

319

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

321

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

323 '''

324 p = self

325 f = Scalar(fraction=fraction)

326 if not isnear0(f):

327 p = p.others(other)

328 if wrap:

329 p = _Wrap.point(p)

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

331 a1, b1 = self.philam

332 a2, b2 = p.philam

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

334 r = vincentys_(a2, a1, db)

335 sr = sin(r)

336 if isnon0(sr):

337 sa1, ca1, sa2, ca2, \

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

339

340 t = f * r

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

342 b = sin( t) # / sr superflous

343

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

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

346 z = a * sa1 + b * sa2

347

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

349 b = atan2d(y, x)

350

351 else: # PYCHOK no cover

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

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

354

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

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

357 return p

358

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

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

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

362

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

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

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

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

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

368 C{degrees360}).

369 @kwarg height: Optional height for intersection point,

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

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

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

373

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

375 intersection point might be the L{antipode} to

376 the returned result.

377

378 @raise IntersectionError: Ambiguous or infinite intersection

379 or colinear, parallel or otherwise

380 non-intersecting lines.

381

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

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

384

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

386 '''

387 try:

388 s2 = self.others(other)

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

390 LatLon=self.classof)

391 except (TypeError, ValueError) as x:

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

393 other=other, end2=end2, wrap=wrap)

394

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

396 height=None, wrap=True):

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

398 by a center point and radius.

399

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

401 see B{C{radius}}).

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

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

404 see B{C{radius}}).

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

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

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

408 B{C{radius}}).

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

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

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

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

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

414

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

416 instance. For abutting circles, both intersection

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

418

419 @raise IntersectionError: Concentric, antipodal, invalid or

420 non-intersecting circles.

421

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

423

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

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

426 '''

427 try:

428 c2 = self.others(other)

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

430 height=height, wrap=wrap,

431 LatLon=self.classof)

432 except (TypeError, ValueError) as x:

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

434 other=other, rad2=rad2, wrap=wrap)

435

436 @deprecated_method

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

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

439 return self.isenclosedBy(points)

440

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

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

443

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

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

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

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

448

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

450 composite, C{False} otherwise.

451

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

453

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

455

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

457

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

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

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

461 '''

462 if _MODS.booleans.isBoolean(points):

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

464

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

466 n0 = self._N_vector

467

468 v2 = Ps[0]._N_vector

469 p1 = Ps[1]

470 v1 = p1._N_vector

471 # check whether this point on same side of all

472 # polygon edges (to the left or right depending

473 # on the anti-/clockwise polygon direction)

474 gc1 = v2.cross(v1)

475 t0 = gc1.angleTo(n0) > PI_2

476 s0 = None

477 # get great-circle vector for each edge

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

479 if wrap and not Ps.looped:

480 p2 = _unrollon(p1, p2)

481 p1 = p2

482 v2 = p2._N_vector

483 gc = v1.cross(v2)

484 t = gc.angleTo(n0) > PI_2

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

486 return False # outside

487

488 # check for convex polygon: angle between

489 # gc vectors, signed by direction of n0

490 # (otherwise the test above is not reliable)

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

492 if s != s0:

493 if s0 is None:

494 s0 = s

495 else:

496 t = _Fmt.SQUARE(points=i)

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

498 gc1, v1 = gc, v2

499

500 return True # inside

501

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

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

504

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

506 @kwarg height: Optional height for midpoint, overriding

507 the mean height (C{meter}).

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

509 may be negative or greater than 1.0.

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

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

512

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

514

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

516

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

518

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

520 '''

521 if fraction is _0_5:

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

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

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

525

526 x = ca2 * cdb + ca1

527 y = ca2 * sdb

528

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

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

531

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

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

534 else:

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

536 return r

537

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

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

540

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

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

543

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

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

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

547 @kwarg wrap_adjust_limit: Optional keyword arguments for functions

548 L{sphericalTrigonometry.nearestOn3} and

549 L{pygeodesy.equirectangular4},

550

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

552

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

554 see function L{pygeodesy.equirectangular4}.

555

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

557 is not (yet) supported.

558

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

560

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

562

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

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

565 '''

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

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

568 if not w:

569 self._notImplemented(within=w)

570

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

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

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

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

575# return point1

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

577# if isnear0(d):

578# return point1 # or point2

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

580# return point2

581# f = a / d

582# if within:

583# if f > EPS1:

584# return point2

585# elif f < EPS:

586# return point1

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

588

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

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

591 **kwds)[0]

592

593 @deprecated_method

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

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

596

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

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

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

600 units of B{C{radius}}.

601 '''

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

603 return r.closest, r.distance

604

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

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

607

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

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

610

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

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

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

614 @kwarg wrap_adjust_limit: Optional keyword arguments for function

615 L{sphericalTrigonometry.nearestOn3} and

616 L{pygeodesy.equirectangular4},

617

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

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

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

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

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

623 function L{pygeodesy.compassAngle}.

624

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

626 see function L{pygeodesy.equirectangular4}.

627

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

629

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

631

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

633

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

635 and L{pygeodesy.nearestOn5}.

636 '''

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

638 LatLon=self.classof, **wrap_adjust_limit)

639

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

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

642

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

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

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

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

647

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

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

650 and C{M} if available.

651

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

653 '''

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

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

656

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

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

659

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

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

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

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

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

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

666

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

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

669

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

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

672 '''

673 B = self.others(otherB=otherB)

674 C = self.others(otherC=otherC)

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

676

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

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

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

680

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

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

683 at this and at the other point.

684

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

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

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

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

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

690

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

692

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

694 '''

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

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

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

698

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

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

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

702 of three corresponding circles.

703

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

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

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

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

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

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

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

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

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

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

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

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

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

717

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

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

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

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

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

723

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

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

726

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

728 largest I{radial} overlap found.

729

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

731 respectively farthest intersection margin.

732

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

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

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

736

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

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

739 no intersection or all (near-)concentric if

740 C{B{area}=False}.

741

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

743

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

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

746 or B{C{radius}}.

747 '''

748 return _trilaterate5(self, distance1,

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

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

751 area=area, radius=radius, eps=eps, wrap=wrap)

752

753

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

755

756

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

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

759 (with the pointsjoined by great circle arcs).

760

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

762 or L{BooleanGH}).

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

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

765 or C{None}.

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

767 (C{bool}).

768

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

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

771

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

773

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

775

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

777

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

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

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

781

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

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

784 L{ellipsoidalKarney.areaOf}.

785 '''

786 if _MODS.booleans.isBoolean(points):

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

788

789 _at2, _t_2 = atan2, tan_2

790 _un, _w180 = unrollPI, wrap180

791

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

793 p1 = p2 = Ps[0]

794 a1, b1 = p1.philam

795 ta1, z1 = _t_2(a1), None

796

797 A = Fsum() # mean phi

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

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

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

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

802 D = Fsum()

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

804 a2, b2 = p2.philam

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

806 A += a2

807 ta2 = _t_2(a2)

808 tdb = _t_2(db, points=i)

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

810 _1_0 + ta1 * ta2)

811 ta1, b1 = ta2, b2

812

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

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

815 if z1 is not None:

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

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

818 p1, z1 = p2, z2

819

820 R = abs(R * _2_0)

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

822 R = abs(R - PI2)

823 if radius:

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

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

826 return float(R)

827

828

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

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

831

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

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

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

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

836

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

838 '''

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

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

841 ca *= sr

842

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

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

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

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

847

848

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

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

851 # and similar logic in .ellipsoidalBaseDI._intersect3

852 a1, b1 = s.philam

853

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

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

856 else: # must be a point

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

858 hs.append(end.height)

859 a2, b2 = end.philam

860 if wrap:

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

862

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

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

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

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

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

868 -(b1 + b2) * _0_5)

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

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

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

872 cb12 * cb21 + sb12 * sb21,

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

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

875

876

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

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

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

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

881

882

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

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

885 two points or as a point and bearing.

886

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

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

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

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

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

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

893 @kwarg radius_exact_height_wrap: Optional keyword arguments, see

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

895

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

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

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

899

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

901

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

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

904

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

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

907 '''

908 c = _T00.others(center=center)

909 p = _T00.others(point=point)

910 try:

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

912 except (TypeError, ValueError) as x:

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

914 **radius_exact_height_wrap)

915

916

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

918 LatLon=None, **LatLon_kwds):

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

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

921

922 s1, s2 = start1, start2

923 if wrap:

924 s2 = _Wrap.point(s2)

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

926

927 a1, b1 = s1.philam

928 a2, b2 = s2.philam

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

930 r12 = vincentys_(a2, a1, db)

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

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

933

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

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

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

937

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

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

940 raise IntersectionError(_parallel_)

941 # handle domain error for equivalent longitudes,

942 # see also functions asin_safe and acos_safe at

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

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

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

946 if sin(db) > 0:

947 t21 = PI2 - t21

948 else:

949 t12 = PI2 - t12

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

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

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

953 raise IntersectionError(_infinite_)

954 sx3 = sx1 * sx2

955# XXX if sx3 < 0:

956# XXX raise ValueError(_ambiguous_)

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

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

959

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

961 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

965

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

967 _N_vector_ = _MODS.nvectorBase._N_vector_

968

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

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

971 x = x1.cross(x2)

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

973 raise IntersectionError(_colinear_)

974 a, b = x.philam

975 # choose intersection similar to sphericalNvector

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

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

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

979

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

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

982 intersection, LatLon, LatLon_kwds)

983

984

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

986 LatLon=LatLon, **LatLon_kwds):

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

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

989

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

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

992 the initial bearing at the first start point

993 (compass C{degrees360}).

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

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

996 the initial bearing at the second start point

997 (compass C{degrees360}).

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

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

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

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

1002 @kwarg LatLon: Optional class to return the intersection

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

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

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

1006

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

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

1009 height)}. An alternate intersection point might

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

1011

1012 @raise IntersectionError: Ambiguous or infinite intersection

1013 or colinear, parallel or otherwise

1014 non-intersecting lines.

1015

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

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

1018

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

1020 '''

1021 s1 = _T00.others(start1=start1)

1022 s2 = _T00.others(start2=start2)

1023 try:

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

1025 LatLon=LatLon, **LatLon_kwds)

1026 except (TypeError, ValueError) as x:

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

1028

1029

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

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

1032 LatLon=LatLon, **LatLon_kwds):

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

1034 by a center point and a radius.

1035

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

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

1038 see B{C{radius}}).

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

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

1041 see B{C{radius}}).

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

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

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

1045 B{C{radius}}).

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

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

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

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

1050 (C{bool}).

1051 @kwarg LatLon: Optional class to return the intersection

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

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

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

1055

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

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

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

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

1060

1061 @raise IntersectionError: Concentric, antipodal, invalid or

1062 non-intersecting circles.

1063

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

1065

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

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

1068

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

1070

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

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

1073 '''

1074 c1 = _T00.others(center1=center1)

1075 c2 = _T00.others(center2=center2)

1076 try:

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

1078 height=height, wrap=wrap,

1079 LatLon=LatLon, **LatLon_kwds)

1080 except (TypeError, ValueError) as x:

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

1082 center2=center2, rad2=rad2, wrap=wrap)

1083

1084

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

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

1087 LatLon=LatLon, **LatLon_kwds):

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

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

1090

1091 def _dest3(bearing, h):

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

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

1094 intersections2, LatLon, LatLon_kwds)

1095

1096 a1, b1 = c1.philam

1097 a2, b2 = c2.philam

1098 if wrap:

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

1100

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

1102 if f: # swapped radii, swap centers

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

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

1105

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

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

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

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

1110

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

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

1113 if r < _0_0:

1114 raise _ValueError(eps=r)

1115

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

1117 if x > max(r, EPS):

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

1119 x = sd * sr1

1120 if isnear0(x):

1121 raise _ValueError(_invalid_)

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

1123

1124 elif x < r: # PYCHOK no cover

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

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

1127

1128 if height is None: # "radical height"

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

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

1131 else:

1132 h = Height(height)

1133

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

1135 if x < EPS4: # externally ...

1136 r = _dest3(b, h)

1137 elif x > _PI_EPS4: # internally ...

1138 r = _dest3(b + PI, h)

1139 else:

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

1141 return r, r # ... abutting circles

1142

1143

1144@deprecated_function

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

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

1147 '''

1148 return ispolar(points, wrap=wrap)

1149

1150

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

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

1153 '''

1154 n = where.__name__

1155 if LatLon is None:

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

1157 else:

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

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

1160 return r

1161

1162

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

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

1165

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

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

1168 (C{meter}).

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

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

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

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

1173

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

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

1176

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

1178

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

1180 '''

1181 def _N_vs(ps, w):

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

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

1184 yield p._N_vector

1185

1186 m = _MODS.nvectorBase

1187 # geographic, vectorial mean

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

1189 lat, lon, h = n.latlonheight

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

1191

1192

1193@deprecated_function

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

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

1196

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

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

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

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

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

1202 '''

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

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

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

1206 return ll, d

1207

1208

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

1210 limit=9, **LatLon_and_kwds):

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

1212

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

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

1215

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

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

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

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

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

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

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

1223 (C{bool}).

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

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

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

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

1228 closest point and optional, additional C{B{LatLon}} keyword

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

1230

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

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

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

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

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

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

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

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

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

1240

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

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

1243

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

1245

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

1247

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

1249

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

1251 L{nearestOn5<pygeodesy.nearestOn5>}.

1252 '''

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

1254 adjust=adjust, limit=limit)

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

1256 h = t.height

1257 n = nearestOn3.__name__

1258

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

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

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

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

1263

1264

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

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

1267 (with great circle arcs joining the points).

1268

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

1270 or L{BooleanGH}).

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

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

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

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

1275

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

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

1278

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

1280

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

1282

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

1284 C{B{points}} a composite.

1285

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

1287

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

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

1290 '''

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

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

1293 a1, b1 = Ps[0].philam

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

1295 a2, b2 = p.philam

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

1297 yield vincentys_(a2, a1, db)

1298 a1, b1 = a2, b2

1299

1300 if _MODS.booleans.isBoolean(points):

1301 if not closed:

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

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

1304 else:

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

1306 return _radians2m(r, radius)

1307

1308

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

1310 excess=excessAbc_,

1311 wrap=False):

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

1313

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

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

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

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

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

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

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

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

1322 or C{None}.

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

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

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

1326 longitudes (C{bool}).

1327

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

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

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

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

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

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

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

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

1336 '''

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

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

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

1340 excess=excess, wrap=wrap)

1341 return _t7Tuple(t, radius)

1342

1343

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

1345 wrap=False):

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

1347 excess} of a (spherical) triangle.

1348

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

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

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

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

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

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

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

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

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

1358 longitudinal deltas (C{bool}).

1359

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

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

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

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

1364 '''

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

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

1367 a = vincentys_(phiC, phiB, d)

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

1369

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

1371 s = sb * sc

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

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

1374

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

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

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

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

1379

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

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

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

1383

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

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

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

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

1388

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

1390

1391

1392 def _t7Tuple(t, radius):

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

1394 '''

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

1396 r = radius if _isRadius(radius) else \

1397 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1400 B, (r * t.b),

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

1402 return t

1403

1404

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

1406 areaOf, # functions

1407 intersecant2, intersection, intersections2, ispolar,

1408 isPoleEnclosedBy, # DEPRECATED, use ispolar

1409 meanOf,

1410 nearestOn2, nearestOn3,

1411 perimeterOf,

1412 sumOf, # XXX == vector3d.sumOf

1413 triangle7, triangle8_)

1414

1415# **) MIT License

1416#

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

1418#

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

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

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

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

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

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

1425#

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

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

1428#

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

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

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

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

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

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

1435# OTHER DEALINGS IN THE SOFTWARE.