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 _parallel_, _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, _r2m, _trilaterate5

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

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

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

50 m2radians, radiansPI2, sincos2_, tan_2, _unrollon, \

51 unrollPI, _unrollon3, _Wrap, wrap180, wrapPI

52 from pygeodesy.vector3d import sumOf, Vector3d

53

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

55

56__all__ = _ALL_LAZY.sphericalTrigonometry

57__version__ = '23.09.09'

58

59_PI_EPS4 = PI - EPS4

60 if _PI_EPS4 >= PI:

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

62

63

64 class Cartesian(CartesianSphericalBase):

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

66 spherical, geodetic L{LatLon}.

67 '''

68

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

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

71

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

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

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

75

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

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

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

79

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

81 '''

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

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

84

85

86 class LatLon(LatLonSphericalBase):

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

88

89 @example:

90

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

92 '''

93

94 def _ab1_ab2_db5(self, other, wrap):

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

96 '''

97 a1, b1 = self.philam

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

99 if wrap:

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

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

102 else: # unrollPI shortcut

103 db = b2 - b1

104 return a1, b1, a2, b2, db

105

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

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

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

109 end point.

110

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

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

113 from the start point to the point where the perpendicular

114 crosses the line.

115

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

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

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

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

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

121

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

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

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

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

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

127 B{C{start}} point.

128

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

130

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

132

133 @example:

134

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

136

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

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

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

140 '''

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

142 cx = cos(x)

143 return _0_0 if isnear0(cx) else \

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

145

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

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

148 '''

149 s = self.others(start=start)

150 e = self.others(end=end)

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

152

153 r = Radius_(radius)

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

155

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

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

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

159 return r, x, -b

160

161 @deprecated_method

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

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

164 '''

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

166

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

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

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

170

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

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

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

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

175

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

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

178 '''

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

180 sa, ca, sa1, ca1, \

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

182 sa1 *= ca2 * ca

183

184 x = sa1 * sdb

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

186 z = ca1 * sdb * ca2 * sa

187

188 h = hypot(x, y)

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

190 return None # great circle doesn't reach latitude

191

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

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

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

195

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

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

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

199

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

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

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

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

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

205

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

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

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

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

210

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

212

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

214

215 @example:

216

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

218

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

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

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

222 '''

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

224 return _r2m(x, radius)

225

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

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

228 travelled the given distance on the given initial bearing.

229

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

231 B{C{radius}}).

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

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

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

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

236

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

238

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

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

241

242 @example:

243

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

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

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

247 '''

248 a, b = self.philam

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

250

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

252 h = self._heigHt(height)

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

254

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

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

257

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

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

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

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

262

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

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

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

266

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

268

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

270

271 @example:

272

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

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

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

276 '''

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

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

279

280# @Property_RO

281# def Ecef(self):

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

283# '''

284# return _MODS.ecef.EcefKarney

285

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

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

288 on the given initial bearing from this point.

289

290 Direction of vector is such that initial bearing vector

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

292

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

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

295 overriding the default L{Vector3d}.

296 @kwarg Vector_kwds: Optional, additional keyword argunents

297 for B{C{Vector}}.

298

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

300

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

302

303 @example:

304

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

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

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

308 '''

309 a, b = self.philam

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

311

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

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

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

315

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

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

318 to an other point.

319

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

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

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

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

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

325 C{sphericalNvector.LatLon.initialBearingTo}.

326

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

328

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

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

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

332

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

334

335 @example:

336

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

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

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

340 '''

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

342 # XXX behavior like sphericalNvector.LatLon.initialBearingTo

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

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

345

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

347

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

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

350 and an other point.

351

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

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

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

355 @kwarg height: Optional height, overriding the intermediate

356 height (C{meter}).

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

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

359

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

361

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

363

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

365

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

367

368 @example:

369

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

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

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

