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

34 _line_, _near_, _not_, _null_, _parallel_, _point_, \

35 _SPACE_, _too_

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

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

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

39 from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, NearestOn3Tuple, \

40 Triangle7Tuple, Triangle8Tuple

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

42 Fmt as _Fmt, notImplemented # XXX shadowed

43 from pygeodesy.props import deprecated_function, deprecated_method

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

45 LatLonSphericalBase, _rads3, _r2m, _trilaterate5

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

47 from pygeodesy.units import Bearing_, Height, Lam_, Phi_, Radius_, Scalar

48 from pygeodesy.utily import acos1, asin1, atan1d, atan2d, degrees90, degrees180, \

49 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \

50 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI

51 from pygeodesy.vector3d import sumOf, Vector3d

52

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

54

55__all__ = _ALL_LAZY.sphericalTrigonometry

56__version__ = '23.09.29'

57

58_PI_EPS4 = PI - EPS4

59 if _PI_EPS4 >= PI:

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

61

62

63 class Cartesian(CartesianSphericalBase):

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

65 spherical, geodetic L{LatLon}.

66 '''

67

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

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

70

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

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

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

74

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

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

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

78

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

80 '''

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

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

83

84

85 class LatLon(LatLonSphericalBase):

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

87

88 @example:

89

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

91 '''

92

93 def _ab1_ab2_db5(self, other, wrap):

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

95 '''

96 a1, b1 = self.philam

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

98 if wrap:

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

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

101 else: # unrollPI shortcut

102 db = b2 - b1

103 return a1, b1, a2, b2, db

104

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

106 '''Compute the (signed) distance from the start to the closest

107 point on the great circle line defined by a start and an

108 end point.

109

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

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

112 from the start point to the point where the perpendicular

113 crosses the line.

114

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

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

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

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

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

120

121 @return: Distance along the great circle line (C{radians}

122 if C{B{radius} is None} or C{meter}, same units

123 as B{C{radius}}), positive if I{after} the

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

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

126 B{C{start}} point.

127

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

129

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

131

132 @example:

133

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

135

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

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

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

139 '''

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

141 cx = cos(x)

142 return _0_0 if isnear0(cx) else \

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

144

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

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

147 '''

148 s = self.others(start=start)

149 e = self.others(end=end)

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

151

152 r = Radius_(radius)

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

154

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

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

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

158 return r, x, -b

159

160 @deprecated_method

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

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

163 '''

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

165

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

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

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

169

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

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

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

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

174

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

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

177 '''

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

179 sa, ca, sa1, ca1, \

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

181 sa1 *= ca2 * ca

182

183 x = sa1 * sdb

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

185 z = ca1 * sdb * ca2 * sa

186

187 h = hypot(x, y)

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

189 return None # great circle doesn't reach latitude

190

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

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

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

194

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

196 '''Compute the (signed) distance from this point to

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

198

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

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

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

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

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

204

205 @return: Distance to the great circle (C{radians} if

206 B{C{radius}} or C{meter}, same units as

207 B{C{radius}}), I{negative} if to the left or

208 I{positive} if to the right of the line.

209

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

211

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

213

214 @example:

215

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

217

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

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

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

221 '''

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

223 return _r2m(x, radius)

224

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

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

227 travelled the given distance on the given initial bearing.

228

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

230 B{C{radius}}).

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

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

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

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

235

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

237

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

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

240

241 @example:

242

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

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

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

246 '''

247 a, b = self.philam

248 r, t = _angular(distance, radius, low=None), Bearing_(bearing)

249

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

251 h = self._heigHt(height)

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

253

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

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

256

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

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

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

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

261

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

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

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

265

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

267

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

269

270 @example:

271

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

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

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

275 '''

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

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

278

279# @Property_RO

280# def Ecef(self):

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

282# '''

283# return _MODS.ecef.EcefKarney

284

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

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

287 on the given initial bearing from this point.

288

289 Direction of vector is such that initial bearing vector

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

291

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

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

294 overriding the default L{Vector3d}.

295 @kwarg Vector_kwds: Optional, additional keyword argunents

296 for B{C{Vector}}.

297

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

299

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

301

302 @example:

303

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

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

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

307 '''

308 a, b = self.philam

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

310

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

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

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

314

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

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

317 to an other point.

318

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

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

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

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

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

324 C{sphericalNvector.LatLon.initialBearingTo}.

325

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

327

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

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

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

331

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

333

334 @example:

335

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

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

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

339 '''

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

