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

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 v1 = v2

544

545 t = gc.angleTo(n0) > PI_2

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

547 return False # outside

548

549 # check for convex polygon: angle between

550 # gc vectors, signed by direction of n0

551 # (otherwise the test above is not reliable)

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

553 if s != s0:

554 if s0 is None:

555 s0 = s

556 else:

557 t = _Fmt.SQUARE(points=i)

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

559 gc1 = gc

560

561 return True # inside

562

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

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

565

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

567 @kwarg height: Optional height for midpoint, overriding

568 the mean height (C{meter}).

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

570 may be negative or greater than 1.0.

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

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

573

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

575

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

577

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

579

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

581

582 @example:

583

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

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

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

587 '''

588 if fraction is _0_5:

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

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

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

592

593 x = ca2 * cdb + ca1

594 y = ca2 * sdb

595

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

597 b += atan2(y, x)

598

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

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

601 else:

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

603 return r

604

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

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

607

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

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

610

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

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

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

614 @kwarg wrap_adjust_limit: Optional keyword arguments for functions

615 L{sphericalTrigonometry.nearestOn3} and

616 L{pygeodesy.equirectangular_},

617

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

619

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

621 see function L{pygeodesy.equirectangular_}.

622

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

624 is not (yet) supported.

625

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

627

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

629

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

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

632 '''

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

634 w = _xkwds_pop(wrap_adjust_limit, within=True)

635 if not w:

636 notImplemented(self, within=w)

637

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

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

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

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

642# return point1

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

644# if isnear0(d):

645# return point1 # or point2

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

647# return point2

648# f = a / d

649# if within:

650# if f > EPS1:

651# return point2

652# elif f < EPS:

653# return point1

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

655

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

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

658 **wrap_adjust_limit)[0]

659

660 @deprecated_method

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

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

663

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

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

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

667 units of B{C{radius}}.

668 '''

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

670 return r.closest, r.distance

671

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

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

674

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

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

677

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

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

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

681 @kwarg wrap_adjust_limit: Optional keyword arguments for function

682 L{sphericalTrigonometry.nearestOn3} and

683 L{pygeodesy.equirectangular_},

684

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

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

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

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

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

690 function L{pygeodesy.compassAngle}.

691

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

693 see function L{pygeodesy.equirectangular_}.

694

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

696

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

698

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

700

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

702 and L{pygeodesy.nearestOn5}.

703 '''

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

705 LatLon=self.classof, **wrap_adjust_limit)

706

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

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

709 coordinates.

710

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

712 and other keyword arguments, ignored

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

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

715 this L{Cartesian} class or specify

716 C{B{Cartesian}=None}.

717

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

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

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

721

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

723 '''

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

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

726

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

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

729

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

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

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

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

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

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

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

737

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

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

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

741

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

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

744 '''

745 B = self.others(otherB=otherB)

746 C = self.others(otherC=otherC)

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

748

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

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

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

752

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

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

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

756 of three corresponding circles.

757

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

759 as B{C{radius}}).

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

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

762 B{C{radius}}).

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

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

765 B{C{radius}}).

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

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

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

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

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

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

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

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

774

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

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

777 the corresponding trilaterated points C{minPoint} and

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

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

780

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

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

783

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

785 respectively largest I{radial} overlap found.

786

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

788 nearest respectively farthest intersection margin.

789

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

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

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

793 largest B{C{distance#}}.

794

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

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

797 no intersection or all (near-)concentric for

798 C{B{area}=False}.

799

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

801

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

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

804 or B{C{radius}}.

805 '''

806 return _trilaterate5(self, distance1,

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

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

809 area=area, radius=radius, eps=eps, wrap=wrap)

810

811

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

813

814

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

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

817 (with the pointsjoined by great circle arcs).

818

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

820 or L{BooleanGH}).

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

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

823 or C{None}.

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

825 (C{bool}).

826

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

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

829

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

831

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

833

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

835

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

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

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

839

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

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

842 L{ellipsoidalKarney.areaOf}.

843

844 @example:

845

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

847 >>> areaOf(b) # 8666058750.718977

848

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

850 >>> areaOf(c) # 6.18e9

851 '''

852 if _MODS.booleans.isBoolean(points):

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

854

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

856 p1 = p2 = Ps[0]

857 a1, b1 = p1.philam

858 ta1, z1 = tan_2(a1), None

859

860 A = Fsum() # mean phi

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

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

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

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

865 D = Fsum()

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

867 a2, b2 = p2.philam

868 db, b2 = unrollPI(b1, b2, wrap=wrap and not Ps.looped)

869 ta2 = tan_2(a2)

870 A += a2

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

872 _1_0 + ta1 * ta2)

873 ta1, b1 = ta2, b2

874

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

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

877 if z1 is not None:

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

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

880 p1, z1 = p2, z2

881

882 R = fabs(E * _2_0)

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

884 R = fabs(R - PI2)

885 if radius:

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

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

888 return float(R)

889

890

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

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

893

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

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

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

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

898

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

900 '''

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

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

903 ca *= sr

904

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

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

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

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

909

910

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

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

913 # and similar logic in .ellipsoidalBaseDI._intersect3

914 a1, b1 = s.philam

915

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

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

918 else: # must be a point

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

920 hs.append(end.height)

921 a2, b2 = end.philam

922 if wrap:

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

924

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

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

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

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

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

930 -(b1 + b2) * _0_5)

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

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

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

934 cb12 * cb21 + sb12 * sb21,

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

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

937

938

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

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

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

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

943

944

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

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

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

948

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

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

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

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

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

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

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

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

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

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

959 circle} methods.

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

961 conventionally) or C{None}.

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

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

964

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

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

967 point C{is} this very instance.

968

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

970

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

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

973

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

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

976 '''

977 c = _T00.others(center=center)

978 p = _T00.others(point=point)

979 try:

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

981 height=height, wrap=wrap)

