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

34 _LatLon_, _line_, _near_, _not_, _null_, _points_, \

35 _SPACE_, _too_

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

37# from pygeodesy.named import notImplemented # from .points

38 from pygeodesy.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.11'

58

59_parallel_ = 'parallel'

60

61_PI_EPS4 = PI - EPS4

62 if _PI_EPS4 >= PI:

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

64

65

66 class Cartesian(CartesianSphericalBase):

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

68 spherical, geodetic L{LatLon}.

69 '''

70

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

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

73

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

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

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

77

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

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

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

81

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

83 '''

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

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

86

87

88 class LatLon(LatLonSphericalBase):

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

90

91 @example:

92

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

94 '''

95

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

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

98 the closest point on the great circle line defined by a

99 start and an end point.

100

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

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

103 from the start point to the point where the perpendicular

104 crosses the line.

105

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

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

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

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

110

111 @return: Distance along the great circle line (C{meter},

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

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

114 toward the B{C{end}} point of the line or negative

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

116

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

118

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

120

121 @example:

122

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

124

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

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

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

128 '''

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

130 cx = cos(x)

131 return _0_0 if isnear0(cx) else \

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

133

134 @deprecated_method

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

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

137 '''

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

139

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

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

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

143

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

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

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

147

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

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

150 '''

151 self.others(other)

152

153 a1, b1 = self.philam

154 a2, b2 = other.philam

155

156 a = radians(lat)

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

158

159 sa, ca, sa1, ca1, \

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

161 sa1 *= ca2 * ca

162

163 x = sa1 * sdb

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

165 z = ca1 * sdb * ca2 * sa

166

167 h = hypot(x, y)

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

169 return None # great circle doesn't reach latitude

170

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

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

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

174

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

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

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

178

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

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

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

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

183

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

185 or positive if to the right of the line) (C{meter},

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

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

188

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

190

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

192

193 @example:

194

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

196

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

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

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

200 '''

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

202 return _r2m(x, radius)

203

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

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

206 travelled the given distance on the given initial bearing.

207

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

209 B{C{radius}}).

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

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

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

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

214

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

216

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

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

219

220 @example:

221

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

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

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

225 '''

226 a, b = self.philam

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

228

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

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

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

232

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

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

235

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

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

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

239

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

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

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

243

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

245

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

247

248 @example:

249

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

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

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

253 '''

254 self.others(other)

255

256 a1, b1 = self.philam

257 a2, b2 = other.philam

258

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

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

261

262# @Property_RO

263# def Ecef(self):

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

265# '''

266# return _MODS.ecef.EcefKarney

267

268 def greatCircle(self, bearing):

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

270 on the given initial bearing from this point.

271

272 Direction of vector is such that initial bearing vector

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

274

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

276

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

278

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

280

281 @example:

282

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

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

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

286 '''

287 a, b = self.philam

288

289 t = Bearing_(bearing)

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

291

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

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

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

295

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

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

298 to an other point.

299

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

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

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

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

304 C{sphericalNvector.LatLon.initialBearingTo}.

305

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

307

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

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

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

311

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

313

314 @example:

315

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

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

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

319 '''

320 self.others(other)

321

322 a1, b1 = self.philam

323 a2, b2 = other.philam

324

325 # XXX behavior like sphericalNvector.LatLon.initialBearingTo

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

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

328

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

330

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

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

333 and an other point.

334

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

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

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

338 @kwarg height: Optional height, overriding the intermediate

339 height (C{meter}).

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

341

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

343

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

345

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

347

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

349

350 @example:

351

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

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

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

355 '''

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

357 if isnear0(f): # PYCHOK no cover

358 r = self

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

360 r = self.others(other)

361 else:

362 a1, b1 = self.philam

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

364

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

366 r = vincentys_(a2, a1, db)

367 sr = sin(r)

368 if isnon0(sr):

369 sa1, ca1, sa2, ca2, \

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

371

372 t = f * r

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

374 b = sin( t) # / sr superflous

375

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

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

378 z = a * sa1 + b * sa2

379

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

381 b = atan2(y, x)

382

383 else: # PYCHOK no cover

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

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

386

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

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

389 return r

390

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

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

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

394

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

396 initial bearing at this point (compass

397 C{degrees360}).

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

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

400 the initial bearing at the other start point

401 (compass C{degrees360}).

402 @kwarg height: Optional height for intersection point,

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

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

405

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

407 intersection point might be the L{antipode} to

408 the returned result.

409

410 @raise IntersectionError: Ambiguous or infinite intersection

411 or colinear, parallel or otherwise

412 non-intersecting lines.

413

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

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

416

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

418

419 @example:

420

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

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

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

424 '''

425 try:

426 s2 = self.others(other)

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

428 LatLon=self.classof)

