1

2# -*- coding: utf-8 -*-

3

4 u'''Spherical, C{N-vector}-based geodesy.

5

6 N-vector-based classes geodetic (lat-/longitude) L{LatLon}, geocentric

7(ECEF) L{Cartesian} and L{Nvector} and functions L{areaOf}, L{intersection},

8 L{meanOf}, L{nearestOn3}, L{perimeterOf}, L{sumOf}, L{triangulate} and

9 L{trilaterate}, I{all spherical}.

10

11 Pure Python implementation of n-vector-based spherical geodetic (lat-/longitude)

12 methods, transcoded from JavaScript originals by I{(C) Chris Veness 2011-2016},

13 published under the same MIT Licence**. See U{Vector-based geodesy

14<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>} and

15 U{Module latlon-nvector-spherical

16<https://www.Movable-Type.co.UK/scripts/geodesy/docs/module-latlon-nvector-spherical.html>}.

17

18 Tools for working with points and lines on (a spherical model of) the

19 earth’s surface using using n-vectors rather than the more common

20 spherical trigonometry. N-vectors make many calculations much simpler,

21 and easier to follow, compared with the trigonometric equivalents.

22

23 Based on Kenneth Gade’s U{‘Non-singular Horizontal Position Representation’

24<https://www.NavLab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf>},

25 The Journal of Navigation (2010), vol 63, nr 3, pp 395-417.

26

27 Note that the formulations below take x => 0°N,0°E, y => 0°N,90°E and

28 z => 90°N while Gade uses x => 90°N, y => 0°N,90°E, z => 0°N,0°E.

29

30 Also note that on a spherical earth model, an n-vector is equivalent

31 to a normalised version of an (ECEF) cartesian coordinate.

32'''

33# make sure int/int division yields float quosient, see .basics

34 from __future__ import division as _; del _ # PYCHOK semicolon

35

36 from pygeodesy.basics import isscalar, _xinstanceof

37 from pygeodesy.constants import EPS, EPS0, PI, PI2, PI_2, R_M, \

38 _0_0, _0_5, _1_0

39# from pygeodesy.datums import Datums # from .sphericalBase

40 from pygeodesy.errors import PointsError, VectorError, _xError, _xkwds

41 from pygeodesy.fmath import fmean, fsum

42# from pygeodesy.fsums import fsum # from .fmath

43 from pygeodesy.interns import _composite_, _end_, _Nv00_, _other_, \

44 _point_, _pole_

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

46# from pygeodesy.named import notImplemented # from .nvectorBase

47# from pygeodesy.namedTuples import NearestOn3Tuple # from .points

48 from pygeodesy.nvectorBase import LatLonNvectorBase, NorthPole, \

49 notImplemented, NvectorBase, _nsumOf, \

50 _triangulate, _trilaterate

51 from pygeodesy.points import NearestOn3Tuple, ispolar # PYCHOK exported

52 from pygeodesy.props import deprecated_function, deprecated_method, \

53 property_RO

54 from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \

55 _intersecant2, LatLonSphericalBase, \

56 _radians2m, Datums

57 from pygeodesy.units import Bearing, Bearing_, Radius, Scalar

58 from pygeodesy.utily import atan2, degrees360, fabs, sincos2, sincos2_, \

59 sincos2d, _unrollon, _Wrap

60

61# from math import atan2, fabs # from utily

62

63__all__ = _ALL_LAZY.sphericalNvector

64__version__ = '23.10.24'

65

66_lines_ = 'lines'

67

68

69 class Cartesian(CartesianSphericalBase):

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

71 L{Nvector} and n-vector-based, spherical L{LatLon}.

72 '''

73

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

75 '''Convert this cartesian to an C{Nvector}-based geodetic point.

76

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

78 arguments, like C{datum}. Use C{B{LatLon}=...}

79 to override this L{LatLon} class or specify

80 C{B{LatLon}=None}.

81

82 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set

83 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,

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

85

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

87 '''

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

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

90

91 def toNvector(self, **Nvector_and_kwds): # PYCHOK Datums.WGS84

92 '''Convert this cartesian to L{Nvector} components, I{including height}.

93

94 @kwarg Nvector_and_kwds: Optional L{Nvector} and L{Nvector} keyword

95 arguments, like C{datum}. Use C{B{Nvector}=...}

96 to override this L{Nvector} class or specify

97 C{B{Nvector}=None}.

98

99 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}}

100 is set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}

101

102 @raise TypeError: Invalid B{C{Nvector_and_kwds}} argument.

103 '''

104 # ll = CartesianBase.toLatLon(self, LatLon=LatLon,

105 # datum=datum or self.datum)

106 # kwds = _xkwds(kwds, Nvector=Nvector)

107 # return ll.toNvector(**kwds)

108 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector, datum=self.datum)

109 return CartesianSphericalBase.toNvector(self, **kwds)

110

111

112 class LatLon(LatLonNvectorBase, LatLonSphericalBase):

113 '''New n-vector based point on a spherical earth model.

114

115 Tools for working with points, lines and paths on (a spherical

116 model of) the earth's surface using vector-based methods.

117

118 @example:

119

120 >>> from sphericalNvector import LatLon

121 >>> p = LatLon(52.205, 0.119)

122 '''

123 _Nv = None # cached_toNvector L{Nvector})

124

125 def _update(self, updated, *attrs, **setters): # PYCHOK args

126 '''(INTERNAL) Zap cached attributes if updated.

127 '''

128 if updated: # reset caches

129 LatLonNvectorBase._update(self, updated, _Nv=self._Nv) # special case

130 LatLonSphericalBase._update(self, updated, *attrs, **setters)

131

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

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

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

135 end point.

136

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

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

139 from the start point to the point where the perpendicular

140 crosses the line.

141

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

143 @arg end: End point of great circle line (L{LatLon}) or

144 initial bearing from start point (compass

145 C{degrees360}).

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

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

148 the B{C{start}} and B{C{end}} points (C{bool}).

149

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

151 if C{B{radius} is None} else C{meter}, same units

152 as B{C{radius}}), positive if "after" the start

153 toward the end point of the line or negative if

154 "before" the start point.

155

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

157

158 @raise Valuerror: Some points coincide.

159

160 @example:

161

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

163

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

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

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

167 '''

