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 from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, \

54 _intersecant2, LatLonSphericalBase, \

55 Datums

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

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

58 sincos2d, _unrollon, _Wrap

59

60# from math import atan2, fabs # from utily

61

62__all__ = _ALL_LAZY.sphericalNvector

63__version__ = '23.05.05'

64

65_lines_ = 'lines'

66

67

68 class Cartesian(CartesianSphericalBase):

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

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

71 '''

72

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

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

75

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

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

78 to override this L{LatLon} class or specify

79 C{B{LatLon}=None}.

80

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

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

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

84

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

86 '''

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

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

89

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

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

92

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

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

95 to override this L{Nvector} class or specify

96 C{B{Nvector}=None}.

97

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

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

100

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

102 '''

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

104 # datum=datum or self.datum)

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

106 # return ll.toNvector(**kwds)

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

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

109

110

111 class LatLon(LatLonNvectorBase, LatLonSphericalBase):

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

113

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

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

116

117 @example:

118

119 >>> from sphericalNvector import LatLon

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

121 '''

122 _Nv = None # cached_toNvector L{Nvector})

123

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

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

126 '''

127 if updated: # reset caches

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

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

130

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

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

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

134 end point.

135

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

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

138 from the start point to the point where the perpendicular

139 crosses the line.

140

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

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

143 initial bearing from start point (compass

144 C{degrees360}).

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

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

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

148

149 @return: Distance along the great circle line (positive if

150 after the start toward the end point of the line

151 or negative if before the start point).

152

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

154

155 @raise Valuerror: Some points coincide.

156

157 @example:

158

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

160

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

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

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

164 '''

165 p = self.others(start=start)

166 n = self.toNvector()

167

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

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

170 return start.toNvector().angleTo(a, vSign=gc) * radius

171

172 @deprecated_method

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

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

175 '''

176 return self.initialBearingTo(other)

177

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

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

180 defined by a start and end point.

181

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

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

184 initial bearing from start point (compass

185 C{degrees360}).

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

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

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

189

190 @return: Distance to great circle (negative if to the

191 left or positive if to the right of the line).

192

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

194

195 @raise Valuerror: Some points coincide.

196

197 @example:

198

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

200

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

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

203

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

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

206 '''

207 p = self.others(start=start)

208 n = self.toNvector()

209

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

211 return (gc.angleTo(n) - PI_2) * radius

212

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

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

215 travelled the given distance on the given bearing.

216

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

218 as B{C{radius}}).

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

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

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

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

223

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

225

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

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

228

229 @example:

230

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

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

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

234 '''

235 b = Bearing_(bearing)

236 a = _angular(distance, radius)

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

238

239 n = self.toNvector()

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

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

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

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

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

245

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

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

248

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

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

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

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

253

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

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

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

257

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

259

260 @example:

261

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

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

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

265 '''

266 p = self.others(other)

267 if wrap:

268 p = _unrollon(self, p)

269 n = p.toNvector()

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

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

272

273# @Property_RO

274# def Ecef(self):

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

276# '''

277# return _ALL_MODS.ecef.EcefKarney

278

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

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

281 '''

282 s = start.toNvector()

283 if isscalar(end): # bearing

284 gc = s.greatCircle(end)

285 e = None

286 else: # point

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

288 if wrap:

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

290 e = p.toNvector()

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

292 return gc, s, e

293

294 def greatCircle(self, bearing):

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

296 heading on the given bearing from this point.

297

298 Direction of vector is such that initial bearing vector

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

300

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

302

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

304 '''

305 t = Bearing_(bearing)

306 a, b = self.philam

307

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

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

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

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

312

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

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

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

316

317 Direction of vector is such that initial bearing vector

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

319

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

321 this point (compass C{degrees360}).

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

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

324

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

326

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

328

329 @raise Valuerror: Points coincide.

330

331 @example:

332

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

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

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

336

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

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

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

340 '''

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

342 return gc.unit()

343

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

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

346 to an other point.

347

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

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

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

351

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

353

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

355 point or the C{NorthPole}, provided

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

357

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

359

360 @example:

361

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

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

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

365 '''

366 n = self.toNvector()

367 p = self.others(other)

368 if wrap:

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

370 p = p.toNvector()

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

372# gc1 = self.greatCircleTo(other)

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

374# gc2 = self.greatCircleTo(NorthPole)

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

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

377

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

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

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

381

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

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

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

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

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

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

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

389

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