429 except (TypeError, ValueError) as x:

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

431

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

433 height=None, wrap=True):

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

435 by a center point and radius.

436

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

438 see B{C{radius}}).

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

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

441 see B{C{radius}}).

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

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

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

445 B{C{radius}}).

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

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

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

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

450

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

452 instance. For abutting circles, both intersection

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

454

455 @raise IntersectionError: Concentric, antipodal, invalid or

456 non-intersecting circles.

457

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

459

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

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

462 '''

463 try:

464 c2 = self.others(other)

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

466 height=height, wrap=wrap,

467 LatLon=self.classof)

468 except (TypeError, ValueError) as x:

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

470

471 def isenclosedBy(self, points):

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

473

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

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

476

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

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 I{longitudinally}.

489 '''

490 if _MODS.booleans.isBoolean(points):

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

492

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

494 n0 = self._N_vector

495

496 v2 = Ps[0]._N_vector

497 v1 = Ps[1]._N_vector

498 gc1 = v2.cross(v1)

499 # check whether this point on same side of all

500 # polygon edges (to the left or right depending

501 # on anti-/clockwise polygon direction)

502 t0 = gc1.angleTo(n0) > PI_2

503 s0 = None

504 # get great-circle vector for each edge

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

506 v2 = p._N_vector

507 gc = v1.cross(v2)

508 v1 = v2

509

510 t = gc.angleTo(n0) > PI_2

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

512 return False # outside

513

514 # check for convex polygon: angle between

515 # gc vectors, signed by direction of n0

516 # (otherwise the test above is not reliable)

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

518 if s != s0:

519 if s0 is None:

520 s0 = s

521 else:

522 t = _Fmt.SQUARE(points=i)

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

524 gc1 = gc

525

526 return True # inside

527

528 @deprecated_method

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

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

531 return self.isenclosedBy(points)

532

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

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

535

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

537 @kwarg height: Optional height for midpoint, overriding

538 the mean height (C{meter}).

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

540 may be negative or greater than 1.0.

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

542

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

544

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

546

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

548

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

550

551 @example:

552

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

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

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

556 '''

557 if fraction is _0_5:

558 self.others(other)

559

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

561 a1, b1 = self.philam

562 a2, b2 = other.philam

563

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

565

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

567

568 x = ca2 * cdb + ca1

569 y = ca2 * sdb

570

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

572 b = atan2(y, x) + b1

573

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

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

576 else:

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

578 return r

579

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

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

582

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

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

585

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

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

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

589 @kwarg options: Optional keyword arguments for function

590 L{pygeodesy.equirectangular_}.

591

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

593

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

595 see function L{pygeodesy.equirectangular_}.

596

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

598 is not (yet) supported.

599

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

601

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

603

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

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

606 '''

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

608 within = _xkwds_pop(options, within=True)

609 if not within:

610 notImplemented(self, within=within)

611

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

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

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

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

616# return point1

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

618# if isnear0(d):

619# return point1 # or point2

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

621# return point2

622# f = a / d

623# if within:

624# if f > EPS1:

625# return point2

626# elif f < EPS:

627# return point1

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

629

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

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

632 **options)[0]

633

634 @deprecated_method

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

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