168 p = self.others(start=start)

169 n = self.toNvector()

170

171 gc, _, _ = self._gc3(p, end, _end_, wrap=wrap)

172 a = gc.cross(n).cross(gc) # along-track point gc × p × gc

173 return _radians2m(start.toNvector().angleTo(a, vSign=gc), radius)

174

175 @deprecated_method

176 def bearingTo(self, other, **unused): # PYCHOK no cover

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

178 '''

179 return self.initialBearingTo(other)

180

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

182 '''Compute the (signed) distance from this point to great circle

183 defined by a start and end point.

184

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

186 @arg end: End point of great circle line (L{LatLon}) or

187 initial bearing from start point (compass

188 C{degrees360}).

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

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

191 the B{C{start}} and B{C{end}} points (C{bool}).

192

193 @return: Distance to great circle (C{radians} if C{B{radius}

194 is None} else C{meter}, same units as B{C{radius}}),

195 negative if to the left or positive if to the right

196 of the line .

197

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

199

200 @raise Valuerror: Some points coincide.

201

202 @example:

203

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

205

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

207 >>> d = p.crossTrackDistanceTo(s, 96) # -305.7

208

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

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

211 '''

212 p = self.others(start=start)

213 n = self.toNvector()

214

215 gc, _, _ = self._gc3(p, end, _end_, wrap=wrap)

216 return _radians2m(gc.angleTo(n) - PI_2, radius)

217

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

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

220 travelled the given distance on the given bearing.

221

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

223 as B{C{radius}}).

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

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

226 @kwarg height: Optional height at destination, overriding the

227 default height (C{meter}, same units as B{C{radius}}).

228

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

230

231 @raise Valuerror: Polar coincidence or invalid B{C{distance}},

232 B{C{bearing}}, B{C{radius}} or B{C{height}}.

233

234 @example:

235

236 >>> p = LatLon(51.4778, -0.0015)

237 >>> q = p.destination(7794, 300.7)

238 >>> q.toStr() # 51.513546°N, 000.098345°W

239 '''

240 b = Bearing_(bearing)

241 a = _m2radians(distance, radius, low=None)

242 sa, ca, sb, cb = sincos2_(a, b)

243

244 n = self.toNvector()

245 e = NorthPole.cross(n, raiser=_pole_).unit() # east vector at n

246 x = n.cross(e) # north vector at n

247 d = x.times(cb).plus(e.times(sb)) # direction vector @ n

248 n = n.times(ca).plus(d.times(sa))

249 return n.toLatLon(height=height, LatLon=self.classof) # Nvector(n.x, n.y, n.z).toLatLon(...)

250

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

252 '''Compute the distance from this to an other point.

253

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

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

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

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

258

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

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

261 if B{C{radius}} is C{None}).

262

263 @raise TypeError: Invalid B{C{other}} point.

264

265 @example:

266

267 >>> p = LatLon(52.205, 0.119)

268 >>> q = LatLon(48.857, 2.351);

269 >>> d = p.distanceTo(q) # 404.3 Km

270 '''

271 p = self.others(other)

272 if wrap:

273 p = _unrollon(self, p)

274 n = p.toNvector()

275 r = fabs(self.toNvector().angleTo(n, wrap=wrap))

276 return r if radius is None else (Radius(radius) * r)

277

278# @Property_RO

279# def Ecef(self):

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

281# '''

282# return _ALL_MODS.ecef.EcefKarney

283

284 def _gc3(self, start, end, namend, raiser=_point_, wrap=False):

285 '''(INTERNAL) Return great circle, start and end Nvectors.

286 '''

287 s = start.toNvector()

288 if isscalar(end): # bearing

289 gc = s.greatCircle(end)

290 e = None

291 else: # point

292 p = self.others(end, name=namend)

293 if wrap:

294 p = _unrollon(start, p, wrap=wrap)

295 e = p.toNvector()

296 gc = s.cross(e, raiser=raiser) # XXX .unit()?

297 return gc, s, e

298

299 def greatCircle(self, bearing):

300 '''Compute the vector normal to great circle obtained by

301 heading on the given bearing from this point.

302

303 Direction of vector is such that initial bearing vector

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

305

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

307

308 @return: N-vector representing the great circle (L{Nvector}).

309 '''

310 t = Bearing_(bearing)

311 a, b = self.philam

312

313 sa, ca, sb, cb, st, ct = sincos2_(a, b, t)

314 return Nvector(sb * ct - sa * cb * st,

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

316 ca * st, name=self.name) # XXX .unit()

317

318 def greatCircleTo(self, other, wrap=False):

319 '''Compute the vector normal to great circle obtained by

320 heading from this to an other point or on a given bearing.

321

322 Direction of vector is such that initial bearing vector

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

324

325 @arg other: The other point (L{LatLon}) or the bearing from

326 this point (compass C{degrees360}).

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

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

329

330 @return: N-vector representing the great circle (L{Nvector}).

331

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

333

334 @raise Valuerror: Points coincide.

335

336 @example:

337

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

339 >>> gc = p.greatCircle(96.0)

340 >>> gc.toStr() # (-0.79408, 0.12856, 0.59406)

341

342 >>> q = LatLon(53.1887, 0.1334)

343 >>> g = p.greatCircleTo(q)

344 >>> g.toStr() # (-0.79408, 0.12859, 0.59406)

345 '''

346 gc, _, _ = self._gc3(self, other, _other_, wrap=wrap)

347 return gc.unit()

348

349 def initialBearingTo(self, other, wrap=False, **unused): # raiser=...

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

