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 _m2radians, CartesianSphericalBase, \

45 _intersecant2, LatLonSphericalBase, \

46 _rads3, _radians2m, _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, atan1d, atan2d, degrees90, degrees180, \

50 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \

51 unrollPI, _unrollon, _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.10.24'

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 _radians2m(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 _radians2m(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 = _m2radians(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 _radians2m(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 = atan1d(z, hypot(x, y))

399 b = atan2d(y, x)

400

401 else: # PYCHOK no cover

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

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

404

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

406 p = self.classof(a, 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 = atan1d(sa1 + sa2, hypot(x, y))

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

593

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

595 r = self.classof(a, 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, other, **radius_exact_height_wrap):

945 '''Compute the intersections of a circle and a (great circle) line given as

946 two points or as a point and bearing.

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 on the (great circle) line (L{LatLon}).

952 @arg other: An other point on the (great circle) line (L{LatLon}) or

953 the bearing at the B{C{point}} (compass C{degrees360}).

954 @kwarg radius_exact_height_wrap: Optional keyword arguments, see

955 method L{LatLon.intersecant2} for further details.

956

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

958 an instance of the B{C{point}} class. Both points are the same

959 instance if the (great circle) line is tangent to the circle.

960

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

962

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

964 B{C{circle}} or B{C{other}} invalid.

965

966 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},

967 B{C{exact}}, B{C{height}} or B{C{napieradius}}.

968 '''

969 c = _T00.others(center=center)

970 p = _T00.others(point=point)

971 try:

972 return _intersecant2(c, circle, p, other, **radius_exact_height_wrap)

973 except (TypeError, ValueError) as x:

974 raise _xError(x, center=center, circle=circle, point=point, other=other,

975 **radius_exact_height_wrap)

976

977

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

979 LatLon=None, **LatLon_kwds):

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

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

982

983 s1, s2 = start1, start2

984 if wrap:

985 s2 = _Wrap.point(s2)

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

987

988 a1, b1 = s1.philam

989 a2, b2 = s2.philam

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

991 r12 = vincentys_(a2, a1, db)

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

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

994

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

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

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

998

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

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

1001 raise IntersectionError(_parallel_)

1002 # handle domain error for equivalent longitudes,

1003 # see also functions asin_safe and acos_safe at

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

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

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

1007 if sin(db) > 0:

1008 t21 = PI2 - t21

1009 else:

1010 t12 = PI2 - t12

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

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

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

1014 raise IntersectionError(_infinite_)

1015 sx3 = sx1 * sx2

1016# XXX if sx3 < 0:

1017# XXX raise ValueError(_ambiguous_)

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

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

1020

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

1022 # like .ellipsoidalBaseDI,_intersect3, if this intersection

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

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

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

1026

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

1028 _N_vector_ = _MODS.nvectorBase._N_vector_

1029

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

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

1032 x = x1.cross(x2)

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

1034 raise IntersectionError(_colinear_)

1035 a, b = x.philam

1036 # choose intersection similar to sphericalNvector

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

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

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

1040

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

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

1043 intersection, LatLon, LatLon_kwds)

1044

1045

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

1047 LatLon=LatLon, **LatLon_kwds):

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

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

1050

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

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

1053 the initial bearing at the first start point

1054 (compass C{degrees360}).

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

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

1057 the initial bearing at the second start point

1058 (compass C{degrees360}).

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

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

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

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

1063 @kwarg LatLon: Optional class to return the intersection

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

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

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

1067

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

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

1070 height)}. An alternate intersection point might

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

1072

1073 @raise IntersectionError: Ambiguous or infinite intersection

1074 or colinear, parallel or otherwise

1075 non-intersecting lines.

1076

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

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

1079

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

1081

1082 @example:

1083

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

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

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

1087 '''

1088 s1 = _T00.others(start1=start1)

1089 s2 = _T00.others(start2=start2)

1090 try:

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

1092 LatLon=LatLon, **LatLon_kwds)

1093 except (TypeError, ValueError) as x:

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

1095

1096

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

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

1099 LatLon=LatLon, **LatLon_kwds):

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

1101 by a center point and a radius.

1102

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

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

1105 see B{C{radius}}).

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

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

1108 see B{C{radius}}).

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

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

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

1112 B{C{radius}}).

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

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

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

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

1117 (C{bool}).

1118 @kwarg LatLon: Optional class to return the intersection

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

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

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

1122

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

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

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

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

1127

1128 @raise IntersectionError: Concentric, antipodal, invalid or

1129 non-intersecting circles.

1130

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

1132

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

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

1135

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

1137

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

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

1140 '''

1141 c1 = _T00.others(center1=center1)

1142 c2 = _T00.others(center2=center2)

1143 try:

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

1145 height=height, wrap=wrap,

1146 LatLon=LatLon, **LatLon_kwds)

1147 except (TypeError, ValueError) as x:

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

1149 center2=center2, rad2=rad2, wrap=wrap)

1150

1151

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

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

1154 LatLon=LatLon, **LatLon_kwds):

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

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