637

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

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

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

641 units of B{C{radius}}.

642 '''

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

644 return r.closest, r.distance

645

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

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

648

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

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

651

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

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

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

655 @kwarg options: Optional keyword arguments for function

656 L{pygeodesy.equirectangular_}.

657

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

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

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

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

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

663 function L{pygeodesy.compassAngle}.

664

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

666 see function L{pygeodesy.equirectangular_}.

667

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

669

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

671

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

673

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

675 and L{pygeodesy.nearestOn5}.

676 '''

677 if _LatLon_ in options:

678 raise _ValueError(**options)

679

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

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

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

683

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

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

686 coordinates.

687

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

689 and other keyword arguments, ignored

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

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

692 this L{Cartesian} class or specify

693 C{B{Cartesian}=None}.

694

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

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

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

698

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

700 '''

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

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

703

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

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

706 '''

707 self.others(start=start)

708 self.others(end=end)

709

710 r = Radius_(radius)

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

712

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

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

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

716 return r, x, (e - b)

717

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

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

720

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

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

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

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

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

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

727

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

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

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

731

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

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

734 '''

735 self.others(otherB=otherB)

736 self.others(otherC=otherC)

737

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

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

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

741

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

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

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

745 of three corresponding circles.

746

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

748 as B{C{radius}}).

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

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

751 B{C{radius}}).

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

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

754 B{C{radius}}).

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

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

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

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

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

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

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

762

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

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

765 the corresponding trilaterated points C{minPoint} and

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

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

768

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

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

771

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

773 respectively largest I{radial} overlap found.

774

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

776 nearest respectively farthest intersection margin.

777

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

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

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

781 largest B{C{distance#}}.

782

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

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

785 no intersection or all (near-)concentric for

786 C{B{area}=False}.

787

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

789

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

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

792 '''

793 return _trilaterate5(self, distance1,

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

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

796 area=area, radius=radius, eps=eps, wrap=wrap)

797

798

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

800

801

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

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

804 (with the pointsjoined by great circle arcs).

805

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

807 or L{BooleanGH}).

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

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

810 or C{None}.

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

812

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

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

815

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

817

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

819

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

821

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

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

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

825

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

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

828 L{ellipsoidalKarney.areaOf}.

829

830 @example:

831

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

833 >>> areaOf(b) # 8666058750.718977

834

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

836 >>> areaOf(c) # 6.18e9

837 '''

838 if _MODS.booleans.isBoolean(points):

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

840

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

842 p1 = p2 = Ps[0]

843 a1, b1 = p1.philam

844 ta1, z1 = tan_2(a1), None

845

846 A = Fsum() # mean phi

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

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

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

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

851 D = Fsum()

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

853 a2, b2 = p2.philam

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

855 ta2 = tan_2(a2)

856 A += a2

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

858 _1_0 + ta1 * ta2)

859 ta1, b1 = ta2, b2

860

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

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

863 if z1 is not None:

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

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

866 p1, z1 = p2, z2

867

868 R = fabs(E * _2_0)

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

870 R = fabs(R - PI2)

871 if radius:

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

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

874 return float(R)

875

876

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

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

879

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

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

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

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

884

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

886 '''

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

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

889 ca *= sr

890

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

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

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

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

895

896

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

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

899 # and similar logic in .ellipsoidalBaseDI._intersect3

900 a1, b1 = start.philam

901

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

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

904 else: # must be a point

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

906 hs.append(end.height)

907 a2, b2 = end.philam

908

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

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

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

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

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

914 -(b1 + b2) * _0_5)

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

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

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

918 cb12 * cb21 + sb12 * sb21,

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

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

921

922

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

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

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

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

927

928

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

930 height=None, wrap=True):

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

932

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

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

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

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

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

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

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

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

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

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

943 circle} methods.

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

945 conventionally) or C{None}.

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

947

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

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

950 point C{is} this very instance.

951

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

953

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

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

956

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

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