351 to an other point.

352

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

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

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

356

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

358

359 @raise Crosserror: This point coincides with the B{C{other}}

360 point or the C{NorthPole}, provided

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

362

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

364

365 @example:

366

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

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

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

370 '''

371 n = self.toNvector()

372 p = self.others(other)

373 if wrap:

374 p = _unrollon(self, p, wrap=wrap)

375 p = p.toNvector()

376 # see <https://MathForum.org/library/drmath/view/55417.html>

377# gc1 = self.greatCircleTo(other)

378 gc1 = n.cross(p, raiser=_point_) # .unit()

379# gc2 = self.greatCircleTo(NorthPole)

380 gc2 = n.cross(NorthPole, raiser=_pole_) # .unit()

381 return degrees360(gc1.angleTo(gc2, vSign=n))

382

383 def intermediateChordTo(self, other, fraction, height=None, wrap=False):

384 '''Locate the point projected from the point at given fraction

385 on a straight line (chord) between this and an other point.

386

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

388 @arg fraction: Fraction between both points (float, between

389 0.0 for this and 1.0 for the other point).

390 @kwarg height: Optional height at the intermediate point,

391 overriding the fractional height (C{meter}).

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

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

394

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

396

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

398

399 @example:

400

401 >>> p = LatLon(52.205, 0.119)

402 >>> q = LatLon(48.857, 2.351)

403 >>> i = p.intermediateChordTo(q, 0.25) # 51.3723°N, 000.7072°E

404 '''

405 n = self.toNvector()

406 p = self.others(other)

407 if wrap:

408 p = _unrollon(self, p, wrap=wrap)

409

410 f = Scalar(fraction=fraction)

411 i = p.toNvector().times(f).plus(n.times(1 - f))

412# i = p.toNvector() * f + self.toNvector() * (1 - f))

413

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

415 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)

416

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

418 '''Locate the point at a given fraction between this and an

419 other point.

420

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

422 @arg fraction: Fraction between both points (C{float}, between

423 0.0 for this and 1.0 for the other point).

424 @kwarg height: Optional height at the intermediate point,

425 overriding the fractional height (C{meter}).

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

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

428

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

430

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

432

433 @raise Valuerror: Points coincide or invalid B{C{height}}.

434

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

436

437 @example:

438

439 >>> p = LatLon(52.205, 0.119)

440 >>> q = LatLon(48.857, 2.351)

441 >>> i = p.intermediateTo(q, 0.25) # 51.3721°N, 000.7074°E

442 '''

443 n = self.toNvector()

444 p = self.others(other)

445 if wrap:

446 p = _unrollon(self, p, wrap=wrap)

447 p = p.toNvector()

448 f = Scalar(fraction=fraction)

449

450 x = n.cross(p, raiser=_point_)

451 d = x.unit().cross(n) # unit(n × p) × n

452 # angular distance α, tan(α) = |n × p| / n ⋅ p

453 s, c = sincos2(atan2(x.length, n.dot(p)) * f) # interpolated

454 i = n.times(c).plus(d.times(s)) # n * cosα + d * sinα

455

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

457 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)

458

459 def intersection(self, end1, start2, end2, height=None, wrap=False):

460 '''Locate the intersection point of two lines each defined

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

462

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

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

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

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

467 initial bearing at the second point (compass

468 C{degrees}).

469 @kwarg height: Optional height at the intersection point,

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

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

472 start and end points (C{bool}).

473

474 @return: The intersection point (L{LatLon}).

475

476 @raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}}

477 point is not L{LatLon}.

478

479 @raise ValueError: Intersection is ambiguous or infinite or

480 the lines are parallel, coincident or null.

481

482 @see: Function L{sphericalNvector.intersection} and method

483 L{intersection2}.

484

485 @example:

486

487 >>> s = LatLon(51.8853, 0.2545)

488 >>> e = LatLon(49.0034, 2.5735)

489 >>> i = s.intersection(108.55, e, 32.44) # 50.9076°N, 004.5086°E

490 '''

491 return intersection(self, end1, start2, end2, height=height,

492 wrap=wrap, LatLon=self.classof)

493

494 def intersection2(self, end1, start2, end2, height=None, wrap=False):

495 '''Locate the intersections of two (great circle) lines each defined

496 by two points or by a start point and an (initial) bearing.

497

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

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

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

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

502 initial bearing at the second start point (compass

503 C{degrees360}).

504 @kwarg height: Optional height at the intersection and antipodal

505 point, overriding the mean height (C{meter}).

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

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

508

509 @return: 2-Tuple C{(intersection, antipode)}, each a B{C{LatLon}}.

510

511 @raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}}

512 point is not L{LatLon}.

513

514 @raise ValueError: Intersection is ambiguous or infinite or

515 the lines are parallel, coincident or null.

516

517 @see: Function L{sphericalNvector.intersection2} and method

518 L{intersection}.

519 '''

520 return intersection2(self, end1, start2, end2, height=height,

521 wrap=wrap, LatLon=self.classof)

522

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

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

525

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

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

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

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

530

531 @return: C{True} if this point is inside the polygon or composite,

532 C{False} otherwise.

533

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

535

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

537

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

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

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

541 '''

542 if _MODS.booleans.isBoolean(points):

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

544

545 # sum subtended angles of each edge (using n0, the

546 # normal vector to this point for sign of α)

547 def _subtangles(ps, w):

548 Ps = self.PointsIter(ps, loop=1, wrap=w)

549 n0 = self.toNvector()

550 _m0 = n0.minus

551 p1 = Ps[0]

552 vs1 = _m0(p1.toNvector())

553 for p2 in Ps.iterate(closed=True):

554 if w and not Ps.looped:

555 p2 = _unrollon(p1, p2)

556 p1 = p2

557 vs2 = _m0(p2.toNvector())

558 yield vs1.angleTo(vs2, vSign=n0) # PYCHOK false

559 vs1 = vs2