373 '''

374 p = self

375 f = Scalar(fraction=fraction)

376 if not isnear0(f):

377 p = p.others(other)

378 if wrap:

379 p = _Wrap.point(p)

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

381 a1, b1 = self.philam

382 a2, b2 = p.philam

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

384 r = vincentys_(a2, a1, db)

385 sr = sin(r)

386 if isnon0(sr):

387 sa1, ca1, sa2, ca2, \

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

389

390 t = f * r

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

392 b = sin( t) # / sr superflous

393

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

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

396 z = a * sa1 + b * sa2

397

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

399 b = atan2(y, x)

400

401 else: # PYCHOK no cover

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

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

404

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

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

407 return p

408

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

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

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

412

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

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

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

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

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

418 C{degrees360}).

419 @kwarg height: Optional height for intersection point,

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

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

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

423

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

425 intersection point might be the L{antipode} to

426 the returned result.

427

428 @raise IntersectionError: Ambiguous or infinite intersection

429 or colinear, parallel or otherwise

430 non-intersecting lines.

431

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

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

434

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

436

437 @example:

438

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

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

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

442 '''

443 try:

444 s2 = self.others(other)

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

446 LatLon=self.classof)

447 except (TypeError, ValueError) as x:

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

449 other=other, end2=end2, wrap=wrap)

450

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

452 height=None, wrap=True):

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

454 by a center point and radius.

455

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

457 see B{C{radius}}).

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

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

460 see B{C{radius}}).

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

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

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

464 B{C{radius}}).

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

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

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

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

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

470

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

472 instance. For abutting circles, both intersection

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

474

475 @raise IntersectionError: Concentric, antipodal, invalid or

476 non-intersecting circles.

477

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

479

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

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

482 '''

483 try:

484 c2 = self.others(other)

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

486 height=height, wrap=wrap,

487 LatLon=self.classof)

488 except (TypeError, ValueError) as x:

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

490 other=other, rad2=rad2, wrap=wrap)

491

492 @deprecated_method

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

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

495 return self.isenclosedBy(points)

496

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

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

499

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

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

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

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

504

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

506 composite, C{False} otherwise.

507

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

509

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

511

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

513

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

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

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

517 '''

518 if _MODS.booleans.isBoolean(points):

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

520

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

522 n0 = self._N_vector

523

524 v2 = Ps[0]._N_vector

525 p1 = Ps[1]

526 v1 = p1._N_vector

527 # check whether this point on same side of all

528 # polygon edges (to the left or right depending

529 # on the anti-/clockwise polygon direction)

530 gc1 = v2.cross(v1)

531 t0 = gc1.angleTo(n0) > PI_2

532 s0 = None

533 # get great-circle vector for each edge

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

535 if wrap and not Ps.looped:

536 p2 = _unrollon(p1, p2)

537 p1 = p2

538 v2 = p2._N_vector

539 gc = v1.cross(v2)

540 t = gc.angleTo(n0) > PI_2

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

542 return False # outside

543

544 # check for convex polygon: angle between

545 # gc vectors, signed by direction of n0

546 # (otherwise the test above is not reliable)

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

548 if s != s0:

549 if s0 is None:

550 s0 = s

551 else:

552 t = _Fmt.SQUARE(points=i)

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

554 gc1, v1 = gc, v2

555

556 return True # inside

557

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

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

560

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

562 @kwarg height: Optional height for midpoint, overriding

563 the mean height (C{meter}).

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

565 may be negative or greater than 1.0.

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

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

568

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

570

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

572

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

574

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

576

577 @example:

578

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

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

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

582 '''

583 if fraction is _0_5:

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

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

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

587

588 x = ca2 * cdb + ca1

589 y = ca2 * sdb

590

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

592 b += atan2(y, x)

593

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

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

596 else:

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

598 return r

599

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

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

602

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

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

605

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

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

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

609 @kwarg wrap_adjust_limit: Optional keyword arguments for functions

610 L{sphericalTrigonometry.nearestOn3} and

611 L{pygeodesy.equirectangular_},

612

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

614

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

616 see function L{pygeodesy.equirectangular_}.

617

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

619 is not (yet) supported.

620

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

622

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

624

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

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

627 '''

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

629 w = _xkwds_pop(wrap_adjust_limit, within=True)

630 if not w:

631 notImplemented(self, within=w)

632

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

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

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

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

637# return point1

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

639# if isnear0(d):

640# return point1 # or point2

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