959 '''

960 c = _T00.others(center=center)

961 p = _T00.others(point=point)

962 try:

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

964 height=height, wrap=wrap)

965 except (TypeError, ValueError) as x:

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

967

968

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

970 LatLon=None, **LatLon_kwds):

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

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

973

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

975

976 a1, b1 = start1.philam

977 a2, b2 = start2.philam

978

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

980 r12 = vincentys_(a2, a1, db)

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

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

983

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

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

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

987

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

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

990 raise IntersectionError(_parallel_)

991 # handle domain error for equivalent longitudes,

992 # see also functions asin_safe and acos_safe at

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

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

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

996 if sin(db) > 0:

997 t12, t21 = t1, PI2 - t2

998 else:

999 t12, t21 = PI2 - t1, t2

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

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

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

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

1004 raise IntersectionError(_infinite_)

1005 sx3 = sx1 * sx2

1006# XXX if sx3 < 0:

1007# XXX raise ValueError(_ambiguous_)

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

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

1010

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

1012 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

1016

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

1018 _N_vector_ = _MODS.nvectorBase._N_vector_

1019

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

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

1022 x = x1.cross(x2)

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

1024 raise IntersectionError(_colinear_)

1025 a, b = x.philam

1026 # choose intersection similar to sphericalNvector

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

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

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

1030

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

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

1033 intersection, LatLon, LatLon_kwds)

1034

1035

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

1037 LatLon=LatLon, **LatLon_kwds):

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

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

1040

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

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

1043 the initial bearing at the first start point

1044 (compass C{degrees360}).

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

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

1047 the initial bearing at the second start point

1048 (compass C{degrees360}).

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

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

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

1052 @kwarg LatLon: Optional class to return the intersection

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

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

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

1056

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

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

1059 height)}. An alternate intersection point might

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

1061

1062 @raise IntersectionError: Ambiguous or infinite intersection

1063 or colinear, parallel or otherwise

1064 non-intersecting lines.

1065

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

1067

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

1069

1070 @example:

1071

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

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

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

1075 '''

1076 s1 = _T00.others(start1=start1)

1077 s2 = _T00.others(start2=start2)

1078 try:

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

1080 LatLon=LatLon, **LatLon_kwds)

1081 except (TypeError, ValueError) as x:

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

1083

1084

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

1086 height=None, wrap=True,

1087 LatLon=LatLon, **LatLon_kwds):

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

1089 by a center point and a radius.

1090

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

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

1093 see B{C{radius}}).

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

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

1096 see B{C{radius}}).

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

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

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

1100 B{C{radius}}).

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

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

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

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

1105 @kwarg LatLon: Optional class to return the intersection

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

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

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

1109

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

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

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

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

1114

1115 @raise IntersectionError: Concentric, antipodal, invalid or

1116 non-intersecting circles.

1117

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

1119

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

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

1122

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

1124

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

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

1127 '''

1128 c1 = _T00.others(center1=center1)

1129 c2 = _T00.others(center2=center2)

1130 try:

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

1132 height=height, wrap=wrap,

1133 LatLon=LatLon, **LatLon_kwds)

1134 except (TypeError, ValueError) as x:

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

1136

1137

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

1139 height=None, too_d=None, wrap=True,

1140 LatLon=LatLon, **LatLon_kwds):

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

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

1143

1144 def _dest3(bearing, h):

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

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

1147 intersections2, LatLon, LatLon_kwds)

1148

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

1150 if f: # swapped

1151 c1, c2 = c2, c1 # PYCHOK swap

1152

1153 a1, b1 = c1.philam

1154 a2, b2 = c2.philam

1155

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

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

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

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

1160

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

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

1163 if r < _0_0:

1164 raise _ValueError(eps=r)

1165

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

1167 if x > max(r, EPS):

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

1169 x = sd * sr1

1170 if isnear0(x):

1171 raise _ValueError(_invalid_)

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

1173

1174 elif x < r: # PYCHOK no cover

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

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

1177

1178 if height is None: # "radical height"

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

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

1181 else:

1182 h = Height(height)

1183

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

1185 if x < EPS4: # externally ...

1186 r = _dest3(b, h)

1187 elif x > _PI_EPS4: # internally ...

1188 r = _dest3(b + PI, h)

1189 else:

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

1191 return r, r # ... abutting circles

1192

1193

1194@deprecated_function

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

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

1197 '''

