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.04.05'

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 or composite encloses this point.

474

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

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

477

478 @return: C{True} if this point is inside the polygon or composite,

479 C{False} otherwise.

480

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

482

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

484

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

486

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

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

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

490 '''

491 if _MODS.booleans.isBoolean(points):

492 return points._encloses(self.lat, self.lon)

493

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

495 n0 = self._N_vector

496

497 v2 = Ps[0]._N_vector

498 v1 = Ps[1]._N_vector

499 gc1 = v2.cross(v1)

500 # check whether this point on same side of all

501 # polygon edges (to the left or right depending

502 # on anti-/clockwise polygon direction)

503 t0 = gc1.angleTo(n0) > PI_2

504 s0 = None

505 # get great-circle vector for each edge

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

507 v2 = p._N_vector

508 gc = v1.cross(v2)

509 v1 = v2

510

511 t = gc.angleTo(n0) > PI_2

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

513 return False # outside

514

515 # check for convex polygon: angle between

516 # gc vectors, signed by direction of n0

517 # (otherwise the test above is not reliable)

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

519 if s != s0:

520 if s0 is None:

521 s0 = s

522 else:

523 t = _Fmt.SQUARE(points=i)

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

525 gc1 = gc

526

527 return True # inside

528

529 @deprecated_method

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

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

532 return self.isenclosedBy(points)

533

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

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

536

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

538 @kwarg height: Optional height for midpoint, overriding

539 the mean height (C{meter}).

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

541 may be negative or greater than 1.0.

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

543

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

545

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

547

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

549

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

551

552 @example:

553

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

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

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

557 '''

558 if fraction is _0_5:

559 self.others(other)

560

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

562 a1, b1 = self.philam

563 a2, b2 = other.philam

564

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

566

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

568

569 x = ca2 * cdb + ca1

570 y = ca2 * sdb

571

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

573 b = atan2(y, x) + b1

574

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

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

577 else:

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

579 return r

580

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

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

583

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

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

586

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

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

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

590 @kwarg options: Optional keyword arguments for function

591 L{pygeodesy.equirectangular_}.

592

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

594

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

596 see function L{pygeodesy.equirectangular_}.

597

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

599 is not (yet) supported.

600

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

602

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

604

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

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

607 '''

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

609 within = _xkwds_pop(options, within=True)

610 if not within:

611 notImplemented(self, within=within)

612

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

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

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

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

617# return point1

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

619# if isnear0(d):

620# return point1 # or point2

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

622# return point2

623# f = a / d

624# if within:

625# if f > EPS1:

626# return point2

627# elif f < EPS:

628# return point1

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

630

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

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

633 **options)[0]

634

635 @deprecated_method

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

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

638

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

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

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

642 units of B{C{radius}}.

643 '''

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

645 return r.closest, r.distance

646

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

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

649

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

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

652

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

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

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

656 @kwarg options: Optional keyword arguments for function

657 L{pygeodesy.equirectangular_}.

658

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

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

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

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

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

664 function L{pygeodesy.compassAngle}.

665

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

667 see function L{pygeodesy.equirectangular_}.

668

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

670

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

672

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

674

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

676 and L{pygeodesy.nearestOn5}.

677 '''

678 if _LatLon_ in options:

679 raise _ValueError(**options)

680

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

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

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

684

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

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

687 coordinates.

688

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

690 and other keyword arguments, ignored

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

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

693 this L{Cartesian} class or specify

694 C{B{Cartesian}=None}.

695

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

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

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

699

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

701 '''

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

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

704

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

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

707 '''

708 self.others(start=start)

709 self.others(end=end)

710

711 r = Radius_(radius)

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

713

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

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

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

717 return r, x, (e - b)

718

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

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

721

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

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

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

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

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

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

728

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

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

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

732

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

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

735 '''

736 self.others(otherB=otherB)

737 self.others(otherC=otherC)

738

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

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

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

742

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

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

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

746 of three corresponding circles.

747

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

749 as B{C{radius}}).

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

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

752 B{C{radius}}).

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

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

755 B{C{radius}}).

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

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

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

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

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

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

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

763

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

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

766 the corresponding trilaterated points C{minPoint} and

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

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

769

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

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

772

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

774 respectively largest I{radial} overlap found.

775

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

777 nearest respectively farthest intersection margin.

778

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

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

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