642# return point2

643# f = a / d

644# if within:

645# if f > EPS1:

646# return point2

647# elif f < EPS:

648# return point1

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

650

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

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

653 **wrap_adjust_limit)[0]

654

655 @deprecated_method

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

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

658

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

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

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

662 units of B{C{radius}}.

663 '''

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

665 return r.closest, r.distance

666

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

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

669

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

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

672

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

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

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

676 @kwarg wrap_adjust_limit: Optional keyword arguments for function

677 L{sphericalTrigonometry.nearestOn3} and

678 L{pygeodesy.equirectangular_},

679

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

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

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

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

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

685 function L{pygeodesy.compassAngle}.

686

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

688 see function L{pygeodesy.equirectangular_}.

689

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

691

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

693

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

695

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

697 and L{pygeodesy.nearestOn5}.

698 '''

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

700 LatLon=self.classof, **wrap_adjust_limit)

701

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

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

704 coordinates.

705

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

707 and other keyword arguments, ignored

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

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

710 this L{Cartesian} class or specify

711 C{B{Cartesian}=None}.

712

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

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

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

716

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

718 '''

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

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

721

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

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

724

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

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

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

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

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

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

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

732

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

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

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

736

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

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

739 '''

740 B = self.others(otherB=otherB)

741 C = self.others(otherC=otherC)

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

743

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

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

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

747

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

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

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

751 intersection} of three corresponding circles.

752

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

754 as B{C{radius}}).

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

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

757 B{C{radius}}).

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

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

760 B{C{radius}}).

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

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

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

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

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

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

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

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

769

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

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

772 the corresponding trilaterated points C{minPoint} and

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

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

775

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

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

778

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

780 respectively largest I{radial} overlap found.

781

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

783 nearest respectively farthest intersection margin.

784

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

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

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

788 largest B{C{distance#}}.

789

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

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

792 no intersection or all (near-)concentric for

793 C{B{area}=False}.

794

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

796

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

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

799 or B{C{radius}}.

800 '''

801 return _trilaterate5(self, distance1,

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

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

804 area=area, radius=radius, eps=eps, wrap=wrap)

805

806

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

808

809

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

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

812 (with the pointsjoined by great circle arcs).

813

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

815 or L{BooleanGH}).

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

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

818 or C{None}.

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

820 (C{bool}).

821

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

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

824

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

826

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

828

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

830

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

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

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

834

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

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

837 L{ellipsoidalKarney.areaOf}.

838

839 @example:

840

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

842 >>> areaOf(b) # 8666058750.718977

843

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

845 >>> areaOf(c) # 6.18e9

846 '''

847 if _MODS.booleans.isBoolean(points):

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

849

850 _at2, _t_2 = atan2, tan_2

851 _un, _w180 = unrollPI, wrap180

852

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

854 p1 = p2 = Ps[0]

855 a1, b1 = p1.philam

856 ta1, z1 = _t_2(a1), None

857

858 A = Fsum() # mean phi

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

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

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

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

863 D = Fsum()

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

865 a2, b2 = p2.philam

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

867 A += a2

868 ta2 = _t_2(a2)

869 tdb = _t_2(db, points=i)

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

871 _1_0 + ta1 * ta2)

872 ta1, b1 = ta2, b2

873

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

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

876 if z1 is not None:

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

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

879 p1, z1 = p2, z2

880

881 R = abs(R * _2_0)

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

883 R = abs(R - PI2)

884 if radius:

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

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

887 return float(R)

888

889

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

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

892

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

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

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

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

897

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

899 '''

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

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

902 ca *= sr

903

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

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

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

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

908

909

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

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

912 # and similar logic in .ellipsoidalBaseDI._intersect3

913 a1, b1 = s.philam

914

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

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

917 else: # must be a point

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

919 hs.append(end.height)

920 a2, b2 = end.philam

921 if wrap:

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

923

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

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

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

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

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

929 -(b1 + b2) * _0_5)

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

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

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

933 cb12 * cb21 + sb12 * sb21,

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

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

936

937

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

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

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

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

942

943

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

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

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

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 (L{LatLon}).

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

953 a second 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.