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

29 from pygeodesy.formy import antipode_, bearing_, _bearingTo2, excessAbc_, \

30 excessGirard_, excessLHuilier_, opposing_, _radical2, \

31 vincentys_

32 from pygeodesy.interns import _1_, _2_, _coincident_, _composite_, _colinear_, \

33 _concentric_, _convex_, _end_, _infinite_, \

34 _invalid_, _line_, _near_, _not_, _null_, \

35 _point_, _SPACE_, _too_

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

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

38# from pygeodesy.nvectorBase import NvectorBase, sumOf # _MODE

39 from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \

40 NearestOn3Tuple, Triangle7Tuple, \

41 Triangle8Tuple

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

43 notImplemented, Fmt as _Fmt # XXX shadowed

44 from pygeodesy.props import deprecated_function, deprecated_method

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

46 LatLonSphericalBase, _rads3, _trilaterate5

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

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

49 Radius_, Scalar

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

51 m2radians, radiansPI2, sincos2_, tan_2, _unrollon, \

52 unrollPI, _unrollon3, _Wrap, wrap180, wrapPI

53 from pygeodesy.vector3d import sumOf, Vector3d

54

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

56

57__all__ = _ALL_LAZY.sphericalTrigonometry

58__version__ = '23.05.25'

59

60_parallel_ = 'parallel'

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 _ab1_ab2_db5(self, other, wrap):

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

99 '''

100 a1, b1 = self.philam

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

102 if wrap:

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

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

105 else: # unrollPI shortcut

106 db = b2 - b1

107 return a1, b1, a2, b2, db

108

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

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

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

112 start and an end point.

113

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

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

116 from the start point to the point where the perpendicular

117 crosses the line.

118

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

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

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

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

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

124

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

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

127 C{B{radius} is None}, positive if I{after} the

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

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

130 B{C{start}} point.

131

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

133

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

135

136 @example:

137

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

139

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

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

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

143 '''

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

145 cx = cos(x)

146 return _0_0 if isnear0(cx) else \

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

148

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

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

151 '''

152 s = self.others(start=start)

153 e = self.others(end=end)

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

155

156 r = Radius_(radius)

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

158

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

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

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

162 return r, x, -b

163

164 @deprecated_method

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

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

167 '''

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

169

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

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

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

173

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

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

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

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

178

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

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

181 '''

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

183 sa, ca, sa1, ca1, \

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

185 sa1 *= ca2 * ca

186

187 x = sa1 * sdb

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

189 z = ca1 * sdb * ca2 * sa

190

191 h = hypot(x, y)

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

193 return None # great circle doesn't reach latitude

194

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

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

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

198

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

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

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

202

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

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

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

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

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

208

209 @return: Distance to the great circle (I{negative} if to

210 the left or I{positive} if to the right of the

211 line) (C{meter}, same units as B{C{radius}} or

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

213

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

215

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

217

218 @example:

219

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

221

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

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

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

225 '''

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

227 return _r2m(x, radius)

228

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

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

231 travelled the given distance on the given initial bearing.

232

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

234 B{C{radius}}).

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

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

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

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

239

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

241

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

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

244

245 @example:

246

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

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

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

250 '''

251 a, b = self.philam

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

253

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

255 h = self._heigHt(height)

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

257

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

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

260

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

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

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

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

265

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

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

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

269

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

271

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

273

274 @example:

275

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

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

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

279 '''

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

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

282

283# @Property_RO

284# def Ecef(self):

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

286# '''

287# return _MODS.ecef.EcefKarney

288

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

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

291 on the given initial bearing from this point.

292

293 Direction of vector is such that initial bearing vector

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

295

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

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

298 overriding the default L{Vector3d}.

299 @kwarg Vector_kwds: Optional, additional keyword argunents

300 for B{C{Vector}}.

301

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

303

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

305

306 @example:

307

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

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

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

311 '''

312 a, b = self.philam

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

314

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

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

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

318

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

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

321 to an other point.

322

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

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

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

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

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

328 C{sphericalNvector.LatLon.initialBearingTo}.

329

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

331

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

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

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

335

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

337

338 @example:

339

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

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

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

343 '''

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

345 # XXX behavior like sphericalNvector.LatLon.initialBearingTo

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

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

348

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

350

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

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

353 and an other point.

354

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

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

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

358 @kwarg height: Optional height, overriding the intermediate

359 height (C{meter}).

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

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