341 # XXX behavior like sphericalNvector.LatLon.initialBearingTo

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

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

344

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

346

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

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

349 and an other point.

350

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

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

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

354 @kwarg height: Optional height, overriding the intermediate

355 height (C{meter}).

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

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

358

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

360

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

362

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

364

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

366

367 @example:

368

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

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

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

372 '''

373 p = self

374 f = Scalar(fraction=fraction)

375 if not isnear0(f):

376 p = p.others(other)

377 if wrap:

378 p = _Wrap.point(p)

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

380 a1, b1 = self.philam

381 a2, b2 = p.philam

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

383 r = vincentys_(a2, a1, db)

384 sr = sin(r)

385 if isnon0(sr):

386 sa1, ca1, sa2, ca2, \

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

388

389 t = f * r

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

391 b = sin( t) # / sr superflous

392

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

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

395 z = a * sa1 + b * sa2

396

397 a = atan1d(z, hypot(x, y))

398 b = atan2d(y, x)

399

400 else: # PYCHOK no cover

401 a = degrees90( favg(a1, a2, f=f)) # coincident

402 b = degrees180(favg(b1, b2, f=f))

403

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

405 p = self.classof(a, b, height=h)

406 return p

407

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

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

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

411

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

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

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

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

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

417 C{degrees360}).

418 @kwarg height: Optional height for intersection point,

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

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

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

422

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

424 intersection point might be the L{antipode} to

425 the returned result.

426

427 @raise IntersectionError: Ambiguous or infinite intersection

428 or colinear, parallel or otherwise

429 non-intersecting lines.

430

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

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

433

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

435

436 @example:

437

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

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

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

441 '''

442 try:

443 s2 = self.others(other)

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

445 LatLon=self.classof)

446 except (TypeError, ValueError) as x:

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

448 other=other, end2=end2, wrap=wrap)

449

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

451 height=None, wrap=True):

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

453 by a center point and radius.

454

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

456 see B{C{radius}}).

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

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

459 see B{C{radius}}).

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

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

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

463 B{C{radius}}).

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

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

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

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

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

469

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

471 instance. For abutting circles, both intersection

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

473

474 @raise IntersectionError: Concentric, antipodal, invalid or

475 non-intersecting circles.

476

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

478

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

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

481 '''

482 try:

483 c2 = self.others(other)

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

485 height=height, wrap=wrap,

486 LatLon=self.classof)

487 except (TypeError, ValueError) as x:

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

489 other=other, rad2=rad2, wrap=wrap)

490

491 @deprecated_method

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

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

494 return self.isenclosedBy(points)

495

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

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

498

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

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

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

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

503

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

505 composite, C{False} otherwise.

506

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

508

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

510

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

512

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

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

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

516 '''

517 if _MODS.booleans.isBoolean(points):

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

519

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

521 n0 = self._N_vector

522

523 v2 = Ps[0]._N_vector

524 p1 = Ps[1]

525 v1 = p1._N_vector

526 # check whether this point on same side of all

527 # polygon edges (to the left or right depending

528 # on the anti-/clockwise polygon direction)

529 gc1 = v2.cross(v1)

530 t0 = gc1.angleTo(n0) > PI_2

531 s0 = None

532 # get great-circle vector for each edge

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

534 if wrap and not Ps.looped:

535 p2 = _unrollon(p1, p2)

536 p1 = p2

537 v2 = p2._N_vector

538 gc = v1.cross(v2)

539 t = gc.angleTo(n0) > PI_2

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

541 return False # outside

542

543 # check for convex polygon: angle between

544 # gc vectors, signed by direction of n0

545 # (otherwise the test above is not reliable)

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

547 if s != s0:

548 if s0 is None:

549 s0 = s

550 else:

551 t = _Fmt.SQUARE(points=i)

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

553 gc1, v1 = gc, v2

554

555 return True # inside

556

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

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

559

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

561 @kwarg height: Optional height for midpoint, overriding

562 the mean height (C{meter}).

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

564 may be negative or greater than 1.0.

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

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

567

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

569

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

571

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

573

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

575

576 @example:

577

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

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

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