391

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

393

394 @example:

395

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

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

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

399 '''

400 n = self.toNvector()

401 p = self.others(other)

402 if wrap:

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

404

405 f = Scalar(fraction=fraction)

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

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

408

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

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

411

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

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

414 other point.

415

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

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

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

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

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

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

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

423

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

425

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

427

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

429

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

431

432 @example:

433

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

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

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

437 '''

438 n = self.toNvector()

439 p = self.others(other)

440 if wrap:

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

442 p = p.toNvector()

443 f = Scalar(fraction=fraction)

444

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

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

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

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

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

450

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

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

453

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

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

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

457

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

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

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

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

462 initial bearing at the second point (compass

463 C{degrees}).

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

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

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

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

468

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

470

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

472 point is not L{LatLon}.

473

474 @raise ValueError: Intersection is ambiguous or infinite or

475 the lines are parallel, coincident or null.

476

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

478 L{intersection2}.

479

480 @example:

481

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

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

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

485 '''

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

487 wrap=wrap, LatLon=self.classof)

488

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

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

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

492

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

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

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

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

497 initial bearing at the second start point (compass

498 C{degrees360}).

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

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

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

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

503

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

505

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

507 point is not L{LatLon}.

508

509 @raise ValueError: Intersection is ambiguous or infinite or

510 the lines are parallel, coincident or null.

511

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

513 L{intersection}.

514 '''

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

516 wrap=wrap, LatLon=self.classof)

517

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

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

520

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

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

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

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

525

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

527 C{False} otherwise.

528

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

530

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

532

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

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

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

536 '''

537 if _MODS.booleans.isBoolean(points):

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

539

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

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

542 def _subtangles(ps, w):

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

544 n0 = self.toNvector()

545 _m0 = n0.minus

546 p1 = Ps[0]

547 vs1 = _m0(p1.toNvector())

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

549 if w and not Ps.looped:

550 p2 = _unrollon(p1, p2)

551 p1 = p2

552 vs2 = _m0(p2.toNvector())

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

554 vs1 = vs2

555

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

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

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

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

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

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

562 # XXX are winding number optimisations equally applicable to

563 # spherical surface?

564 return fabs(s) > PI

565

566 @deprecated_method

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

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

569 return self.isenclosedBy(points)

570

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

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

573

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

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

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

577 the same hemispere).

578

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

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

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

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

583

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

585 arc, C{False} otherwise.

586

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

588 L{LatLon}.

589 '''

590 p1 = self.others(point1=point1)

591 p2 = self.others(point2=point2)

592 if wrap:

593 p1 = _Wrap.point(p1)

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

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

596

597 # corner case, null arc

598 if n1.isequalTo(n2):

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

600

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

602 return False # PYCHOK returns

603

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

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

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

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

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

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

610

611 @deprecated_method

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

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

614 return self.iswithin(point1, point2)

615

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

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

618

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

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

621 the mean height (C{meter}).

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

623 may be negative or greater than 1.0.

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

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

626

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

628

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

630

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

632

633 @example:

634

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

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

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

638 '''

639 if fraction is _0_5:

640 p = self.others(other)

641 if wrap:

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

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

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

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

646 else:

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

648 return r

649

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

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

652 closest to this point.

653

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

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

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

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

658 to interpolate the height.

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

660 given points, otherwise the closest point

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

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

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

664

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

666

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

668 not supported.

669

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

671

672 @example:

673

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

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

676

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

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

679

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

681

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

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

684 '''

685 p1 = self.others(point1=point1)

686 p2 = self.others(point2=point2)

687 if wrap:

688 p1 = _Wrap.point(p1)

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

690 p0 = self

691

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

693 # closer to arc than to its endpoints,

694 # find the closest point on the arc

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

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

697 n = gc1.cross(gc2)

698

699 elif within: # for backward compatibility, XXX unwrapped

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

701 self.distanceTo(point2)) else point2

702

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

704 n1 = p1.toNvector()

705 n2 = p2.toNvector()

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

707 if n is n1:

708 return p1 # is point1

709 elif n is n2:

710 return p2 # is point2 if not wrap

711

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

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

714 d = p1.distanceTo(p2)

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

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

717 return p

718

719 # @deprecated_method

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

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

722

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

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

725 to that point from this point ...

726 '''

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

728 return r.closest, r.distance

729

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

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

732 arcs joining consecutive points) closest to this point.

733

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

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

736

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

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

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

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

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

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

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

744

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

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

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

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

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

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

751

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

753

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

755 '''

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