362

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

364

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

366

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

368

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

370

371 @example:

372

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

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

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

376 '''

377 p = self

378 f = Scalar(fraction=fraction)

379 if not isnear0(f):

380 p = p.others(other)

381 if wrap:

382 p = _Wrap.point(p)

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

384 a1, b1 = self.philam

385 a2, b2 = p.philam

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

387 r = vincentys_(a2, a1, db)

388 sr = sin(r)

389 if isnon0(sr):

390 sa1, ca1, sa2, ca2, \

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

392

393 t = f * r

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

395 b = sin( t) # / sr superflous

396

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

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

399 z = a * sa1 + b * sa2

400

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

402 b = atan2(y, x)

403

404 else: # PYCHOK no cover

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

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

407

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

409 p = self.classof(degrees90(a), degrees180(b), height=h)

410 return p

411

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

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

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

415

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

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

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

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

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

421 C{degrees360}).

422 @kwarg height: Optional height for intersection point,

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

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

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

426

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

428 intersection point might be the L{antipode} to

429 the returned result.

430

431 @raise IntersectionError: Ambiguous or infinite intersection

432 or colinear, parallel or otherwise

433 non-intersecting lines.

434

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

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

437

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

439

440 @example:

441

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

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

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

445 '''

446 try:

447 s2 = self.others(other)

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

449 LatLon=self.classof)

450 except (TypeError, ValueError) as x:

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

452 other=other, end2=end2, wrap=wrap)

453

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

455 height=None, wrap=True):

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

457 by a center point and radius.

458

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

460 see B{C{radius}}).

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

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

463 see B{C{radius}}).

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

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

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

467 B{C{radius}}).

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

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

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

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

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

473

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

475 instance. For abutting circles, both intersection

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

477

478 @raise IntersectionError: Concentric, antipodal, invalid or

479 non-intersecting circles.

480

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

482

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

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

485 '''

486 try:

487 c2 = self.others(other)

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

489 height=height, wrap=wrap,

490 LatLon=self.classof)

491 except (TypeError, ValueError) as x:

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

493 other=other, rad2=rad2, wrap=wrap)

494

495 @deprecated_method

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

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

498 return self.isenclosedBy(points)

499

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

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

502

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

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

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

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

507

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

509 composite, C{False} otherwise.

510

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

512

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

514

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

516

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

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

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

520 '''

521 if _MODS.booleans.isBoolean(points):

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

523

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

525 n0 = self._N_vector

526

527 v2 = Ps[0]._N_vector

528 p1 = Ps[1]

529 v1 = p1._N_vector

530 # check whether this point on same side of all

531 # polygon edges (to the left or right depending

532 # on the anti-/clockwise polygon direction)

533 gc1 = v2.cross(v1)

534 t0 = gc1.angleTo(n0) > PI_2

535 s0 = None

536 # get great-circle vector for each edge

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

538 if wrap and not Ps.looped:

539 p2 = _unrollon(p1, p2)

540 p1 = p2

541 v2 = p2._N_vector

542 gc = v1.cross(v2)

543 t = gc.angleTo(n0) > PI_2

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

545 return False # outside

546

547 # check for convex polygon: angle between

548 # gc vectors, signed by direction of n0

549 # (otherwise the test above is not reliable)

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

551 if s != s0:

552 if s0 is None:

553 s0 = s

554 else:

555 t = _Fmt.SQUARE(points=i)

556 raise _ValueError(t, p2, wrap=wrap, txt=_not_(_convex_))

557 gc1, v1 = gc, v2

558

559 return True # inside

560

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

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

563

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

565 @kwarg height: Optional height for midpoint, overriding

566 the mean height (C{meter}).

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

568 may be negative or greater than 1.0.

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

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

571

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

573

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

575

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

577

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

579

580 @example:

581

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

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

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

585 '''

586 if fraction is _0_5:

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

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

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

590

591 x = ca2 * cdb + ca1

592 y = ca2 * sdb

593

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

595 b += atan2(y, x)

596

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

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

599 else:

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

601 return r

602

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

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

605

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

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

608

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

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

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

612 @kwarg wrap_adjust_limit: Optional keyword arguments for functions

613 L{sphericalTrigonometry.nearestOn3} and

614 L{pygeodesy.equirectangular_},

615

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

617

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

619 see function L{pygeodesy.equirectangular_}.

620

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

622 is not (yet) supported.

623

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

625

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

627

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

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

630 '''

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