982 except (TypeError, ValueError) as x:

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

984

985

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

987 LatLon=None, **LatLon_kwds):

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

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

990

991 s1, s2 = start1, start2

992 if wrap:

993 s2 = _Wrap.point(s2)

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

995

996 a1, b1 = s1.philam

997 a2, b2 = s2.philam

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

999 r12 = vincentys_(a2, a1, db)

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

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

1002

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

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

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

1006

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

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

1009 raise IntersectionError(_parallel_)

1010 # handle domain error for equivalent longitudes,

1011 # see also functions asin_safe and acos_safe at

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

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

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

1015 if sin(db) > 0:

1016 t12, t21 = t1, PI2 - t2

1017 else:

1018 t12, t21 = PI2 - t1, t2

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

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

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

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

1023 raise IntersectionError(_infinite_)

1024 sx3 = sx1 * sx2

1025# XXX if sx3 < 0:

1026# XXX raise ValueError(_ambiguous_)

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

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

1029

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

1031 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

1035

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

1037 _N_vector_ = _MODS.nvectorBase._N_vector_

1038

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

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

1041 x = x1.cross(x2)

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

1043 raise IntersectionError(_colinear_)

1044 a, b = x.philam

1045 # choose intersection similar to sphericalNvector

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

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

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

1049

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

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

1052 intersection, LatLon, LatLon_kwds)

1053

1054

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

1056 LatLon=LatLon, **LatLon_kwds):

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

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

1059

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

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

1062 the initial bearing at the first start point

1063 (compass C{degrees360}).

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

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

1066 the initial bearing at the second start point

1067 (compass C{degrees360}).

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

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

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

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

1072 @kwarg LatLon: Optional class to return the intersection

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

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

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

1076

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

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

1079 height)}. An alternate intersection point might

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

1081

1082 @raise IntersectionError: Ambiguous or infinite intersection

1083 or colinear, parallel or otherwise

1084 non-intersecting lines.

1085

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

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

1088

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

1090

1091 @example:

1092

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

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

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

1096 '''

1097 s1 = _T00.others(start1=start1)

1098 s2 = _T00.others(start2=start2)

1099 try:

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

1101 LatLon=LatLon, **LatLon_kwds)

1102 except (TypeError, ValueError) as x:

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

1104

1105

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

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

1108 LatLon=LatLon, **LatLon_kwds):

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

1110 by a center point and a radius.

1111

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

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

1114 see B{C{radius}}).

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

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

1117 see B{C{radius}}).

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

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

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

1121 B{C{radius}}).

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

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

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

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

1126 (C{bool}).

1127 @kwarg LatLon: Optional class to return the intersection

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

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

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

1131

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

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

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

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

1136

1137 @raise IntersectionError: Concentric, antipodal, invalid or

1138 non-intersecting circles.

1139

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

1141

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

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

1144

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

1146

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

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

1149 '''

1150 c1 = _T00.others(center1=center1)

1151 c2 = _T00.others(center2=center2)

1152 try:

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

1154 height=height, wrap=wrap,

1155 LatLon=LatLon, **LatLon_kwds)

1156 except (TypeError, ValueError) as x:

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

1158 center2=center2, rad2=rad2, wrap=wrap)

1159

1160

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

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

1163 LatLon=LatLon, **LatLon_kwds):

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

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

1166

1167 def _dest3(bearing, h):

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

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

1170 intersections2, LatLon, LatLon_kwds)

1171

1172 a1, b1 = c1.philam

1173 a2, b2 = c2.philam

1174 if wrap:

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

1176

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

