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, isscalar, 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, fsum_

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_, _LatLon_, _near_, _not_, _null_, \

35 _points_, _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.namedTuples import LatLon2Tuple, LatLon3Tuple, \

39 NearestOn3Tuple, Triangle7Tuple, \

40 Triangle8Tuple

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

42 notImplemented, Fmt as _Fmt # XXX shadowed

43 from pygeodesy.props import deprecated_function, deprecated_method

44 from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, _intersecant2, \

45 LatLonSphericalBase, _rads3, _trilaterate5

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

47 from pygeodesy.units import Bearing_, Height, Lam_, Phi_, Radius, \

48 Radius_, Scalar

49 from pygeodesy.utily import acos1, asin1, degrees90, degrees180, degrees2m, \

50 m2radians, radiansPI2, sincos2_, tan_2, unrollPI, \

51 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__ = '23.03.30'

58

59_parallel_ = 'parallel'

60_path_ = 'path'

61

62_PI_EPS4 = PI - EPS4

63 if _PI_EPS4 >= PI:

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

65

66

67 class Cartesian(CartesianSphericalBase):

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

69 spherical, geodetic L{LatLon}.

70 '''

71

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

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

74

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

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

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

78

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

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

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

82

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

84 '''

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

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

87

88

89 class LatLon(LatLonSphericalBase):

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

91

92 @example:

93

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

95 '''

96

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

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

99 the closest point on the great circle path defined by a

100 start and an end point.

101

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

103 great circle path, the along-track distance is the distance

104 from the start point to the point where the perpendicular

105 crosses the path.

106

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

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

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

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

111

112 @return: Distance along the great circle path (C{meter},

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

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

115 toward the B{C{end}} point of the path or negative

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

117

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

119

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

121

122 @example:

123

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

125

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

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

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

129 '''

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

131 cx = cos(x)

132 return _0_0 if isnear0(cx) else \

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

134

135 @deprecated_method

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

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

138 '''

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

140

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

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

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

144

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

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

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

148

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

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

151 '''

152 self.others(other)

153

154 a1, b1 = self.philam

155 a2, b2 = other.philam

156

157 a = radians(lat)

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

159

160 sa, ca, sa1, ca1, \

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

162 sa1 *= ca2 * ca

163

164 x = sa1 * sdb

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

166 z = ca1 * sdb * ca2 * sa

167

168 h = hypot(x, y)

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

170 return None # great circle doesn't reach latitude

171

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

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

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

175

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

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

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

179

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

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

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

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

184

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

186 or positive if to the right of the path) (C{meter},

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

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

189

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

191

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

193

194 @example:

195

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

197

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

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

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

201 '''

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

203 return _r2m(x, radius)

204

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

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

207 travelled the given distance on the given initial bearing.

208

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

210 B{C{radius}}).

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

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

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

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

215

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

217

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

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

220

221 @example:

222

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

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

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

226 '''

227 a, b = self.philam

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

229

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

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

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

233

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

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

236

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

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

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

240

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

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

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

244

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

246

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

248

249 @example:

250

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

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

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

254 '''

255 self.others(other)

256

257 a1, b1 = self.philam

258 a2, b2 = other.philam

259

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

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

262

263# @Property_RO

264# def Ecef(self):

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

266# '''

267# return _MODS.ecef.EcefKarney

268

269 def greatCircle(self, bearing):

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

271 on the given initial bearing from this point.

272

273 Direction of vector is such that initial bearing vector

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

275

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

277

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

279

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

281

282 @example:

283

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

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

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

287 '''

288 a, b = self.philam

289

290 t = Bearing_(bearing)

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

292

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

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

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

296

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

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

299 to an other point.

300

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

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

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

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

305 C{sphericalNvector.LatLon.initialBearingTo}.

306

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

308

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

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

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

312

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

314

315 @example:

316

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

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

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

320 '''

321 self.others(other)

322

323 a1, b1 = self.philam

324 a2, b2 = other.philam

325

326 # XXX behavior like sphericalNvector.LatLon.initialBearingTo

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

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

329

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

331

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

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

334 and an other point.

335

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

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

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

339 @kwarg height: Optional height, overriding the intermediate

340 height (C{meter}).

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

342

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

344

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

346

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

348

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

350

351 @example:

352

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

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

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

356 '''

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

358 if isnear0(f): # PYCHOK no cover

359 r = self

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

361 r = self.others(other)

362 else:

363 a1, b1 = self.philam

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

365

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

367 r = vincentys_(a2, a1, db)

368 sr = sin(r)

369 if isnon0(sr):

370 sa1, ca1, sa2, ca2, \

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

372

373 t = f * r

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

375 b = sin( t) # / sr superflous

376

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

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

379 z = a * sa1 + b * sa2

380

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

382 b = atan2(y, x)

383

384 else: # PYCHOK no cover

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

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

387

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

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

390 return r

391

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

393 '''Compute the intersection point of two paths, each defined

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

395

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

397 initial bearing at this point (compass

398 C{degrees360}).

399 @arg other: Start point of the other path (L{LatLon}).

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

401 the initial bearing at the other start point

402 (compass C{degrees360}).

403 @kwarg height: Optional height for intersection point,

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

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

406

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

408 intersection point might be the L{antipode} to

409 the returned result.

410

411 @raise IntersectionError: Ambiguous or infinite intersection

412 or colinear, parallel or otherwise

413 non-intersecting paths.

414

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

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

417

418 @raise ValueError: Invalid B{C{height}} or C{null} path.

419

420 @example:

421

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

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

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

425 '''

426 try:

427 s2 = self.others(other)

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

429 LatLon=self.classof)

430 except (TypeError, ValueError) as x:

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

432

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

434 height=None, wrap=True):

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

436 by a center point and radius.

437

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

439 see B{C{radius}}).

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

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

442 see B{C{radius}}).

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

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

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

446 B{C{radius}}).

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

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

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

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

451

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

453 instance. For abutting circles, both intersection

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

455

456 @raise IntersectionError: Concentric, antipodal, invalid or

457 non-intersecting circles.

458

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

460

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

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

463 '''

464 try:

465 c2 = self.others(other)

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

467 height=height, wrap=wrap,

468 LatLon=self.classof)