632 w = _xkwds_pop(wrap_adjust_limit, within=True)

633 if not w:

634 notImplemented(self, within=w)

635

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

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

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

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

640# return point1

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

642# if isnear0(d):

643# return point1 # or point2

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

645# return point2

646# f = a / d

647# if within:

648# if f > EPS1:

649# return point2

650# elif f < EPS:

651# return point1

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

653

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

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

656 **wrap_adjust_limit)[0]

657

658 @deprecated_method

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

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

661

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

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

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

665 units of B{C{radius}}.

666 '''

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

668 return r.closest, r.distance

669

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

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

672

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

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

675

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

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

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

679 @kwarg wrap_adjust_limit: Optional keyword arguments for function

680 L{sphericalTrigonometry.nearestOn3} and

681 L{pygeodesy.equirectangular_},

682

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

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

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

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

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

688 function L{pygeodesy.compassAngle}.

689

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

691 see function L{pygeodesy.equirectangular_}.

692

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

694

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

696

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

698

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

700 and L{pygeodesy.nearestOn5}.

701 '''

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

703 LatLon=self.classof, **wrap_adjust_limit)

704

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

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

707 coordinates.

708

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

710 and other keyword arguments, ignored

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

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

713 this L{Cartesian} class or specify

714 C{B{Cartesian}=None}.

715

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

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

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

719

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

721 '''

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

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

724

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

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

727

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

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

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

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

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

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

734 B{C{otherB}} and B{C{otherC}} points (C{bool}).

735

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

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

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

739

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

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

742 '''

743 B = self.others(otherB=otherB)

744 C = self.others(otherC=otherC)

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

746

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

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

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

750

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

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

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

754 of three corresponding circles.

755

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

757 as B{C{radius}}).

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

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

760 B{C{radius}}).

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

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

763 B{C{radius}}).

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

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

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

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

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

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

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

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

772

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

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

775 the corresponding trilaterated points C{minPoint} and

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

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

778

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

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

781

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

783 respectively largest I{radial} overlap found.

784

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

786 nearest respectively farthest intersection margin.

787

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

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

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

791 largest B{C{distance#}}.

792

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

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

795 no intersection or all (near-)concentric for

796 C{B{area}=False}.

797

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

799

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

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

802 or B{C{radius}}.

803 '''

804 return _trilaterate5(self, distance1,

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

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

807 area=area, radius=radius, eps=eps, wrap=wrap)

808

809

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

811

812

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

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

815 (with the pointsjoined by great circle arcs).

816

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

818 or L{BooleanGH}).

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

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

821 or C{None}.

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

823 (C{bool}).

824

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

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

827

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

829

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

831

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

833

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

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

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

837

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

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

840 L{ellipsoidalKarney.areaOf}.

841

842 @example:

843

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

845 >>> areaOf(b) # 8666058750.718977

846

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

848 >>> areaOf(c) # 6.18e9

849 '''

850 if _MODS.booleans.isBoolean(points):

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

852

853 _at2, _t_2 = atan2, tan_2

854 _un, _w180 = unrollPI, wrap180

855

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

857 p1 = p2 = Ps[0]

858 a1, b1 = p1.philam

859 ta1, z1 = _t_2(a1), None

860

861 A = Fsum() # mean phi

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

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

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

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

866 D = Fsum()

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

868 a2, b2 = p2.philam

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

870 A += a2

871 ta2 = _t_2(a2)

872 tdb = _t_2(db, points=i)

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

874 _1_0 + ta1 * ta2)

875 ta1, b1 = ta2, b2

876

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

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

879 if z1 is not None:

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

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

882 p1, z1 = p2, z2

883

884 R = abs(R * _2_0)

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

886 R = abs(R - PI2)

887 if radius:

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

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

890 return float(R)

891

892

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

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

895

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

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

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

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

900

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

902 '''

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

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

905 ca *= sr

906

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

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

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

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

911

912

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

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

915 # and similar logic in .ellipsoidalBaseDI._intersect3

916 a1, b1 = s.philam

917

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

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

920 else: # must be a point

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

922 hs.append(end.height)

923 a2, b2 = end.philam

924 if wrap:

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

926

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

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

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

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

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

932 -(b1 + b2) * _0_5)

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

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

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

936 cb12 * cb21 + sb12 * sb21,

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

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

939

940

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

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

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

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

945

946

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

948 height=None, wrap=False): # was=True

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

950

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

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

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

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

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

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

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

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

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

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

961 circle} methods.

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

963 conventionally) or C{None}.

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

965 and the B{C{circle}} and B{C{bearing}} points (C{bool}).

966

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

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

969 point C{is} this very instance.

970

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

972

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

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

975

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

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

978 '''

979 c = _T00.others(center=center)

980 p = _T00.others(point=point)

981 try:

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

983 height=height, wrap=wrap)

