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.

946

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

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

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

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

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

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

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

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

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

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

957 circle} methods.

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

959 conventionally) or C{None}.

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

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

962

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

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

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

966 or I{normalized}.

967

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

969

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

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

972

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

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

975 '''

976 c = _T00.others(center=center)

977 p = _T00.others(point=point)

978 try:

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

980 height=height, wrap=wrap)

981 except (TypeError, ValueError) as x:

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

983

984

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

986 LatLon=None, **LatLon_kwds):

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

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

989

990 s1, s2 = start1, start2

991 if wrap:

992 s2 = _Wrap.point(s2)

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

994

995 a1, b1 = s1.philam

996 a2, b2 = s2.philam

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

998 r12 = vincentys_(a2, a1, db)

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

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

1001

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

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

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

1005

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

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

1008 raise IntersectionError(_parallel_)

1009 # handle domain error for equivalent longitudes,

1010 # see also functions asin_safe and acos_safe at

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

1012 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1)

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

1014 if sin(db) > 0:

1015 t21 = PI2 - t21

1016 else:

1017 t12 = PI2 - t12

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

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

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

1021 raise IntersectionError(_infinite_)

1022 sx3 = sx1 * sx2

1023# XXX if sx3 < 0:

1024# XXX raise ValueError(_ambiguous_)

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

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

1027

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

1029 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

1033

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

1035 _N_vector_ = _MODS.nvectorBase._N_vector_

1036

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

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

1039 x = x1.cross(x2)

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

1041 raise IntersectionError(_colinear_)

1042 a, b = x.philam

1043 # choose intersection similar to sphericalNvector

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

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

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

1047

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

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

1050 intersection, LatLon, LatLon_kwds)

1051

1052

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

1054 LatLon=LatLon, **LatLon_kwds):

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

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

1057

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

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

1060 the initial bearing at the first start point

1061 (compass C{degrees360}).

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

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

1064 the initial bearing at the second start point

1065 (compass C{degrees360}).

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

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

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

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

1070 @kwarg LatLon: Optional class to return the intersection

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

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

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

1074

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

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

1077 height)}. An alternate intersection point might

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

1079

1080 @raise IntersectionError: Ambiguous or infinite intersection

1081 or colinear, parallel or otherwise

1082 non-intersecting lines.

1083

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

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

1086

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

1088

1089 @example:

1090

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

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

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

1094 '''

1095 s1 = _T00.others(start1=start1)

1096 s2 = _T00.others(start2=start2)

1097 try:

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

1099 LatLon=LatLon, **LatLon_kwds)

1100 except (TypeError, ValueError) as x:

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

1102

1103

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

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

1106 LatLon=LatLon, **LatLon_kwds):

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

1108 by a center point and a radius.

1109

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

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

1112 see B{C{radius}}).

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

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

1115 see B{C{radius}}).

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

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

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

1119 B{C{radius}}).

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

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

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

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

1124 (C{bool}).

1125 @kwarg LatLon: Optional class to return the intersection

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

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

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

1129

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

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

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

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

1134

1135 @raise IntersectionError: Concentric, antipodal, invalid or

1136 non-intersecting circles.

1137

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

1139

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

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

1142

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

1144

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

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

1147 '''

1148 c1 = _T00.others(center1=center1)

1149 c2 = _T00.others(center2=center2)

1150 try:

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

1152 height=height, wrap=wrap,

1153 LatLon=LatLon, **LatLon_kwds)

1154 except (TypeError, ValueError) as x:

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

1156 center2=center2, rad2=rad2, wrap=wrap)

1157

1158

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

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

1161 LatLon=LatLon, **LatLon_kwds):

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

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

1164

1165 def _dest3(bearing, h):

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

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

1168 intersections2, LatLon, LatLon_kwds)

1169

1170 a1, b1 = c1.philam

1171 a2, b2 = c2.philam

1172 if wrap:

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

1174

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

1176 if f: # swapped radii, swap centers

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

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

1179

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

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

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

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

1184

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

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

1187 if r < _0_0:

1188 raise _ValueError(eps=r)

1189

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

1191 if x > max(r, EPS):

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

1193 x = sd * sr1

1194 if isnear0(x):

1195 raise _ValueError(_invalid_)

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

1197

1198 elif x < r: # PYCHOK no cover

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

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

1201

1202 if height is None: # "radical height"

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

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

1205 else:

1206 h = Height(height)

1207

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

1209 if x < EPS4: # externally ...

1210 r = _dest3(b, h)

1211 elif x > _PI_EPS4: # internally ...

1212 r = _dest3(b + PI, h)

1213 else:

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

1215 return r, r # ... abutting circles

1216

1217

1218@deprecated_function

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

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

1221 '''