782 largest B{C{distance#}}.

783

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

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

786 no intersection or all (near-)concentric for

787 C{B{area}=False}.

788

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

790

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

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

793 '''

794 return _trilaterate5(self, distance1,

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

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

797 area=area, radius=radius, eps=eps, wrap=wrap)

798

799

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

801

802

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

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

805 (with the pointsjoined by great circle arcs).

806

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

808 or L{BooleanGH}).

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

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

811 or C{None}.

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

813

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

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

816

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

818

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

820

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

822

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

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

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

826

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

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

829 L{ellipsoidalKarney.areaOf}.

830

831 @example:

832

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

834 >>> areaOf(b) # 8666058750.718977

835

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

837 >>> areaOf(c) # 6.18e9

838 '''

839 if _MODS.booleans.isBoolean(points):

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

841

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

843 p1 = p2 = Ps[0]

844 a1, b1 = p1.philam

845 ta1, z1 = tan_2(a1), None

846

847 A = Fsum() # mean phi

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

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

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

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

852 D = Fsum()

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

854 a2, b2 = p2.philam

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

856 ta2 = tan_2(a2)

857 A += a2

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

859 _1_0 + ta1 * ta2)

860 ta1, b1 = ta2, b2

861

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

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

864 if z1 is not None:

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

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

867 p1, z1 = p2, z2

868

869 R = fabs(E * _2_0)

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

871 R = fabs(R - PI2)

872 if radius:

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

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

875 return float(R)

876

877

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

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

880

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

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

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

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

885

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

887 '''

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

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

890 ca *= sr

891

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

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

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

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

896

897

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

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

900 # and similar logic in .ellipsoidalBaseDI._intersect3

901 a1, b1 = start.philam

902

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

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

905 else: # must be a point

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

907 hs.append(end.height)

908 a2, b2 = end.philam

909

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

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

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

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

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

915 -(b1 + b2) * _0_5)

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

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

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

919 cb12 * cb21 + sb12 * sb21,

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

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

922

923

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

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

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

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

928

929

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

931 height=None, wrap=True):

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

933

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

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

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

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

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

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

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

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

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

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

944 circle} methods.

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

946 conventionally) or C{None}.

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

948

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

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

951 point C{is} this very instance.

952

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

954

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

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

957

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

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

960 '''

961 c = _T00.others(center=center)

962 p = _T00.others(point=point)

963 try:

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

965 height=height, wrap=wrap)

966 except (TypeError, ValueError) as x:

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

968

969

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

971 LatLon=None, **LatLon_kwds):

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

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

974

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

976

977 a1, b1 = start1.philam

978 a2, b2 = start2.philam

979

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

981 r12 = vincentys_(a2, a1, db)

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

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

984

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

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

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

988

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

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

991 raise IntersectionError(_parallel_)

992 # handle domain error for equivalent longitudes,

993 # see also functions asin_safe and acos_safe at

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

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

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

997 if sin(db) > 0:

998 t12, t21 = t1, PI2 - t2

999 else:

1000 t12, t21 = PI2 - t1, t2

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

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

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

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

1005 raise IntersectionError(_infinite_)

1006 sx3 = sx1 * sx2

1007# XXX if sx3 < 0:

1008# XXX raise ValueError(_ambiguous_)

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

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

1011

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

1013 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

1017

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

1019 _N_vector_ = _MODS.nvectorBase._N_vector_

1020

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

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

1023 x = x1.cross(x2)

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

1025 raise IntersectionError(_colinear_)

1026 a, b = x.philam

1027 # choose intersection similar to sphericalNvector

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

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

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

1031

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

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

1034 intersection, LatLon, LatLon_kwds)

1035

1036

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

1038 LatLon=LatLon, **LatLon_kwds):

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

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

1041

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

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

1044 the initial bearing at the first start point

1045 (compass C{degrees360}).

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

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

1048 the initial bearing at the second start point

1049 (compass C{degrees360}).

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

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

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

1053 @kwarg LatLon: Optional class to return the intersection

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

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

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

1057

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

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

1060 height)}. An alternate intersection point might

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

1062

1063 @raise IntersectionError: Ambiguous or infinite intersection

1064 or colinear, parallel or otherwise

1065 non-intersecting paths.

1066

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

1068

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

1070

1071 @example:

1072

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

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

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

1076 '''

1077 s1 = _T00.others(start1=start1)

1078 s2 = _T00.others(start2=start2)

1079 try:

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

1081 LatLon=LatLon, **LatLon_kwds)

1082 except (TypeError, ValueError) as x:

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

1084

1085

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

1087 height=None, wrap=True,

1088 LatLon=LatLon, **LatLon_kwds):

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

1090 by a center point and a radius.

1091

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

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

1094 see B{C{radius}}).

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

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

1097 see B{C{radius}}).

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

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

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

1101 B{C{radius}}).

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

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

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

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

1106 @kwarg LatLon: Optional class to return the intersection

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

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

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

1110

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

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

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

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

1115

1116 @raise IntersectionError: Concentric, antipodal, invalid or

1117 non-intersecting circles.

1118

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

1120

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

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

1123

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

1125

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

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

1128 '''

1129 c1 = _T00.others(center1=center1)

1130 c2 = _T00.others(center2=center2)

1131 try:

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

1133 height=height, wrap=wrap,

1134 LatLon=LatLon, **LatLon_kwds)