984 except (TypeError, ValueError) as x:

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

986

987

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

989 LatLon=None, **LatLon_kwds):

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

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

992

993 s1, s2 = start1, start2

994 if wrap:

995 s2 = _Wrap.point(s2)

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

997

998 a1, b1 = s1.philam

999 a2, b2 = s2.philam

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

1001 r12 = vincentys_(a2, a1, db)

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

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

1004

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

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

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

1008

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

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

1011 raise IntersectionError(_parallel_)

1012 # handle domain error for equivalent longitudes,

1013 # see also functions asin_safe and acos_safe at

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

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

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

1017 if sin(db) > 0:

1018 t12, t21 = t1, PI2 - t2

1019 else:

1020 t12, t21 = PI2 - t1, t2

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

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

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

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

1025 raise IntersectionError(_infinite_)

1026 sx3 = sx1 * sx2

1027# XXX if sx3 < 0:

1028# XXX raise ValueError(_ambiguous_)

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

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

1031

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

1033 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

1037

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

1039 _N_vector_ = _MODS.nvectorBase._N_vector_

1040

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

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

1043 x = x1.cross(x2)

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

1045 raise IntersectionError(_colinear_)

1046 a, b = x.philam

1047 # choose intersection similar to sphericalNvector

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

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

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

1051

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

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

1054 intersection, LatLon, LatLon_kwds)

1055

1056

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

1058 LatLon=LatLon, **LatLon_kwds):

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

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

1061

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

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

1064 the initial bearing at the first start point

1065 (compass C{degrees360}).

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

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

1068 the initial bearing at the second start point

1069 (compass C{degrees360}).

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

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

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

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

1074 @kwarg LatLon: Optional class to return the intersection

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

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

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

1078

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

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

1081 height)}. An alternate intersection point might

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

1083

1084 @raise IntersectionError: Ambiguous or infinite intersection

1085 or colinear, parallel or otherwise

1086 non-intersecting lines.

1087

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

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

1090

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

1092

1093 @example:

1094

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

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

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

1098 '''

1099 s1 = _T00.others(start1=start1)

1100 s2 = _T00.others(start2=start2)

1101 try:

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

1103 LatLon=LatLon, **LatLon_kwds)

1104 except (TypeError, ValueError) as x:

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

1106

1107

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

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

1110 LatLon=LatLon, **LatLon_kwds):

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

1112 by a center point and a radius.

1113

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

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

1116 see B{C{radius}}).

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

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

1119 see B{C{radius}}).

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

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

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

1123 B{C{radius}}).

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

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

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

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

1128 (C{bool}).

1129 @kwarg LatLon: Optional class to return the intersection

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

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

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

1133

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

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

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

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

1138

1139 @raise IntersectionError: Concentric, antipodal, invalid or

1140 non-intersecting circles.

1141

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

1143

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

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

1146

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

1148

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

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

1151 '''

1152 c1 = _T00.others(center1=center1)

1153 c2 = _T00.others(center2=center2)

1154 try:

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

1156 height=height, wrap=wrap,

1157 LatLon=LatLon, **LatLon_kwds)

1158 except (TypeError, ValueError) as x:

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

1160 center2=center2, rad2=rad2, wrap=wrap)

1161

1162

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

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

1165 LatLon=LatLon, **LatLon_kwds):

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

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

1168

1169 def _dest3(bearing, h):

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

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

1172 intersections2, LatLon, LatLon_kwds)

1173

1174 a1, b1 = c1.philam

1175 a2, b2 = c2.philam

1176 if wrap:

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

1178

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

1180 if f: # swapped radii, swap centers

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

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

1183

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

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

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

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

1188

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

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

1191 if r < _0_0:

1192 raise _ValueError(eps=r)

1193

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

1195 if x > max(r, EPS):

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

1197 x = sd * sr1

1198 if isnear0(x):

1199 raise _ValueError(_invalid_)

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

1201