469 except (TypeError, ValueError) as x:

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

471

472 def isenclosedBy(self, points):

473 '''Check whether a (convex) polygon encloses this point.

474

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

476

477 @return: C{True} if the polygon encloses this point,

478 C{False} otherwise.

479

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

481

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

483

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

485

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

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

488 enclose a pole or wrap around the earth longitudinally.

489

490 @example:

491

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

493 >>> p = LatLon(45,1, 1.1)

494 >>> inside = p.isEnclosedBy(b) # True

495 '''

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

497 n0 = self._N_vector

498

499 v2 = Ps[0]._N_vector

500 v1 = Ps[1]._N_vector

501 gc1 = v2.cross(v1)

502 # check whether this point on same side of all

503 # polygon edges (to the left or right depending

504 # on anti-/clockwise polygon direction)

505 t0 = gc1.angleTo(n0) > PI_2

506 s0 = None

507 # get great-circle vector for each edge

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

509 v2 = p._N_vector

510 gc = v1.cross(v2)

511 v1 = v2

512

513 t = gc.angleTo(n0) > PI_2

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

515 return False # outside

516

517 # check for convex polygon: angle between

518 # gc vectors, signed by direction of n0

519 # (otherwise the test above is not reliable)

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

521 if s != s0:

522 if s0 is None:

523 s0 = s

524 else:

525 t = _Fmt.SQUARE(points=i)

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

527 gc1 = gc

528

529 return True # inside

530

531 @deprecated_method

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

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

534 return self.isenclosedBy(points)

535

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

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

538

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

540 @kwarg height: Optional height for midpoint, overriding

541 the mean height (C{meter}).

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

543 may be negative or greater than 1.0.

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

545

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

547

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

549

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

551

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

553

554 @example:

555

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

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

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

559 '''

560 if fraction is _0_5:

561 self.others(other)

562

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

564 a1, b1 = self.philam

565 a2, b2 = other.philam

566

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

568

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

570

571 x = ca2 * cdb + ca1

572 y = ca2 * sdb

573

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

575 b = atan2(y, x) + b1

576

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

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

579 else:

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

581 return r

582

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

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

585

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

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

588

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

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

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

592 @kwarg options: Optional keyword arguments for function

593 L{pygeodesy.equirectangular_}.

594

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

596

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

598 see function L{pygeodesy.equirectangular_}.

599

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

601 is not (yet) supported.

602

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

604

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

606

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

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

609 '''

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

611 within = _xkwds_pop(options, within=True)

612 if not within:

613 notImplemented(self, within=within)