581 '''

582 if fraction is _0_5:

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

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

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

586

587 x = ca2 * cdb + ca1

588 y = ca2 * sdb

589

590 a = atan1d(sa1 + sa2, hypot(x, y))

591 b = degrees180(b + atan2(y, x))

592

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

594 r = self.classof(a, b, height=h)

595 else:

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

597 return r

598

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

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

601

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

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

604

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

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

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

608 @kwarg wrap_adjust_limit: Optional keyword arguments for functions

609 L{sphericalTrigonometry.nearestOn3} and

610 L{pygeodesy.equirectangular_},

611

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

613

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

615 see function L{pygeodesy.equirectangular_}.

616

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

618 is not (yet) supported.

619

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

621

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

623

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

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

626 '''

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

628 w = _xkwds_pop(wrap_adjust_limit, within=True)

629 if not w:

630 notImplemented(self, within=w)

631

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

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

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

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

636# return point1

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

638# if isnear0(d):

639# return point1 # or point2

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

641# return point2

642# f = a / d

643# if within:

644# if f > EPS1:

645# return point2

646# elif f < EPS:

647# return point1

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

649

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

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

652 **wrap_adjust_limit)[0]

653

654 @deprecated_method

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

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

657

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

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

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

661 units of B{C{radius}}.

662 '''

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

664 return r.closest, r.distance

665

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

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

668

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

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

671

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

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

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

675 @kwarg wrap_adjust_limit: Optional keyword arguments for function

676 L{sphericalTrigonometry.nearestOn3} and

677 L{pygeodesy.equirectangular_},

678

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

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

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

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

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

684 function L{pygeodesy.compassAngle}.

685

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

687 see function L{pygeodesy.equirectangular_}.

688

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

690

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

692

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

694

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

696 and L{pygeodesy.nearestOn5}.

697 '''

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

699 LatLon=self.classof, **wrap_adjust_limit)

700

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

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

703 coordinates.

704

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

706 and other keyword arguments, ignored

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

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

709 this L{Cartesian} class or specify

710 C{B{Cartesian}=None}.

711

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

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

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

715

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

717 '''

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

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

720

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

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

723

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

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

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

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

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

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

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

731

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

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

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

735

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

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

738 '''

739 B = self.others(otherB=otherB)

740 C = self.others(otherC=otherC)

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

742

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

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

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

746

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

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

749 '''Trilaterate three points by I{area overlap} or I{perimeter

750 intersection} of three corresponding circles.

751

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

753 as B{C{radius}}).

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

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

756 B{C{radius}}).

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

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

759 B{C{radius}}).

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

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

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

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

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

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

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

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

768

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

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

771 the corresponding trilaterated points C{minPoint} and

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

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

774

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

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

777

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

779 respectively largest I{radial} overlap found.

780

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

782 nearest respectively farthest intersection margin.

783

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

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

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

787 largest B{C{distance#}}.

788

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

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

791 no intersection or all (near-)concentric for

792 C{B{area}=False}.

793

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

795

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

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

798 or B{C{radius}}.

799 '''

800 return _trilaterate5(self, distance1,

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

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

803 area=area, radius=radius, eps=eps, wrap=wrap)

804

805

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

807

808

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

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

811 (with the pointsjoined by great circle arcs).

812

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

814 or L{BooleanGH}).

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

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

817 or C{None}.

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

819 (C{bool}).

820

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

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

823

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

825

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

827

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

829

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

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

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

833

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

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

836 L{ellipsoidalKarney.areaOf}.

837

838 @example:

839

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

841 >>> areaOf(b) # 8666058750.718977

842

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

844 >>> areaOf(c) # 6.18e9

845 '''

846 if _MODS.booleans.isBoolean(points):

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

848

849 _at2, _t_2 = atan2, tan_2

850 _un, _w180 = unrollPI, wrap180

851

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

853 p1 = p2 = Ps[0]

854 a1, b1 = p1.philam

855 ta1, z1 = _t_2(a1), None

856

857 A = Fsum() # mean phi

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

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

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

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

862 D = Fsum()

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

864 a2, b2 = p2.philam

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

866 A += a2

867 ta2 = _t_2(a2)

868 tdb = _t_2(db, points=i)

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

870 _1_0 + ta1 * ta2)

871 ta1, b1 = ta2, b2

872

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

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

875 if z1 is not None:

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

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

878 p1, z1 = p2, z2

879

880 R = abs(R * _2_0)

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

882 R = abs(R - PI2)

883 if radius:

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

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

886 return float(R)

887

888

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

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

891

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

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

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

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

896

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

898 '''

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

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