560

561 # Note, this method uses angle summation test: on a plane,

562 # angles for an enclosed point will sum to 360°, angles for

563 # an exterior point will sum to 0°. On a sphere, enclosed

564 # point angles will sum to less than 360° (due to spherical

565 # excess), exterior point angles will be small but non-zero.

566 s = fsum(_subtangles(points, wrap), floats=True) # normal vector

567 # XXX are winding number optimisations equally applicable to

568 # spherical surface?

569 return fabs(s) > PI

570

571 @deprecated_method

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

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

574 return self.isenclosedBy(points)

575

576 def iswithin(self, point1, point2, wrap=False):

577 '''Check whether this point is between two other points.

578

579 If this point is not on the great circle arc defined by

580 both points, return whether it is within the area bound

581 by perpendiculars to the great circle at each point (in

582 the same hemispere).

583

584 @arg point1: Start point of the arc (L{LatLon}).

585 @arg point2: End point of the arc (L{LatLon}).

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

587 B{C{point1}} and B{C{point2}} (C{bool}).

588

589 @return: C{True} if this point is within the (great circle)

590 arc, C{False} otherwise.

591

592 @raise TypeError: If B{C{point1}} or B{C{point2}} is not

593 L{LatLon}.

594 '''

595 p1 = self.others(point1=point1)

596 p2 = self.others(point2=point2)

597 if wrap:

598 p1 = _Wrap.point(p1)

599 p2 = _unrollon(p1, p2, wrap=wrap)

600 n, n1, n2 = (_.toNvector() for _ in (self, p1, p2))

601

602 # corner case, null arc

603 if n1.isequalTo(n2):

604 return n.isequalTo(n1) or n.isequalTo(n2) # PYCHOK returns

605

606 if n.dot(n1) < 0 or n.dot(n2) < 0: # different hemisphere

607 return False # PYCHOK returns

608

609 # get vectors representing d0=p0->p1 and d2=p2->p1 and the

610 # dot product d0⋅d2 tells us if p0 is on the p2 side of p1 or

611 # on the other side (similarly for d0=p0->p2 and d1=p1->p2

612 # and dot product d0⋅d1 and p0 on the p1 side of p2 or not)

613 return n.minus(n1).dot(n2.minus(n1)) >= 0 and \

614 n.minus(n2).dot(n1.minus(n2)) >= 0

615

616 @deprecated_method

617 def isWithin(self, point1, point2): # PYCHOK no cover

618 '''DEPRECATED, use method C{iswithin}.'''

619 return self.iswithin(point1, point2)

620

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

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

623

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

625 @kwarg height: Optional height at the midpoint, overriding

626 the mean height (C{meter}).

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

628 may be negative or greater than 1.0.

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

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

631

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

633

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

635

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

637

638 @example:

639

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

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

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

643 '''

644 if fraction is _0_5:

645 p = self.others(other)

646 if wrap:

647 p = _unrollon(self, p, wrap=wrap)

648 m = self.toNvector().plus(p.toNvector())

649 h = self._havg(other, f=fraction, h=height)

650 r = m.toLatLon(height=h, LatLon=self.classof)

651 else:

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

653 return r

654

655 def nearestOn(self, point1, point2, height=None, within=True, wrap=False):

656 '''Locate the point on the great circle arc between two points

657 closest to this point.

658

659 @arg point1: Start point of the arc (L{LatLon}).

660 @arg point2: End point of the arc (L{LatLon}).

661 @kwarg height: Optional height, overriding the mean height for

662 the point within the arc (C{meter}), or C{None}

663 to interpolate the height.

664 @kwarg within: If C{True}, return the closest point between both

665 given points, otherwise the closest point

666 elsewhere on the great circle arc (C{bool}).

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

668 B{C{point1}} and B{C{point2}} (C{bool}).

669

670 @return: Closest point on the arc (L{LatLon}).

671

672 @raise NotImplementedError: Keyword argument C{B{wrap}=True}

673 not supported.

674

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

676

677 @example:

678

679 >>> s1 = LatLon(51.0, 1.0)

680 >>> s2 = LatLon(51.0, 2.0)

681

682 >>> s = LatLon(51.0, 1.9)

683 >>> p = s.nearestOn(s1, s2) # 51.0004°N, 001.9000°E

684

685 >>> d = p.distanceTo(s) # 42.71 m

686

687 >>> s = LatLon(51.0, 2.1)

688 >>> p = s.nearestOn(s1, s2) # 51.0000°N, 002.0000°E

689 '''

690 p1 = self.others(point1=point1)

691 p2 = self.others(point2=point2)

692 if wrap:

693 p1 = _Wrap.point(p1)

694 p2 = _unrollon(p1, p2, wrap=wrap)

695 p0 = self

696

697 if p0.iswithin(p1, p2) and not p1.isequalTo(p2, EPS):

698 # closer to arc than to its endpoints,

699 # find the closest point on the arc

700 gc1 = p1.toNvector().cross(p2.toNvector())

701 gc2 = p0.toNvector().cross(gc1)

702 n = gc1.cross(gc2)

703

704 elif within: # for backward compatibility, XXX unwrapped

705 return point1 if (self.distanceTo(point1) <

706 self.distanceTo(point2)) else point2

707

708 else: # handle beyond arc extent by .vector3d.nearestOn

709 n1 = p1.toNvector()

710 n2 = p2.toNvector()

711 n = p0.toNvector().nearestOn(n1, n2, within=False)

712 if n is n1:

713 return p1 # is point1

714 elif n is n2:

715 return p2 # is point2 if not wrap

716

717 p = n.toLatLon(height=height or 0, LatLon=self.classof)

718 if height in (None, False): # interpolate height within extent

719 d = p1.distanceTo(p2)

720 f = (p1.distanceTo(p) / d) if d > EPS0 else _0_5

721 p.height = p1._havg(p2, f=max(_0_0, min(f, _1_0)))

722 return p

723

724 # @deprecated_method

725 def nearestOn2(self, points, **closed_radius_height): # PYCHOK no cover

726 '''DEPRECATED, use method L{sphericalNvector.LatLon.nearestOn3}.