614

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

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

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

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

619# return point1

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

621# if isnear0(d):

622# return point1 # or point2

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

624# return point2

625# f = a / d

626# if within:

627# if f > EPS1:

628# return point2

629# elif f < EPS:

630# return point1

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

632

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

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

635 **options)[0]

636

637 @deprecated_method

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

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

640

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

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

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

644 units of B{C{radius}}.

645 '''

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

647 return r.closest, r.distance

648

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

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

651

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

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

654

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

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

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

658 @kwarg options: Optional keyword arguments for function

659 L{pygeodesy.equirectangular_}.

660

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

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

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

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

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

666 function L{pygeodesy.compassAngle}.

667

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

669 see function L{pygeodesy.equirectangular_}.

670

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

672

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

674

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

676

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

678 and L{pygeodesy.nearestOn5}.

679 '''

680 if _LatLon_ in options:

681 raise _ValueError(**options)

682

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

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

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

686

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

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

689 coordinates.

690

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

692 and other keyword arguments, ignored

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

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

695 this L{Cartesian} class or specify

696 C{B{Cartesian}=None}.

697

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

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

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

701

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

703 '''

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

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

706

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

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

709 '''

710 self.others(start=start)

711 self.others(end=end)

712

713 r = Radius_(radius)

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

715

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

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

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

719 return r, x, (e - b)

720

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

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

723

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

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

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

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

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

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

730

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

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

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

734

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

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

737 '''

738 self.others(otherB=otherB)

739 self.others(otherC=otherC)

740

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

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

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

744

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

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

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

748 of three corresponding circles.

749

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

751 as B{C{radius}}).

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

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

754 B{C{radius}}).

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

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

757 B{C{radius}}).

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

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

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

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

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

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

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

765

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

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

768 the corresponding trilaterated points C{minPoint} and

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

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

771

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

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

774

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

776 respectively largest I{radial} overlap found.

777

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

779 nearest respectively farthest intersection margin.

780

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

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

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

784 largest B{C{distance#}}.

785

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

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

788 no intersection or all (near-)concentric for

789 C{B{area}=False}.

790

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

792

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

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

795 '''

796 return _trilaterate5(self, distance1,

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

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

799 area=area, radius=radius, eps=eps, wrap=wrap)

800

801

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

803

804

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

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

807 (with the pointsjoined by great circle arcs).

808

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

810 or L{BooleanGH}).

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

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

813 or C{None}.

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

815

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

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

818

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

820

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

822

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

824

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

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

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

828

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

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

831 L{ellipsoidalKarney.areaOf}.

832

833 @example:

834

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

836 >>> areaOf(b) # 8666058750.718977

837

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

839 >>> areaOf(c) # 6.18e9

840 '''

841 if _MODS.booleans.isBoolean(points):

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

843

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

845 p1 = p2 = Ps[0]

846 a1, b1 = p1.philam

847 ta1, z1 = tan_2(a1), None

848

849 A = Fsum() # mean phi

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

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

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

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

854 D = Fsum()

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

856 a2, b2 = p2.philam

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

858 ta2 = tan_2(a2)

859 A += a2

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

861 _1_0 + ta1 * ta2)

862 ta1, b1 = ta2, b2

863

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

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

866 if z1 is not None:

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

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

869 p1, z1 = p2, z2

870

871 R = fabs(E * _2_0)

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

873 R = fabs(R - PI2)

874 if radius:

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

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

877 return float(R)

878

879

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

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

882

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

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

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

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

887

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

889 '''

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

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

892 ca *= sr

893

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

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

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

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

898

899

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

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

902 # and similar logic in .ellipsoidalBaseDI._intersect3

903 a1, b1 = start.philam

904

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

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

907 else: # must be a point

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

909 hs.append(end.height)

910 a2, b2 = end.philam

911

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

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

914 raise _ValueError(_SPACE_(_path_ + _i_, _null_))

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

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

917 -(b1 + b2) * _0_5)

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

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

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

921 cb12 * cb21 + sb12 * sb21,

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

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

924

925

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

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

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

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

930

931

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

933 height=None, wrap=True):

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

935

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

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

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

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

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

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

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

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

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

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

946 circle} methods.

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

948 conventionally) or C{None}.

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

950

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

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

953 point C{is} this very instance.

954

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

956

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

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

959

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

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