901 ca *= sr

902

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

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

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

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

907

908

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

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

911 # and similar logic in .ellipsoidalBaseDI._intersect3

912 a1, b1 = s.philam

913

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

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

916 else: # must be a point

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

918 hs.append(end.height)

919 a2, b2 = end.philam

920 if wrap:

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

922

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

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

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

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

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

928 -(b1 + b2) * _0_5)

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

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

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

932 cb12 * cb21 + sb12 * sb21,

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

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

935

936

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

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

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

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

941

942

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

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

945 '''Compute the intersections of a circle and a line given as a point and

946 bearing or as two points.

947

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

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

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

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

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

953 an other point on the line (L{LatLon}).

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

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

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

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

958 circle} methods.

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

960 conventionally) or C{None}.

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

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

963

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

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

966 points are the same instance, the B{C{point}} or wrapped

967 or I{normalized}.

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 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1)

1014 t21, t23 = acos1((sa1 - sa2 * cr12) / x2), radiansPI2(end2)

1015 if sin(db) > 0:

1016 t21 = PI2 - t21

1017 else:

1018 t12 = PI2 - t12

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

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

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

1022 raise IntersectionError(_infinite_)

1023 sx3 = sx1 * sx2

1024# XXX if sx3 < 0:

1025# XXX raise ValueError(_ambiguous_)

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

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

1028

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

1030 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

1034

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

1036 _N_vector_ = _MODS.nvectorBase._N_vector_

1037

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

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

1040 x = x1.cross(x2)

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

1042 raise IntersectionError(_colinear_)

1043 a, b = x.philam

1044 # choose intersection similar to sphericalNvector

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

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

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

1048

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

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

1051 intersection, LatLon, LatLon_kwds)

1052

1053

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

1055 LatLon=LatLon, **LatLon_kwds):

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

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

1058

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

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

1061 the initial bearing at the first start point

1062 (compass C{degrees360}).

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

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

1065 the initial bearing at the second start point

1066 (compass C{degrees360}).

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

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

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

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

1071 @kwarg LatLon: Optional class to return the intersection

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

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

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

1075

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

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

1078 height)}. An alternate intersection point might

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

1080

1081 @raise IntersectionError: Ambiguous or infinite intersection

1082 or colinear, parallel or otherwise

1083 non-intersecting lines.

1084

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

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

1087

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

1089

1090 @example:

1091

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

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

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

1095 '''

1096 s1 = _T00.others(start1=start1)

1097 s2 = _T00.others(start2=start2)

1098 try:

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

1100 LatLon=LatLon, **LatLon_kwds)

1101 except (TypeError, ValueError) as x:

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

1103

1104

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

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

1107 LatLon=LatLon, **LatLon_kwds):

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

1109 by a center point and a radius.

1110

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

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

1113 see B{C{radius}}).

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

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

1116 see B{C{radius}}).

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

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

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

1120 B{C{radius}}).

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

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

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

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

1125 (C{bool}).

1126 @kwarg LatLon: Optional class to return the intersection

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

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

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

1130

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

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

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

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

1135

1136 @raise IntersectionError: Concentric, antipodal, invalid or

1137 non-intersecting circles.

1138

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

1140

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

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

1143

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

1145

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

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

1148 '''

1149 c1 = _T00.others(center1=center1)

1150 c2 = _T00.others(center2=center2)

1151 try:

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

1153 height=height, wrap=wrap,

1154 LatLon=LatLon, **LatLon_kwds)

1155 except (TypeError, ValueError) as x:

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

1157 center2=center2, rad2=rad2, wrap=wrap)

1158

1159

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

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

1162 LatLon=LatLon, **LatLon_kwds):

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

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

1165

1166 def _dest3(bearing, h):

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

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

1169 intersections2, LatLon, LatLon_kwds)

1170

1171 a1, b1 = c1.philam

1172 a2, b2 = c2.philam

1173 if wrap:

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

1175

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

1177 if f: # swapped radii, swap centers

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

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

1180

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

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

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

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

1185

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

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

1188 if r < _0_0:

1189 raise _ValueError(eps=r)

1190

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

1192 if x > max(r, EPS):

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

1194 x = sd * sr1

1195 if isnear0(x):

1196 raise _ValueError(_invalid_)

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

1198

1199 elif x < r: # PYCHOK no cover

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

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

1202

1203 if height is None: # "radical height"

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

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

1206 else:

1207 h = Height(height)

1208

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

1210 if x < EPS4: # externally ...

1211 r = _dest3(b, h)

1212 elif x > _PI_EPS4: # internally ...

1213 r = _dest3(b + PI, h)

1214 else:

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

1216 return r, r # ... abutting circles

1217

1218

1219@deprecated_function

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

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

1222 '''