727

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

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

730 to that point from this point ...

731 '''

732 r = self.nearestOn3(points, **closed_radius_height)

733 return r.closest, r.distance

734

735 def nearestOn3(self, points, closed=False, radius=R_M, height=None, wrap=False):

736 '''Locate the point on a path or polygon (with great circle

737 arcs joining consecutive points) closest to this point.

738

739 The closest point is either on within the extent of any great

740 circle arc or the nearest of the arc's end points.

741

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

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

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

745 @kwarg height: Optional height, overriding the mean height

746 for a point within the arc (C{meter}).

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

748 the B{C{points}} (C{bool}).

749

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

751 the C{closest} point (L{LatLon}), the C{distance}

752 between this and the C{closest} point in C{meter},

753 same units as B{C{radius}} or in C{radians} if

754 B{C{radius}} is C{None} and the C{angle} from this

755 to the C{closest} point in compass C{degrees360}.

756

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

758

759 @raise ValueError: No B{C{points}}.

760 '''

761 Ps = self.PointsIter(points, loop=1, wrap=wrap)

762 _r = self.distanceTo

763 _n = self.nearestOn

764

765 c = p1 = Ps[0]

766 r = _r(c, radius=None) # radians

767 for p2 in Ps.iterate(closed=closed):

768 if wrap and not Ps.looped:

769 p2 = _unrollon(p1, p2)

770 p = _n(p1, p2, height=height)

771 d = _r(p, radius=None) # radians

772 if d < r:

773 c, r = p, d

774 p1 = p2

775 d = r if radius is None else (Radius(radius) * r)

776 return NearestOn3Tuple(c, d, degrees360(r))

777

778 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian, datum=None

779 '''Convert this point to C{Nvector}-based cartesian (ECEF) coordinates.

780

781 @kwarg Cartesian_and_kwds: Optional L{Cartesian} and L{Cartesian} keyword

782 arguments, like C{datum}. Use C{B{Cartesian}=...}

783 to override this L{Cartesian} class or specify

784 C{B{Cartesian}=None}.

785

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

787 set to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,

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

789

790 @raise TypeError: Invalid B{C{Cartesian_and_kwds}} argument.

791 '''

792 kwds = _xkwds(Cartesian_and_kwds, Cartesian=Cartesian, datum=self.datum)

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

794

795 def toNvector(self, **Nvector_and_kwds): # PYCHOK signature

796 '''Convert this point to L{Nvector} components, I{including height}.

797

798 @kwarg Nvector_and_kwds: Optional L{Nvector} and L{Nvector} keyword

799 arguments. Specify C{B{Nvector}=...} to

800 override this L{Nvector} class or use

801 C{B{Nvector}=None}.

802

803 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}} is

804 set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}.

805

806 @raise TypeError: Invalid B{C{Nvector_and_kwds}} argument.

807

808 @example:

809

810 >>> p = LatLon(45, 45)

811 >>> n = p.toNvector()

812 >>> n.toStr() # [0.50, 0.50, 0.70710]

813 '''

814 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector)

815 return LatLonNvectorBase.toNvector(self, **kwds)

816

817

818 class Nvector(NvectorBase):

819 '''An n-vector is a position representation using a (unit) vector

820 normal to the earth's surface. Unlike lat-/longitude points,

821 n-vectors have no singularities or discontinuities.

822

823 For many applications, n-vectors are more convenient to work

824 with than other position representations like lat-/longitude,

825 earth-centred earth-fixed (ECEF) vectors, UTM coordinates, etc.

826

827 On a spherical model earth, an n-vector is equivalent to an

828 earth-centred earth-fixed (ECEF) vector.

829

830 Note commonality with L{ellipsoidalNvector.Nvector}.

831 '''

832 _datum = Datums.Sphere # default datum (L{Datum})

833

834 @property_RO

835 def sphericalNvector(self):

836 '''Get this C{Nvector}'s spherical class.

837 '''

838 return type(self)

839

840 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian

841 '''Convert this n-vector to C{Nvector}-based cartesian

842 (ECEF) coordinates.

843

844 @kwarg Cartesian_and_kwds: Optional L{Cartesian} and L{Cartesian} keyword

845 arguments, like C{h}. Use C{B{Cartesian}=...}

846 to override this L{Cartesian} class or specify

847 C{B{Cartesian}=None}.

848

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

850 set to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,

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

852

853 @raise TypeError: Invalid B{C{Cartesian_and_kwds}} argument.

854 '''

855 kwds = _xkwds(Cartesian_and_kwds, h=self.h, Cartesian=Cartesian)

856 return NvectorBase.toCartesian(self, **kwds) # class or .classof

857

858 def toLatLon(self, **LatLon_and_kwds): # PYCHOK height=None, LatLon=LatLon

859 '''Convert this n-vector to an C{Nvector}-based geodetic point.

860

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

862 arguments, like C{height}. Use C{B{LatLon}=...}

863 to override this L{LatLon} class or specify

864 C{B{LatLon}=None}.

865

866 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set

867 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,

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

869

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

871

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

873 '''

874 kwds = _xkwds(LatLon_and_kwds, height=self.h, LatLon=LatLon)

875 return NvectorBase.toLatLon(self, **kwds) # class or .classof

876

877 def greatCircle(self, bearing):

878 '''Compute the n-vector normal to great circle obtained by

879 heading on given compass bearing from this point as its

880 n-vector.

881

882 Direction of vector is such that initial bearing vector

883 b = c × p.

884

885 @arg bearing: Initial compass bearing (C{degrees}).

886

887 @return: N-vector representing great circle (L{Nvector}).

888

889 @raise Valuerror: Polar coincidence.

890

891 @example:

892