1178 if f: # swapped radii, swap centers

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

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

1181

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

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

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

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

1186

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

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

1189 if r < _0_0:

1190 raise _ValueError(eps=r)

1191

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

1193 if x > max(r, EPS):

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

1195 x = sd * sr1

1196 if isnear0(x):

1197 raise _ValueError(_invalid_)

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

1199

1200 elif x < r: # PYCHOK no cover

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

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

1203

1204 if height is None: # "radical height"

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

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

1207 else:

1208 h = Height(height)

1209

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

1211 if x < EPS4: # externally ...

1212 r = _dest3(b, h)

1213 elif x > _PI_EPS4: # internally ...

1214 r = _dest3(b + PI, h)

1215 else:

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

1217 return r, r # ... abutting circles

1218

1219

1220@deprecated_function

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

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

1223 '''

1224 return ispolar(points, wrap=wrap)

1225

1226

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

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

1229 '''

1230 n = func.__name__

1231 if LatLon is None:

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

1233 else:

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

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

1236 return r

1237

1238

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

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

1241

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

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

1244 height (C{meter}).

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

1246 (C{bool}).

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

1248 or C{None}.

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

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

1251

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

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

1254 is C{None}.

1255

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

1257

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

1259 '''

1260 def _N_vs(ps, w):

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

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

1263 yield p._N_vector

1264

1265 m = _MODS.nvectorBase

1266 # geographic, vectorial mean

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

1268 lat, lon, h = n.latlonheight

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

1270

1271

1272@deprecated_function

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

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

1275

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

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

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

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

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

1281 '''

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

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

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

1285 return ll, d

1286

1287

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

1289 limit=9, **LatLon_and_kwds):

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

1291

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

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

1294

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

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

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

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

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

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

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

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

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

1304 spherical earth radius L{R_KM}).

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

1306 or C{None}.

1307 @kwarg options: Optional keyword arguments for function

1308 L{pygeodesy.equirectangular_}.

1309

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

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

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

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

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

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

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

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

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

1319

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

1321 see function L{pygeodesy.equirectangular_}.

1322

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

1324

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

1326

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

1328

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

1330 '''

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

1332 adjust=adjust, limit=limit)

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

1334 h = t.height

1335 n = nearestOn3.__name__

1336

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

1338 LL = _xkwds_pop(kwds, LatLon=LatLon)

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

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

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

1342

1343

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

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

1346 (with great circle arcs joining the points).

1347

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

1349 or L{BooleanGH}).

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

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

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

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

1354

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

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

1357

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

1359

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

1361

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

1363 C{B{points}} a composite.

1364

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

1366

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

1368 and L{ellipsoidalKarney.perimeterOf}.

1369 '''

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

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

1372 a1, b1 = Ps[0].philam

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

1374 a2, b2 = p.philam

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

1376 yield vincentys_(a2, a1, db)

1377 a1, b1 = a2, b2

1378

1379 if _MODS.booleans.isBoolean(points):

1380 if not closed:

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

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

1383 else:

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

1385 return _r2m(r, radius)

1386

1387

1388 def _r2m(r, radius):

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

1390 '''

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

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

1393 return r

1394

1395

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

1397 excess=excessAbc_,

1398 wrap=False):

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

1400

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

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

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

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

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

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

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

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

1409 or C{None}.

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

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

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

1413 longitudes (C{bool}).

1414

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

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

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

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

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

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

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

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

1423 '''

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

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

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

1427 excess=excess, wrap=wrap)

1428 return _t7Tuple(t, radius)

1429

1430

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

1432 wrap=False):

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

1434 excess} of a (spherical) triangle.

1435

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

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

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

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

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

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

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

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

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

1445 longitudinal deltas (C{bool}).

1446

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

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

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

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

1451 '''

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

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

1454 a = vincentys_(phiC, phiB, d)

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

1456

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

1458 s = sb * sc

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

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

1461

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

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

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

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

1466

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

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

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

1470

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

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

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

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

1475

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

1477

1478

1479 def _t7Tuple(t, radius):

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

1481 '''

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

1483 r = radius if isscalar(radius) else \

1484 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1487 B, (r * t.b),

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

1489 return t

1490

1491

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

1493 areaOf, # functions

1494 intersecant2, intersection, intersections2, ispolar,

1495 isPoleEnclosedBy, # DEPRECATED, use ispolar

1496 meanOf,

1497 nearestOn2, nearestOn3,

1498 perimeterOf,

1499 sumOf, # XXX == vector3d.sumOf

1500 triangle7, triangle8_)

1501

1502# **) MIT License

1503#

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

1505#

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

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

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

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

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

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

1512#

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

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

1515#

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

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

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

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

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

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

1522# OTHER DEALINGS IN THE SOFTWARE.