962 '''

963 c = _T00.others(center=center)

964 p = _T00.others(point=point)

965 try:

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

967 height=height, wrap=wrap)

968 except (TypeError, ValueError) as x:

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

970

971

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

973 LatLon=None, **LatLon_kwds):

974 # (INTERNAL) Intersect two (spherical) path, see L{intersection}

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

976

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

978

979 a1, b1 = start1.philam

980 a2, b2 = start2.philam

981

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

983 r12 = vincentys_(a2, a1, db)

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

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

986

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

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

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

990

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

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

993 raise IntersectionError(_parallel_)

994 # handle domain error for equivalent longitudes,

995 # see also functions asin_safe and acos_safe at

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

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

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

999 if sin(db) > 0:

1000 t12, t21 = t1, PI2 - t2

1001 else:

1002 t12, t21 = PI2 - t1, t2

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

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

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

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

1007 raise IntersectionError(_infinite_)

1008 sx3 = sx1 * sx2

1009# XXX if sx3 < 0:

1010# XXX raise ValueError(_ambiguous_)

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

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

1013

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

1015 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

1019

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

1021 _N_vector_ = _MODS.nvectorBase._N_vector_

1022

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

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

1025 x = x1.cross(x2)

1026 if x.length < EPS: # [nearly] colinear or parallel paths

1027 raise IntersectionError(_colinear_)

1028 a, b = x.philam

1029 # choose intersection similar to sphericalNvector

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

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

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

1033

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

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

1036 intersection, LatLon, LatLon_kwds)

1037

1038

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

1040 LatLon=LatLon, **LatLon_kwds):

1041 '''Compute the intersection point of two paths, each defined

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

1043

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

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

1046 the initial bearing at the first start point

1047 (compass C{degrees360}).

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

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

1050 the initial bearing at the second start point

1051 (compass C{degrees360}).

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

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

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

1055 @kwarg LatLon: Optional class to return the intersection

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

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

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

1059

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

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

1062 height)}. An alternate intersection point might

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

1064

1065 @raise IntersectionError: Ambiguous or infinite intersection

1066 or colinear, parallel or otherwise

1067 non-intersecting paths.

1068

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

1070

1071 @raise ValueError: Invalid B{C{height}} or C{null} path.

1072

1073 @example:

1074

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

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

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

1078 '''

1079 s1 = _T00.others(start1=start1)

1080 s2 = _T00.others(start2=start2)

1081 try:

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

1083 LatLon=LatLon, **LatLon_kwds)

1084 except (TypeError, ValueError) as x:

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

1086

1087

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

1089 height=None, wrap=True,

1090 LatLon=LatLon, **LatLon_kwds):

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

1092 by a center point and a radius.

1093

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

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

1096 see B{C{radius}}).

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

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

1099 see B{C{radius}}).

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

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

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

1103 B{C{radius}}).

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

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

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

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

1108 @kwarg LatLon: Optional class to return the intersection

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

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

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

1112

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

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

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

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

1117

1118 @raise IntersectionError: Concentric, antipodal, invalid or

1119 non-intersecting circles.

1120

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

1122

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

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

1125

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

1127

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

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

1130 '''

1131 c1 = _T00.others(center1=center1)

1132 c2 = _T00.others(center2=center2)

1133 try:

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

1135 height=height, wrap=wrap,

1136 LatLon=LatLon, **LatLon_kwds)

1137 except (TypeError, ValueError) as x:

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

1139

1140

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

1142 height=None, too_d=None, wrap=True,

1143 LatLon=LatLon, **LatLon_kwds):

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

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

1146

1147 def _dest3(bearing, h):

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

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

1150 intersections2, LatLon, LatLon_kwds)

1151

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

1153 if f: # swapped

1154 c1, c2 = c2, c1 # PYCHOK swap

1155

1156 a1, b1 = c1.philam

1157 a2, b2 = c2.philam

1158

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

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

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

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

1163

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

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

1166 if r < _0_0:

1167 raise _ValueError(eps=r)

1168

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

1170 if x > max(r, EPS):

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

1172 x = sd * sr1

1173 if isnear0(x):

1174 raise _ValueError(_invalid_)

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

1176

1177 elif x < r: # PYCHOK no cover

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

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

1180

1181 if height is None: # "radical height"

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

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

1184 else:

1185 h = Height(height)

1186

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

1188 if x < EPS4: # externally ...

1189 r = _dest3(b, h)

1190 elif x > _PI_EPS4: # internally ...