893 >>> n = LatLon(53.3206, -1.7297).toNvector()

894 >>> gc = n.greatCircle(96.0) # [-0.794, 0.129, 0.594]

895 '''

896 s, c = sincos2d(Bearing(bearing))

897

898 e = NorthPole.cross(self, raiser=_pole_) # easting

899 n = self.cross(e, raiser=_point_) # northing

900

901 e = e.times(c / e.length)

902 n = n.times(s / n.length)

903 return n.minus(e)

904

905

906_Nvll = LatLon(_0_0, _0_0, name=_Nv00_) # reference instance (L{LatLon})

907

908

909 def areaOf(points, radius=R_M, wrap=False):

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

911 great circle arcs joining consecutive points).

912

913 @arg points: The polygon points or clips (C{LatLon}[],

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

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

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

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

918

919 @return: Polygon area (C{meter} I{squared} , same units as

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

921

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

923

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

925

926 @see: Functions L{pygeodesy.areaOf}, L{sphericalTrigonometry.areaOf}

927 and L{ellipsoidalKarney.areaOf}.

928

929 @example:

930

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

932 >>> areaOf(b) # 8666058750.718977

933 '''

934 def _interangles(ps, w): # like .karney._polygon

935 Ps = _Nvll.PointsIter(ps, loop=2, wrap=w)

936 # use vector to 1st point as plane normal for sign of α

937 n0 = Ps[0].toNvector()

938

939 v2 = Ps[0]._N_vector # XXX v2 == no?

940 p1 = Ps[1]

941 v1 = p1._N_vector

942 gc = v2.cross(v1)

943 for p2 in Ps.iterate(closed=True):

944 if w and not Ps.looped:

945 p2 = _unrollon(p1, p2)

946 p1 = p2

947 v2 = p2._N_vector

948 gc1 = v1.cross(v2)

949 v1 = v2

950 yield gc.angleTo(gc1, vSign=n0)

951 gc = gc1

952

953 if _MODS.booleans.isBoolean(points):

954 r = points._sum2(LatLon, areaOf, radius=None, wrap=wrap)

955 else:

956 # sum interior angles: depending on whether polygon is cw or ccw,

957 # angle between edges is π−α or π+α, where α is angle between

958 # great-circle vectors; so sum α, then take n·π − |Σα| (cannot

959 # use Σ(π−|α|) as concave polygons would fail)

960 s = fsum(_interangles(points, wrap), floats=True)

961 # using Girard’s theorem: A = [Σθᵢ − (n−2)·π]·R²

962 # (PI2 - abs(s) == (n*PI - abs(s)) - (n-2)*PI)

963 r = fabs(PI2 - fabs(s))

964 return r if radius is None else (r * Radius(radius)**2)

965

966

967 def intersecant2(center, circle, point, other, **radius_exact_height_wrap):

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

969 two points or as a point and bearing.

970

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

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

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

974 @arg point: A point on the (great circle) line (L{LatLon}).

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

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

977 @kwarg radius_exact_height_wrap: Optional keyword arguments, see

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

979

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

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

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

983

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

985

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

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

988

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

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

991 '''

992 c = _Nvll.others(center=center)

993 p = _Nvll.others(point=point)

994 try:

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

996 except (TypeError, ValueError) as x:

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

998 **radius_exact_height_wrap)

999

1000

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

1002 LatLon=LatLon, **LatLon_kwds):

1003 '''Locate the intersections of two (great circle) lines each defined

1004 by two points or by a start point and an (initial) bearing.

1005

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

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

1008 bearing at the first start point (compass C{degrees360}).

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

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

1011 bearing at the second start point (compass C{degrees360}).

1012 @kwarg height: Optional height at the intersection point,

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

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

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

1016 @kwarg LatLon: Optional class to return the intersection point

1017 (L{LatLon}).

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

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

1020