757 _r = self.distanceTo

758 _n = self.nearestOn

759

760 c = p1 = Ps[0]

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

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

763 if wrap and not Ps.looped:

764 p2 = _unrollon(p1, p2)

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

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

767 if d < r:

768 c, r = p, d

769 p1 = p2

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

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

772

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

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

775

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

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

778 to override this L{Cartesian} class or specify

779 C{B{Cartesian}=None}.

780

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

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

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

784

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

786 '''

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

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

789

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

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

792

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

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

795 override this L{Nvector} class or use

796 C{B{Nvector}=None}.

797

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

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

800

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

802

803 @example:

804

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

806 >>> n = p.toNvector()

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

808 '''

809 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector)

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

811

812

813 class Nvector(NvectorBase):

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

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

816 n-vectors have no singularities or discontinuities.

817

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

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

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

821

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

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

824

825 Note commonality with L{ellipsoidalNvector.Nvector}.

826 '''

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

828

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

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

831 (ECEF) coordinates.

832

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

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

835 to override this L{Cartesian} class or specify

836 C{B{Cartesian}=None}.

837

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

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

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

841

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

843 '''

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

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

846

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

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

849

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

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

852 to override this L{LatLon} class or specify

853 C{B{LatLon}=None}.

854

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

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

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

858

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

860

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

862 '''

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

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

865

866 def greatCircle(self, bearing):

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

868 heading on given compass bearing from this point as its

869 n-vector.

870

871 Direction of vector is such that initial bearing vector

872 b = c × p.

873

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

875

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

877

878 @raise Valuerror: Polar coincidence.

879

880 @example:

881

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

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

884 '''

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

886

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

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

889

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

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

892 return n.minus(e)

893

894

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

896

897

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

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

900 great circle arcs joining consecutive points).

901

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

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

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

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

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

907

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

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

910

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

912

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

914

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

916 and L{ellipsoidalKarney.areaOf}.

917

918 @example:

919

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

921 >>> areaOf(b) # 8666058750.718977

922 '''

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

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

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

926 n0 = Ps[0].toNvector()

927

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

929 p1 = Ps[1]

930 v1 = p1._N_vector

931 gc = v2.cross(v1)

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

933 if w and not Ps.looped:

934 p2 = _unrollon(p1, p2)

935 p1 = p2

936 v2 = p2._N_vector

937 gc1 = v1.cross(v2)

938 v1 = v2

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

940 gc = gc1

941

942 if _MODS.booleans.isBoolean(points):

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

944 else:

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

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

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

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

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

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

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

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

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

954

955

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

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

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

959

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

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

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

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

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

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

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

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

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

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

970 circle} methods.

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

972 conventionally) or C{None}.

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

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

975

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

977 an instance of this class. For a tangent line, each point C{is}

978 this very instance.

979

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

981

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

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

984

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

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

987 '''

988 c = _Nvll.others(center=center)

989 p = _Nvll.others(point=point)

990 try:

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

992 height=height, wrap=wrap)

993 except (TypeError, ValueError) as x:

994 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing,

995 exact=exact, wrap=wrap)

996

997

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

999 LatLon=LatLon, **LatLon_kwds):

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

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

1002

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

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

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

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

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

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

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

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

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

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

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

1014 (L{LatLon}).

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

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

1017

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

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

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

1021

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

1023

1024 @raise ValueError: Intersection is ambiguous or infinite or

1025 the lines are parallel, coincident or null.

1026

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

1028

1029 @example:

1030

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

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

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

1034 '''

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

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

1037 return i.toLatLon(**kwds)

1038

1039

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

1041 LatLon=LatLon, **LatLon_kwds):

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

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

1044

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

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

1047 initial bearing at the first start point

1048 (compass C{degrees360}).

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

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

1051 initial bearing at the second start point

1052 (compass C{degrees360}).

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

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

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

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

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

1058 antipodal points (L{LatLon}).

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

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

1061

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

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

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

1065 if available.

1066

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

1068

1069 @raise ValueError: Intersection is ambiguous or infinite or

1070 the lines are parallel, coincident or null.

1071

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

1073 '''

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

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

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

1077

1078

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

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

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