1191 r = _dest3(b + PI, h)

1192 else:

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

1194 return r, r # ... abutting circles

1195

1196

1197@deprecated_function

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

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

1200 '''

1201 return ispolar(points, wrap=wrap)

1202

1203

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

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

1206 '''

1207 n = func.__name__

1208 if LatLon is None:

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

1210 else:

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

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

1213 return r

1214

1215

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

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

1218

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

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

1221 height (C{meter}).

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

1223 or C{None}.

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

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

1226

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

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

1229 is C{None}.

1230

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

1232

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

1234 '''

1235 # geographic mean

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

1237 lat, lon = m._N_vector.latlon

1238

1239 if height is None:

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

1241 else: # PYCHOK no cover

1242 h = Height(height)

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

1244

1245

1246@deprecated_function

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

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

1249

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

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

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

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

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

1255 '''

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

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

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

1259 return ll, d

1260

1261

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

1263 LatLon=LatLon, **options):

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

1265

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

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

1268

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

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

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

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

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

1274 or C{None}.

1275 @kwarg options: Optional keyword arguments for function

1276 L{pygeodesy.equirectangular_}.

1277

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

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

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

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

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

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

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

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

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

1287

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

1289 see function L{pygeodesy.equirectangular_}.

1290

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

1292

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

1294

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

1296

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

1298 '''

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

1300 LatLon=None, **options)

1301 d = degrees2m(d, radius=radius)

1302 n = nearestOn3.__name__

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

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

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

1306

1307

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

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

1310 (with great circle arcs joining the points).

1311

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

1313 or L{BooleanGH}).

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

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

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

1317

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

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

1320

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

1322

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

1324

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

1326 C{B{points}} a composite.

1327

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

1329

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

1331 and L{ellipsoidalKarney.perimeterOf}.

1332 '''

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

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

1335 a1, b1 = Ps[0].philam

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

1337 a2, b2 = p.philam

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

1339 yield vincentys_(a2, a1, db)

1340 a1, b1 = a2, b2

1341

1342 if _MODS.booleans.isBoolean(points):

1343 if not closed:

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

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

1346 else:

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

1348 return _r2m(r, radius)

1349

1350

1351 def _r2m(r, radius):

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

1353 '''

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

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

1356 return r

1357

1358

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

1360 excess=excessAbc,

1361 wrap=False):

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

1363

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

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

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

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

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

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

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

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

1372 or C{None}.

1373 @kwarg excess: I{Spherical excess} callable (L{excessAbc},

1374 L{excessGirard} or L{excessLHuilier}).

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

1376

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

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

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

1380 in I{square} C{meter} or units of B{C{radius}}

1381 I{squared} or if B{C{radius}} is C{None}, a

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

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

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

1385 '''

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

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

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

1389 excess=excess, wrap=wrap)

1390 return _t7Tuple(t, radius)

1391

1392

1393 def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc,

1394 wrap=False):

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

1396 excess} of a (spherical) triangle.

1397

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

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

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

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

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

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

1404 @kwarg excess: I{Spherical excess} callable (L{excessAbc},

1405 L{excessGirard} or L{excessLHuilier}).

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

1407

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

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

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

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

1412 '''

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

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

1415 a = vincentys_(phiC, phiB, d)

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

1417

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

1419 s = sb * sc

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

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

1422

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

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

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

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

1427

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

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

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

1431

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

1433 E = excessGirard(A, B, C) if excess in (excessGirard, True) else (

1434 excessLHuilier(a, b, c) if excess in (excessLHuilier, False) else

1435 excessAbc(*max((A, b, c), (B, c, a), (C, a, b))))

1436

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

1438

1439

1440 def _t7Tuple(t, radius):

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

1442 '''

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

1444 r = radius if isscalar(radius) else \

1445 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1448 B, (r * t.b),

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

1450 return t

1451

1452

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

1454 areaOf, # functions

1455 intersecant2, intersection, intersections2, ispolar,

1456 isPoleEnclosedBy, # DEPRECATED, use ispolar

1457 meanOf,

1458 nearestOn2, nearestOn3,

1459 perimeterOf,

1460 sumOf, # == vector3d.sumOf

1461 triangle7, triangle8_)

1462

1463# **) MIT License

1464#

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

1466#

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

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

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

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

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

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

1473#

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

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

1476#

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

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

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

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

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

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

1483# OTHER DEALINGS IN THE SOFTWARE.