1198 return ispolar(points, wrap=wrap)

1199

1200

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

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

1203 '''

1204 n = func.__name__

1205 if LatLon is None:

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

1207 else:

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

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

1210 return r

1211

1212

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

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

1215

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

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

1218 height (C{meter}).

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

1220 or C{None}.

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

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

1223

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

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

1226 is C{None}.

1227

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

1229

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

1231 '''

1232 # geographic mean

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

1234 lat, lon = m._N_vector.latlon

1235

1236 if height is None:

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

1238 else: # PYCHOK no cover

1239 h = Height(height)

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

1241

1242

1243@deprecated_function

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

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

1246

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

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

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

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

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

1252 '''

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

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

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

1256 return ll, d

1257

1258

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

1260 LatLon=LatLon, **options):

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

1262

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

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

1265

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

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

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

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

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

1271 or C{None}.

1272 @kwarg options: Optional keyword arguments for function

1273 L{pygeodesy.equirectangular_}.

1274

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

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

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

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

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

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

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

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

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

1284

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

1286 see function L{pygeodesy.equirectangular_}.

1287

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

1289

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

1291

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

1293

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

1295 '''

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

1297 LatLon=None, **options)

1298 d = degrees2m(d, radius=radius)

1299 n = nearestOn3.__name__

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

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

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

1303

1304

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

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

1307 (with great circle arcs joining the points).

1308

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

1310 or L{BooleanGH}).

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

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

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

1314

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

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

1317

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

1319

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

1321

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

1323 C{B{points}} a composite.

1324

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

1326

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

1328 and L{ellipsoidalKarney.perimeterOf}.

1329 '''

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

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

1332 a1, b1 = Ps[0].philam

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

1334 a2, b2 = p.philam

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

1336 yield vincentys_(a2, a1, db)

1337 a1, b1 = a2, b2

1338

1339 if _MODS.booleans.isBoolean(points):

1340 if not closed:

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

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

1343 else:

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

1345 return _r2m(r, radius)

1346

1347

1348 def _r2m(r, radius):

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

1350 '''

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

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

1353 return r

1354

1355

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

1357 excess=excessAbc_,

1358 wrap=False):

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

1360

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

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

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

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

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

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

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

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

1369 or C{None}.

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

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

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

1373

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

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

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

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

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

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

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

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

1382 '''

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

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

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

1386 excess=excess, wrap=wrap)

1387 return _t7Tuple(t, radius)

1388

1389

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

1391 wrap=False):

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

1393 excess} of a (spherical) triangle.

1394

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

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

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

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

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

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

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

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

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

1404

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

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

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

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

1409 '''

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

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

1412 a = vincentys_(phiC, phiB, d)

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

1414

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

1416 s = sb * sc

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

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

1419

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

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

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

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

1424

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

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

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

1428

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

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

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

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

1433

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

1435

1436

1437 def _t7Tuple(t, radius):

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

1439 '''

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

1441 r = radius if isscalar(radius) else \

1442 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1445 B, (r * t.b),

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

1447 return t

1448

1449

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

1451 areaOf, # functions

1452 intersecant2, intersection, intersections2, ispolar,

1453 isPoleEnclosedBy, # DEPRECATED, use ispolar

1454 meanOf,

1455 nearestOn2, nearestOn3,

1456 perimeterOf,

1457 sumOf, # == vector3d.sumOf

1458 triangle7, triangle8_)

1459

1460# **) MIT License

1461#

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

1463#

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

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

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

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

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

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

1470#

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

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

1473#

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

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

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

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

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

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

1480# OTHER DEALINGS IN THE SOFTWARE.