1157

1158 def _dest3(bearing, h):

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

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

1161 intersections2, LatLon, LatLon_kwds)

1162

1163 a1, b1 = c1.philam

1164 a2, b2 = c2.philam

1165 if wrap:

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

1167

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

1169 if f: # swapped radii, swap centers

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

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

1172

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

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

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

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

1177

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

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

1180 if r < _0_0:

1181 raise _ValueError(eps=r)

1182

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

1184 if x > max(r, EPS):

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

1186 x = sd * sr1

1187 if isnear0(x):

1188 raise _ValueError(_invalid_)

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

1190

1191 elif x < r: # PYCHOK no cover

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

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

1194

1195 if height is None: # "radical height"

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

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

1198 else:

1199 h = Height(height)

1200

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

1202 if x < EPS4: # externally ...

1203 r = _dest3(b, h)

1204 elif x > _PI_EPS4: # internally ...

1205 r = _dest3(b + PI, h)

1206 else:

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

1208 return r, r # ... abutting circles

1209

1210

1211@deprecated_function

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

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

1214 '''

1215 return ispolar(points, wrap=wrap)

1216

1217

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

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

1220 '''

1221 n = where.__name__

1222 if LatLon is None:

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

1224 else:

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

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

1227 return r

1228

1229

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

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

1232

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

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

1235 height (C{meter}).

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

1237 (C{bool}).

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

1239 or C{None}.

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

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

1242

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

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

1245 is C{None}.

1246

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

1248

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

1250 '''

1251 def _N_vs(ps, w):

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

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

1254 yield p._N_vector

1255

1256 m = _MODS.nvectorBase

1257 # geographic, vectorial mean

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

1259 lat, lon, h = n.latlonheight

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

1261

1262

1263@deprecated_function

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

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

1266

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

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

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

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

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

1272 '''

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

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

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

1276 return ll, d

1277

1278

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

1280 limit=9, **LatLon_and_kwds):

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

1282

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

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

1285

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

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

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

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

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

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

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

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

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

1295 spherical earth radius L{R_KM}).

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

1297 or C{None}.

1298 @kwarg options: Optional keyword arguments for function

1299 L{pygeodesy.equirectangular_}.

1300

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

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

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

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

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

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

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

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

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

1310

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

1312 see function L{pygeodesy.equirectangular_}.

1313

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

1315

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

1317

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

1319

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

1321 '''

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

1323 adjust=adjust, limit=limit)

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

1325 h = t.height

1326 n = nearestOn3.__name__

1327

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

1329 LL = _xkwds_pop(kwds, LatLon=LatLon)

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

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

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

1333

1334

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

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

1337 (with great circle arcs joining the points).

1338

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

1340 or L{BooleanGH}).

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

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

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

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

1345

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

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

1348

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

1350

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

1352

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

1354 C{B{points}} a composite.

1355

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

1357

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

1359 and L{ellipsoidalKarney.perimeterOf}.

1360 '''

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

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

1363 a1, b1 = Ps[0].philam

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

1365 a2, b2 = p.philam

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

1367 yield vincentys_(a2, a1, db)

1368 a1, b1 = a2, b2

1369

1370 if _MODS.booleans.isBoolean(points):

1371 if not closed:

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

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

1374 else:

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

1376 return _radians2m(r, radius)

1377

1378

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

1380 excess=excessAbc_,

1381 wrap=False):

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

1383

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

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

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

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

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

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

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

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

1392 or C{None}.

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

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

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

1396 longitudes (C{bool}).

1397

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

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

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

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

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

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

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

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

1406 '''

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

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

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

1410 excess=excess, wrap=wrap)

1411 return _t7Tuple(t, radius)

1412

1413

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

1415 wrap=False):

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

1417 excess} of a (spherical) triangle.

1418

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

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

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

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

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

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

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

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

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

1428 longitudinal deltas (C{bool}).

1429

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

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

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

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

1434 '''

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

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

1437 a = vincentys_(phiC, phiB, d)

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

1439

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

1441 s = sb * sc

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

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

1444

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

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

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

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

1449

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

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

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

1453

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

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

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

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

1458

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

1460

1461

1462 def _t7Tuple(t, radius):

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

1464 '''

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

1466 r = radius if isscalar(radius) else \

1467 _ellipsoidal_datum(radius).ellipsoid.Rmean

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

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

1470 B, (r * t.b),

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

1472 return t

1473

1474

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

1476 areaOf, # functions

1477 intersecant2, intersection, intersections2, ispolar,

1478 isPoleEnclosedBy, # DEPRECATED, use ispolar

1479 meanOf,

1480 nearestOn2, nearestOn3,

1481 perimeterOf,

1482 sumOf, # XXX == vector3d.sumOf

1483 triangle7, triangle8_)

1484

1485# **) MIT License

1486#

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

1488#

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

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

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

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

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

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

1495#

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

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

1498#

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

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

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

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

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

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

1505# OTHER DEALINGS IN THE SOFTWARE.