1135 except (TypeError, ValueError) as x:

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

1137

1138

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

1140 height=None, too_d=None, wrap=True,

1141 LatLon=LatLon, **LatLon_kwds):

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

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

1144

1145 def _dest3(bearing, h):

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

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

1148 intersections2, LatLon, LatLon_kwds)

1149

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

1151 if f: # swapped

1152 c1, c2 = c2, c1 # PYCHOK swap

1153

1154 a1, b1 = c1.philam

1155 a2, b2 = c2.philam

1156

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

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

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

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

1161

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

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

1164 if r < _0_0:

1165 raise _ValueError(eps=r)

1166

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

1168 if x > max(r, EPS):

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

1170 x = sd * sr1

1171 if isnear0(x):

1172 raise _ValueError(_invalid_)

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

1174

1175 elif x < r: # PYCHOK no cover

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

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

1178

1179 if height is None: # "radical height"

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

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

1182 else:

1183 h = Height(height)

1184

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

1186 if x < EPS4: # externally ...

1187 r = _dest3(b, h)

1188 elif x > _PI_EPS4: # internally ...

1189 r = _dest3(b + PI, h)

1190 else:

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

1192 return r, r # ... abutting circles

1193

1194

1195@deprecated_function

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

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

1198 '''

1199 return ispolar(points, wrap=wrap)

1200

1201

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

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

1204 '''

1205 n = func.__name__

1206 if LatLon is None:

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

1208 else:

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

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

1211 return r

1212

1213

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

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

1216

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

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

1219 height (C{meter}).

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

1221 or C{None}.

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

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

1224

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

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

1227 is C{None}.

1228

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

1230

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

1232 '''

1233 # geographic mean

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

1235 lat, lon = m._N_vector.latlon

1236

1237 if height is None:

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

1239 else: # PYCHOK no cover

1240 h = Height(height)

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

1242

1243

1244@deprecated_function

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

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

1247

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

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

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

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

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

1253 '''

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

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

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

1257 return ll, d

1258

1259

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

1261 LatLon=LatLon, **options):

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

1263

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

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

1266

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

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

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

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

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

1272 or C{None}.

1273 @kwarg options: Optional keyword arguments for function

1274 L{pygeodesy.equirectangular_}.

1275

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

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

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

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

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

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

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

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

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

1285

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

1287 see function L{pygeodesy.equirectangular_}.

1288

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

1290

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

1292

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

1294

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

1296 '''

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

1298 LatLon=None, **options)

1299 d = degrees2m(d, radius=radius)

1300 n = nearestOn3.__name__

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

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

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

1304

1305

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

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

1308 (with great circle arcs joining the points).

1309

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

1311 or L{BooleanGH}).

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

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

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

1315

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

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

1318

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

1320

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

1322

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

1324 C{B{points}} a composite.

1325

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

1327

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

1329 and L{ellipsoidalKarney.perimeterOf}.

1330 '''

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

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

1333 a1, b1 = Ps[0].philam

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

1335 a2, b2 = p.philam

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

1337 yield vincentys_(a2, a1, db)

1338 a1, b1 = a2, b2

1339

1340 if _MODS.booleans.isBoolean(points):

1341 if not closed:

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

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

1344 else:

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

1346 return _r2m(r, radius)

1347

1348

1349 def _r2m(r, radius):

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

1351 '''

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

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

1354 return r

1355

1356

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

1358 excess=excessAbc_,

1359 wrap=False):

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

1361

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

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

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

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

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

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

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

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

1370 or C{None}.

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

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

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

1374

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

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

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

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

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

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

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

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

1383 '''

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

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

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

1387 excess=excess, wrap=wrap)

1388 return _t7Tuple(t, radius)

1389

1390

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

1392 wrap=False):

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

1394 excess} of a (spherical) triangle.

1395

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

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

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

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

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

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

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

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

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

1405

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

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

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

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

1410 '''

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

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

1413 a = vincentys_(phiC, phiB, d)

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

1415

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

1417 s = sb * sc

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

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

1420

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

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

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

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

1425

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

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

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

1429

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

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

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

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

1434

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

1436

1437

1438 def _t7Tuple(t, radius):

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

1440 '''

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

1442 r = radius if isscalar(radius) else \

1443 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1446 B, (r * t.b),

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

1448 return t

1449

1450

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

1452 areaOf, # functions

1453 intersecant2, intersection, intersections2, ispolar,

1454 isPoleEnclosedBy, # DEPRECATED, use ispolar

1455 meanOf,

1456 nearestOn2, nearestOn3,

1457 perimeterOf,

1458 sumOf, # == vector3d.sumOf

1459 triangle7, triangle8_)

1460

1461# **) MIT License

1462#

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

1464#

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

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

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

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

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

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

1471#

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

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

1474#

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

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

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

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

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

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

1481# OTHER DEALINGS IN THE SOFTWARE.