1021 @return: The intersection point (B{C{LatLon}}) or if C{B{LatLon}

1022 is None}, a cartesian L{Ecef9Tuple}C{(x, y, z, lat, lon,

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

1024

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

1026

1027 @raise ValueError: Intersection is ambiguous or infinite or

1028 the lines are parallel, coincident or null.

1029

1030 @see: Function L{sphericalNvector.intersection2}.

1031

1032 @example:

1033

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

1035 >>> q = LatLon(49.0034, 2.5735)

1036 >>> i = intersection(p, 108.55, q, 32.44) # 50.9076°N, 004.5086°E

1037 '''

1038 i, _, h = _intersect3(start1, end1, start2, end2, height, wrap)

1039 kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)

1040 return i.toLatLon(**kwds)

1041

1042

1043 def intersection2(start1, end1, start2, end2, height=None, wrap=False,

1044 LatLon=LatLon, **LatLon_kwds):

1045 '''Locate the intersections of two (great circle) lines each defined

1046 by two points or by a start point and an (initial) bearing.

1047

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

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

1050 initial bearing at the first start point

1051 (compass C{degrees360}).

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

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

1054 initial bearing at the second start point

1055 (compass C{degrees360}).

1056 @kwarg height: Optional height at the intersection and antipodal

1057 point, overriding the mean height (C{meter}).

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

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

1060 @kwarg LatLon: Optional class to return the intersection and

1061 antipodal points (L{LatLon}).

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

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

1064

1065 @return: 2-Tuple C{(intersection, antipode)}, each a (B{C{LatLon}})

1066 or if C{B{LatLon} is None}, a cartesian L{Ecef9Tuple}C{(x,

1067 y, z, lat, lon, height, C, M, datum)} with C{C} and C{M}

1068 if available.

1069

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

1071

1072 @raise ValueError: Intersection is ambiguous or infinite or

1073 the lines are parallel, coincident or null.

1074

1075 @see: Function L{sphericalNvector.intersection}.

1076 '''

1077 i, a, h = _intersect3(start1, end1, start2, end2, height, wrap)

1078 kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)

1079 return i.toLatLon(**kwds), a.toLatLon(**kwds)

1080

1081

1082 def _intersect3(start1, end1, start2, end2, height, wrap):

1083 '''(INTERNAL) Return the intersection and antipodal points for

1084 functions C{intersection} and C{intersection2}.

1085 '''

1086 p1 = _Nvll.others(start1=start1)

1087 p2 = _Nvll.others(start2=start2)

1088 if wrap:

1089 p2 = _unrollon(p1, p2, wrap=wrap)

1090 # If gc1 and gc2 are great circles through start and end points

1091 # (or defined by start point and bearing), then the candidate

1092 # intersections are simply gc1 × gc2 and gc2 × gc1. Most of the

1093 # work is deciding the correct intersection point to select! If

1094 # bearing is given, that determines the intersection, but if both

1095 # lines are defined by start/end points, take closer intersection.

1096 gc1, s1, e1 = _Nvll._gc3(p1, end1, 'end1', wrap=wrap)

1097 gc2, s2, e2 = _Nvll._gc3(p2, end2, 'end2', wrap=wrap)

1098

1099 hs = start1.height, start2.height

1100 # there are two (antipodal) candidate intersection

1101 # points ... we have to choose the one to return

1102 i1 = gc1.cross(gc2, raiser=_lines_)

1103 i2 = gc2.cross(gc1, raiser=_lines_)

1104

1105 # selection of intersection point depends on how

1106 # lines are defined (by bearings or endpoints)

1107 if e1 and e2: # endpoint+endpoint

1108 d = sumOf((s1, s2, e1, e2)).dot(i1)

1109 hs += end1.height, end2.height

1110 elif e1 and not e2: # endpoint+bearing

1111 # gc2 x v2 . i1 +ve means v2 bearing points to i1

1112 d = gc2.cross(s2).dot(i1)

1113 hs += end1.height,

1114 elif e2 and not e1: # bearing+endpoint

1115 # gc1 x v1 . i1 +ve means v1 bearing points to i1

1116 d = gc1.cross(s1).dot(i1)

1117 hs += end2.height,

1118 else: # bearing+bearing

1119 # if gc x v . i1 is +ve, initial bearing is

1120 # towards i1, otherwise towards antipodal i2

1121 d1 = gc1.cross(s1).dot(i1) # +ve means p1 bearing points to i1

1122 d2 = gc2.cross(s2).dot(i1) # +ve means p2 bearing points to i1

1123 if d1 > 0 and d2 > 0:

1124 d = 1 # both point to i1

1125 elif d1 < 0 and d2 < 0:

1126 d = -1 # both point to i2

1127 else: # d1, d2 opposite signs

1128 # intersection is at further-away intersection point,

1129 # take opposite intersection from mid- point of v1

1130 # and v2 [is this always true?] XXX changed to always

1131 # get intersection p1 bearing points to, aka being

1132 # located "after" p1 along the bearing at p1, like

1133 # function .sphericalTrigonometry._intersect and

1134 # .ellipsoidalBaseDI._intersect3

1135 d = d1 # neg(s1.plus(s2).dot(i1))

1136

1137 h = fmean(hs) if height is None else height

1138 return (i1, i2, h) if d > 0 else (i2, i1, h)

1139

1140

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

1142 '''Compute the I{geographic} mean of the supplied points.

1143

1144 @arg points: Array of points to be averaged (L{LatLon}[]).

1145 @kwarg height: Optional height, overriding the mean height

1146 (C{meter}).

1147 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{points}} (C{bool}).

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

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

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

1151

1152 @return: Point at geographic mean and mean height (B{C{LatLon}}).

1153

1154 @raise PointsError: Insufficient number of B{C{points}} or

1155 some B{C{points}} are not C{LatLon}.

1156 '''

1157 def _N_vs(ps, w):

1158 Ps = _Nvll.PointsIter(ps, wrap=w)

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

1160 yield p._N_vector

1161

1162 try:

1163 # geographic mean

1164 n = _nsumOf(_N_vs(points, wrap), height, Nvector, {})

1165 except (TypeError, ValueError) as x:

1166 raise PointsError(points=points, wrap=wrap, LatLon=LatLon, cause=x)

1167 return n.toLatLon(**_xkwds(LatLon_kwds, LatLon=LatLon, height=n.h,

1168 name=meanOf.__name__))

1169

1170

1171@deprecated_function

1172 def nearestOn2(point, points, **closed_radius_height): # PYCHOK no cover

1173 '''DEPRECATED, use method L{sphericalNvector.nearestOn3}.

1174

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

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

1177 between the C{closest} and the given B{C{point}} ...

1178 '''

1179 r = nearestOn3(point, points, **closed_radius_height)

1180 return r.closest, r.distance

1181

1182

1183 def nearestOn3(point, points, closed=False, radius=R_M, height=None, wrap=False):

1184 '''Locate the point on a polygon (with great circle arcs joining

1185 consecutive points) closest to an other point.

1186

1187 If the given point is between the end points of a great circle

1188 arc, the closest point is on that arc. Otherwise, the closest

1189 point is the nearest of the arc's end points.

1190

1191 @arg point: The other, reference point (L{LatLon}).

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

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

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

1195 @kwarg height: Optional height, overriding the mean height for

1196 a point within the (great circle) arc (C{meter}).

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

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

1199

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

1201 the C{closest} point (L{LatLon}) on the polygon, the

1202 C{distance} and the C{angle} between the C{closest}

1203 and the given B{C{point}}. The C{distance} is in

1204 C{meter}, same units as B{C{radius}} or in C{radians}

1205 if B{C{radius}} is C{None}, the C{angle} is in compass

1206 C{degrees360}.

1207

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

1209

1210 @raise TypeError: Some B{C{points}} or B{C{point}} not C{LatLon}.

1211 '''

1212 _xinstanceof(LatLon, point=point)

1213

1214 return point.nearestOn3(points, closed=closed, radius=radius,

1215 height=height, wrap=wrap)

1216

1217

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

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

1220 great circle arcs joining consecutive points).

1221

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

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

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

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

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

1227

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

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

1230

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

1232

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

1234

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

1236 C{B{points}} a composite.

1237

1238 @see: Functions L{pygeodesy.perimeterOf}, L{ellipsoidalKarney.perimeterOf}

1239 and L{sphericalTrigonometry.perimeterOf}.

1240 '''

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

1242 Ps = _Nvll.PointsIter(ps, loop=1, wrap=w)

1243 p1 = Ps[0]

1244 v1 = p1._N_vector

1245 for p2 in Ps.iterate(closed=c):

1246 if w and not (c and Ps.looped):

1247 p2 = _unrollon(p1, p2)

1248 p1 = p2

1249 v2 = p2._N_vector

1250 yield v1.angleTo(v2)

1251 v1 = v2

1252

1253 if _MODS.booleans.isBoolean(points):

1254 if not closed:

1255 notImplemented(None, closed=closed, points=_composite_)

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

1257 else:

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

1259 return r if radius is None else (Radius(radius) * r)

1260

1261

1262 def sumOf(nvectors, Vector=Nvector, h=None, **Vector_kwds):

1263 '''Return the I{vectorial} sum of two or more n-vectors.

1264

1265 @arg nvectors: Vectors to be added (L{Nvector}[]).

1266 @kwarg Vector: Optional class for the vectorial sum (L{Nvector}).

1267 @kwarg h: Optional height, overriding the mean height (C{meter}).

1268 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword arguments.

1269

1270 @return: Vectorial sum (B{C{Vector}}).

1271

1272 @raise VectorError: No B{C{nvectors}}.

1273 '''

1274 try:

1275 return _nsumOf(nvectors, h, Vector, Vector_kwds)

1276 except (TypeError, ValueError) as x:

1277 raise VectorError(nvectors=nvectors, Vector=Vector, cause=x)

1278

1279

1280 def triangulate(point1, bearing1, point2, bearing2,

1281 height=None, wrap=False,

1282 LatLon=LatLon, **LatLon_kwds):

1283 '''Locate a point given two known points and the initial bearings

1284 from those points.

1285

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

1287 @arg bearing1: Bearing at the first point (compass C{degrees360}).

1288 @arg point2: Second reference point (L{LatLon}).

1289 @arg bearing2: Bearing at the second point (compass C{degrees360}).

1290 @kwarg height: Optional height at the triangulated point, overriding

1291 the mean height (C{meter}).

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

1293 (C{bool}).

1294 @kwarg LatLon: Optional class to return the triangulated point (L{LatLon}).

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

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

1297

1298 @return: Triangulated point (B{C{LatLon}}).

1299

1300 @raise TypeError: If B{C{point1}} or B{C{point2}} is not L{LatLon}.

1301

1302 @raise Valuerror: Points coincide.

1303

1304 @example:

1305

1306 >>> p = LatLon("47°18.228'N","002°34.326'W") # Basse Castouillet

1307 >>> q = LatLon("47°18.664'N","002°31.717'W") # Basse Hergo

1308 >>> t = triangulate(p, 7, q, 295) # 47.323667°N, 002.568501°W'

1309 '''

1310 return _triangulate(_Nvll.others(point1=point1), bearing1,

1311 _Nvll.others(point2=point2), bearing2,

1312 height=height, wrap=wrap,

1313 LatLon=LatLon, **LatLon_kwds)

1314

1315

1316 def trilaterate(point1, distance1, point2, distance2, point3, distance3, # PYCHOK args

1317 radius=R_M, height=None, useZ=False, wrap=False,

1318 LatLon=LatLon, **LatLon_kwds):

1319 '''Locate a point at given distances from three other points.

1320

1321 @arg point1: First point (L{LatLon}).

1322 @arg distance1: Distance to the first point (C{meter}, same units

1323 as B{C{radius}}).

1324 @arg point2: Second point (L{LatLon}).

1325 @arg distance2: Distance to the second point (C{meter}, same units

1326 as B{C{radius}}).

1327 @arg point3: Third point (L{LatLon}).

1328 @arg distance3: Distance to the third point (C{meter}, same units

1329 as B{C{radius}}).

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

1331 @kwarg height: Optional height at the trilaterated point, overriding

1332 the IDW height (C{meter}, same units as B{C{radius}}).

1333 @kwarg useZ: Include Z component iff non-NaN, non-zero (C{bool}).

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

1335 and B{C{point3}} (C{bool}).

1336 @kwarg LatLon: Optional class to return the trilaterated point (L{LatLon}).

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

1338 ignored if C{B{LatLon} is None}.

1339

1340 @return: Trilaterated point (B{C{LatLon}}).

1341

1342 @raise IntersectionError: No intersection, trilateration failed.

1343

1344 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or B{C{point3}}.

1345

1346 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}},

1347 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.

1348

1349 @see: U{Trilateration<https://WikiPedia.org/wiki/Trilateration>}.

1350 '''

1351 return _trilaterate(_Nvll.others(point1=point1), distance1,

1352 _Nvll.others(point2=point2), distance2,

1353 _Nvll.others(point3=point3), distance3,

1354 radius=radius, height=height, useZ=useZ,

1355 wrap=wrap, LatLon=LatLon, **LatLon_kwds)

1356

1357

1358__all__ += _ALL_OTHER(Cartesian, LatLon, Nvector, # classes

1359 areaOf, # functions

1360 intersecant2, intersection, intersection2, ispolar,

1361 meanOf,

1362 nearestOn2, nearestOn3,

1363 perimeterOf,

1364 sumOf,

1365 triangulate, trilaterate)

1366

1367# **) MIT License

1368#

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

1370#

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

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

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

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

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

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

1377#

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

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

1380#

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

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

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

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

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

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

1387# OTHER DEALINGS IN THE SOFTWARE.