1202 elif x < r: # PYCHOK no cover

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

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

1205

1206 if height is None: # "radical height"

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

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

1209 else:

1210 h = Height(height)

1211

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

1213 if x < EPS4: # externally ...

1214 r = _dest3(b, h)

1215 elif x > _PI_EPS4: # internally ...

1216 r = _dest3(b + PI, h)

1217 else:

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

1219 return r, r # ... abutting circles

1220

1221

1222@deprecated_function

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

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

1225 '''

1226 return ispolar(points, wrap=wrap)

1227

1228

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

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

1231 '''

1232 n = where.__name__

1233 if LatLon is None:

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

1235 else:

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

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

1238 return r

1239

1240

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

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

1243

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

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

1246 height (C{meter}).

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

1248 (C{bool}).

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

1250 or C{None}.

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

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

1253

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

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

1256 is C{None}.

1257

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

1259

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

1261 '''

1262 def _N_vs(ps, w):

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

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

1265 yield p._N_vector

1266

1267 m = _MODS.nvectorBase

1268 # geographic, vectorial mean

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

1270 lat, lon, h = n.latlonheight

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

1272

1273

1274@deprecated_function

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

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

1277

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

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

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

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

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

1283 '''

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

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

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

1287 return ll, d

1288

1289

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

1291 limit=9, **LatLon_and_kwds):

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

1293

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

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

1296

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

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

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

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

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

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

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

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

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

1306 spherical earth radius L{R_KM}).

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

1308 or C{None}.

1309 @kwarg options: Optional keyword arguments for function

1310 L{pygeodesy.equirectangular_}.

1311

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

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

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

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

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

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

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

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

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

1321

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

1323 see function L{pygeodesy.equirectangular_}.

1324

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

1326

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

1328

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

1330

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

1332 '''

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

1334 adjust=adjust, limit=limit)

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

1336 h = t.height

1337 n = nearestOn3.__name__

1338

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

1340 LL = _xkwds_pop(kwds, LatLon=LatLon)

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

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

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

1344

1345

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

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

1348 (with great circle arcs joining the points).

1349

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

1351 or L{BooleanGH}).

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

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

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

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

1356

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

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

1359

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

1361

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

1363

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

1365 C{B{points}} a composite.

1366

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

1368

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

1370 and L{ellipsoidalKarney.perimeterOf}.

1371 '''

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

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

1374 a1, b1 = Ps[0].philam

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

1376 a2, b2 = p.philam

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

1378 yield vincentys_(a2, a1, db)

1379 a1, b1 = a2, b2

1380

1381 if _MODS.booleans.isBoolean(points):

1382 if not closed:

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

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

1385 else:

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

1387 return _r2m(r, radius)

1388

1389

1390 def _r2m(r, radius):

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

1392 '''

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

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

1395 return r

1396

1397

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

1399 excess=excessAbc_,

1400 wrap=False):

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

1402

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

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

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

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

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

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

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

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

1411 or C{None}.

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

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

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

1415 longitudes (C{bool}).

1416

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

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

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

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

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

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

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

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

1425 '''

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

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

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

1429 excess=excess, wrap=wrap)

1430 return _t7Tuple(t, radius)

1431

1432

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

1434 wrap=False):

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

1436 excess} of a (spherical) triangle.

1437

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

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

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

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

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

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

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

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

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

1447 longitudinal deltas (C{bool}).

1448

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

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

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

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

1453 '''

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

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

1456 a = vincentys_(phiC, phiB, d)

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

1458

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

1460 s = sb * sc

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

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

1463

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

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

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

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

1468

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

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

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

1472

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

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

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

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

1477

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

1479

1480

1481 def _t7Tuple(t, radius):

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

1483 '''

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

1485 r = radius if isscalar(radius) else \

1486 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1489 B, (r * t.b),

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

1491 return t

1492

1493

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

1495 areaOf, # functions

1496 intersecant2, intersection, intersections2, ispolar,

1497 isPoleEnclosedBy, # DEPRECATED, use ispolar

1498 meanOf,

1499 nearestOn2, nearestOn3,

1500 perimeterOf,

1501 sumOf, # XXX == vector3d.sumOf

1502 triangle7, triangle8_)

1503

1504# **) MIT License

1505#

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

1507#

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

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

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

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

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

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

1514#

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

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

1517#

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

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

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

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

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

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

1524# OTHER DEALINGS IN THE SOFTWARE.