1222 return ispolar(points, wrap=wrap)

1223

1224

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

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

1227 '''

1228 n = where.__name__

1229 if LatLon is None:

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

1231 else:

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

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

1234 return r

1235

1236

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

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

1239

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

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

1242 height (C{meter}).

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

1244 (C{bool}).

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

1246 or C{None}.

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

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

1249

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

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

1252 is C{None}.

1253

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

1255

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

1257 '''

1258 def _N_vs(ps, w):

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

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

1261 yield p._N_vector

1262

1263 m = _MODS.nvectorBase

1264 # geographic, vectorial mean

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

1266 lat, lon, h = n.latlonheight

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

1268

1269

1270@deprecated_function

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

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

1273

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

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

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

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

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

1279 '''

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

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

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

1283 return ll, d

1284

1285

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

1287 limit=9, **LatLon_and_kwds):

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

1289

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

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

1292

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

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

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

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

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

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

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

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

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

1302 spherical earth radius L{R_KM}).

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

1304 or C{None}.

1305 @kwarg options: Optional keyword arguments for function

1306 L{pygeodesy.equirectangular_}.

1307

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

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

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

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

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

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

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

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

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

1317

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

1319 see function L{pygeodesy.equirectangular_}.

1320

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

1322

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

1324

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

1326

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

1328 '''

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

1330 adjust=adjust, limit=limit)

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

1332 h = t.height

1333 n = nearestOn3.__name__

1334

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

1336 LL = _xkwds_pop(kwds, LatLon=LatLon)

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

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

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

1340

1341

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

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

1344 (with great circle arcs joining the points).

1345

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

1347 or L{BooleanGH}).

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

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

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

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

1352

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

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

1355

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

1357

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

1359

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

1361 C{B{points}} a composite.

1362

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

1364

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

1366 and L{ellipsoidalKarney.perimeterOf}.

1367 '''

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

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

1370 a1, b1 = Ps[0].philam

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

1372 a2, b2 = p.philam

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

1374 yield vincentys_(a2, a1, db)

1375 a1, b1 = a2, b2

1376

1377 if _MODS.booleans.isBoolean(points):

1378 if not closed:

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

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

1381 else:

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

1383 return _r2m(r, radius)

1384

1385

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

1387 excess=excessAbc_,

1388 wrap=False):

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

1390

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

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

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

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

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

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

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

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

1399 or C{None}.

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

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

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

1403 longitudes (C{bool}).

1404

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

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

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

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

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

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

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

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

1413 '''

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

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

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

1417 excess=excess, wrap=wrap)

1418 return _t7Tuple(t, radius)

1419

1420

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

1422 wrap=False):

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

1424 excess} of a (spherical) triangle.

1425

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

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

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

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

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

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

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

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

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

1435 longitudinal deltas (C{bool}).

1436

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

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

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

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

1441 '''

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

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

1444 a = vincentys_(phiC, phiB, d)

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

1446

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

1448 s = sb * sc

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

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

1451

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

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

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

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

1456

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

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

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

1460

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

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

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

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

1465

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

1467

1468

1469 def _t7Tuple(t, radius):

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

1471 '''

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

1473 r = radius if isscalar(radius) else \

1474 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1477 B, (r * t.b),

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

1479 return t

1480

1481

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

1483 areaOf, # functions

1484 intersecant2, intersection, intersections2, ispolar,

1485 isPoleEnclosedBy, # DEPRECATED, use ispolar

1486 meanOf,

1487 nearestOn2, nearestOn3,

1488 perimeterOf,

1489 sumOf, # XXX == vector3d.sumOf

1490 triangle7, triangle8_)

1491

1492# **) MIT License

1493#

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

1495#

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

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

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

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

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

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

1502#

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

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

1505#

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

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

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

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

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

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

1512# OTHER DEALINGS IN THE SOFTWARE.