1223 return ispolar(points, wrap=wrap)

1224

1225

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

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

1228 '''

1229 n = where.__name__

1230 if LatLon is None:

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

1232 else:

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

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

1235 return r

1236

1237

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

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

1240

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

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

1243 height (C{meter}).

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

1245 (C{bool}).

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

1247 or C{None}.

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

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

1250

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

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

1253 is C{None}.

1254

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

1256

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

1258 '''

1259 def _N_vs(ps, w):

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

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

1262 yield p._N_vector

1263

1264 m = _MODS.nvectorBase

1265 # geographic, vectorial mean

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

1267 lat, lon, h = n.latlonheight

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

1269

1270

1271@deprecated_function

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

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

1274

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

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

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

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

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

1280 '''

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

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

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

1284 return ll, d

1285

1286

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

1288 limit=9, **LatLon_and_kwds):

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

1290

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

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

1293

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

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

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

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

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

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

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

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

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

1303 spherical earth radius L{R_KM}).

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

1305 or C{None}.

1306 @kwarg options: Optional keyword arguments for function

1307 L{pygeodesy.equirectangular_}.

1308

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

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

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

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

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

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

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

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

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

1318

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

1320 see function L{pygeodesy.equirectangular_}.

1321

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

1323

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

1325

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

1327

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

1329 '''

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

1331 adjust=adjust, limit=limit)

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

1333 h = t.height

1334 n = nearestOn3.__name__

1335

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

1337 LL = _xkwds_pop(kwds, LatLon=LatLon)

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

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

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

1341

1342

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

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

1345 (with great circle arcs joining the points).

1346

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

1348 or L{BooleanGH}).

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

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

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

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

1353

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

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

1356

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

1358

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

1360

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

1362 C{B{points}} a composite.

1363

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

1365

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

1367 and L{ellipsoidalKarney.perimeterOf}.

1368 '''

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

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

1371 a1, b1 = Ps[0].philam

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

1373 a2, b2 = p.philam

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

1375 yield vincentys_(a2, a1, db)

1376 a1, b1 = a2, b2

1377

1378 if _MODS.booleans.isBoolean(points):

1379 if not closed:

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

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

1382 else:

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

1384 return _r2m(r, radius)

1385

1386

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

1388 excess=excessAbc_,

1389 wrap=False):

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

1391

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

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

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

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

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

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

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

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

1400 or C{None}.

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

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

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

1404 longitudes (C{bool}).

1405

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

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

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

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

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

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

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

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

1414 '''

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

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

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

1418 excess=excess, wrap=wrap)

1419 return _t7Tuple(t, radius)

1420

1421

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

1423 wrap=False):

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

1425 excess} of a (spherical) triangle.

1426

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

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

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

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

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

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

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

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

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

1436 longitudinal deltas (C{bool}).

1437

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

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

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

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

1442 '''

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

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

1445 a = vincentys_(phiC, phiB, d)

1446 return a, (phiB, lamB, phiC, lamC, phiA, lamA) # rotate A, B, C

1447

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

1449 s = sb * sc

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

1451 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2_'s

1452

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

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

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

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

1457

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

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

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

1461

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

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

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

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

1466

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

1468

1469

1470 def _t7Tuple(t, radius):

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

1472 '''

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

1474 r = radius if isscalar(radius) else \

1475 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1478 B, (r * t.b),

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

1480 return t

1481

1482

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

1484 areaOf, # functions

1485 intersecant2, intersection, intersections2, ispolar,

1486 isPoleEnclosedBy, # DEPRECATED, use ispolar

1487 meanOf,

1488 nearestOn2, nearestOn3,

1489 perimeterOf,

1490 sumOf, # XXX == vector3d.sumOf

1491 triangle7, triangle8_)

1492

1493# **) MIT License

1494#

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

1496#

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

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

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

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

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

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

1503#

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

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

1506#

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

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

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

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

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

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

1513# OTHER DEALINGS IN THE SOFTWARE.