1082 '''

1083 p1 = _Nvll.others(start1=start1)

1084 p2 = _Nvll.others(start2=start2)

1085 if wrap:

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

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

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

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

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

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

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

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

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

1095

1096 hs = start1.height, start2.height

1097 # there are two (antipodal) candidate intersection

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

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

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

1101

1102 # selection of intersection point depends on how

1103 # lines are defined (by bearings or endpoints)

1104 if e1 and e2: # endpoint+endpoint

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

1106 hs += end1.height, end2.height

1107 elif e1 and not e2: # endpoint+bearing

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

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

1110 hs += end1.height,

1111 elif e2 and not e1: # bearing+endpoint

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

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

1114 hs += end2.height,

1115 else: # bearing+bearing

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

1117 # towards i1, otherwise towards antipodal i2

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

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

1120 if d1 > 0 and d2 > 0:

1121 d = 1 # both point to i1

1122 elif d1 < 0 and d2 < 0:

1123 d = -1 # both point to i2

1124 else: # d1, d2 opposite signs

1125 # intersection is at further-away intersection point,

1126 # take opposite intersection from mid- point of v1

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

1128 # get intersection p1 bearing points to, aka being

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

1130 # function .sphericalTrigonometry._intersect and

1131 # .ellipsoidalBaseDI._intersect3

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

1133

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

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

1136

1137

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

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

1140

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

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

1143 (C{meter}).

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

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

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

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

1148

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

1150

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

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

1153 '''

1154 def _N_vs(ps, w):

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

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

1157 yield p._N_vector

1158

1159 try:

1160 # geographic mean

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

1162 print(n)

1163 except (TypeError, ValueError) as x:

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

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

1166 name=meanOf.__name__))

1167

1168

1169@deprecated_function

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

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

1172

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

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

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

1176 '''

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

1178 return r.closest, r.distance

1179

1180

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

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

1183 consecutive points) closest to an other point.

1184

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

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

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

1188

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

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

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

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

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

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

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

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

1197

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

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

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

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

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

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

1204 C{degrees360}.

1205

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

1207

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

1209 '''

1210 _xinstanceof(LatLon, point=point)

1211

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

1213 height=height, wrap=wrap)

1214

1215

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

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

1218 great circle arcs joining consecutive points).

1219

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

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

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

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

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

1225

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

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

1228

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

1230

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

1232

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

1234 C{B{points}} a composite.

1235

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

1237 and L{sphericalTrigonometry.perimeterOf}.

1238 '''

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

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

1241 p1 = Ps[0]

1242 v1 = p1._N_vector

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

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

1245 p2 = _unrollon(p1, p2)

1246 p1 = p2

1247 v2 = p2._N_vector

1248 yield v1.angleTo(v2)

1249 v1 = v2

1250

1251 if _MODS.booleans.isBoolean(points):

1252 if not closed:

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

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

1255 else:

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

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

1258

1259

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

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

1262

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

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

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

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

1267

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

1269

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

1271 '''

1272 try:

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

1274 except (TypeError, ValueError) as x:

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

1276

1277

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

1279 height=None, wrap=False,

1280 LatLon=LatLon, **LatLon_kwds):

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

1282 from those points.

1283

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

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

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

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

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

1289 the mean height (C{meter}).

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

1291 (C{bool}).

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

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

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

1295

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

1297

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

1299

1300 @raise Valuerror: Points coincide.

1301

1302 @example:

1303

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

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

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

1307 '''

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

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

1310 height=height, wrap=wrap,

1311 LatLon=LatLon, **LatLon_kwds)

1312

1313

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

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

1316 LatLon=LatLon, **LatLon_kwds):

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

1318

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

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

1321 as B{C{radius}}).

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

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

1324 as B{C{radius}}).

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

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

1327 as B{C{radius}}).

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

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

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

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

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

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

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

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

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

1337

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

1339

1340 @raise IntersectionError: No intersection, trilateration failed.

1341

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

1343

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

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

1346

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

1348 '''

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

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

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

1352 radius=radius, height=height, useZ=useZ,

1353 wrap=wrap, LatLon=LatLon, **LatLon_kwds)

1354

1355

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

1357 areaOf, # functions

1358 intersecant2, intersection, intersection2, ispolar,

1359 meanOf,

1360 nearestOn2, nearestOn3,

1361 perimeterOf,

1362 sumOf,

1363 triangulate, trilaterate)

1364

1365# **) MIT License

1366#

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

1368#

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

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

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

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

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

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

1375#

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

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

1378#

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

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

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

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

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

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

1385# OTHER DEALINGS IN THE SOFTWARE.