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 C{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-2024},

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 _xinstanceof # _MODS

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 fdot_, 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 .points

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

48 from pygeodesy.nvectorBase import LatLonNvectorBase, NorthPole, _nsumOf, \

49 NvectorBase, _triangulate, _trilaterate

50 from pygeodesy.points import NearestOn3Tuple, notImplemented, \

51 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_, _isDegrees, 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__ = '24.11.07'

65

66_lines_ = 'lines'

67

68

69 class Cartesian(CartesianSphericalBase):

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

71 C{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 C{LatLon} class and C{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: A C{LatLon} or if C{LatLon is None}, an L{Ecef9Tuple}C{(x, y, z,

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

84

85 @raise TypeError: Invalid C{LatLon} or other B{C{LatLon_and_kwds}} item.

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 C{Nvector} components, I{including height}.

92

93 @kwarg Nvector_and_kwds: Optional C{Nvector} class and C{Nvector} keyword

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

95 to override this C{Nvector} class or specify

96 C{B{Nvector}=None}.

97

98 @return: An C{Nvector}) or if C{Nvector is None}, a L{Vector4Tuple}C{(x, y, z, h)}.

99

100 @raise TypeError: Invalid C{Nvector} or other B{C{Nvector_and_kwds}} item.

101 '''

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

103 # datum=datum or self.datum)

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

105 # return ll.toNvector(**kwds)

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

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

108

109

110 class LatLon(LatLonNvectorBase, LatLonSphericalBase):

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

112

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

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

115 '''

116 _Nv = None # cached_toNvector C{Nvector})

117

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

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

120 '''

121 if updated: # reset caches

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

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

124

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

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

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

128 end point.

129

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

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

132 from the start point to the point where the perpendicular

133 crosses the line.

134

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

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

137 initial bearing from start point (compass

138 C{degrees360}).

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

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

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

142

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

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

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

146 toward the end point of the line or negative if

147 "before" the start point.

148

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

150

151 @raise Valuerror: Some points coincide.

152 '''

153 p = self.others(start=start)

154 n = self.toNvector()

155

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

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

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

159

160 @deprecated_method

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

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

163 '''

164 return self.initialBearingTo(other)

165

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

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

168 defined by a start and end point.

169

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

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

172 bearing from start point (compass C{degrees360}).

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

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

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

176

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

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

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

180 of the line .

181

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

183

184 @raise Valuerror: Some points coincide.

185 '''

186 p = self.others(start=start)

187 n = self.toNvector()

188

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

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

191

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

193 '''Locate the destination from this point after having travelled

194 the given distance on the given bearing.

195

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

197 B{C{radius}}).

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

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

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

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

202

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

204

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

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

207 '''

208 b = Bearing_(bearing)

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

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

211

212 n = self.toNvector()

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

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

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

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

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

218

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

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

221

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

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

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

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

226

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

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

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

230

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

232 '''

233 p = self.others(other)

234 if wrap:

235 p = _unrollon(self, p)

236 n = p.toNvector()

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

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

239

240# @Property_RO

241# def Ecef(self):

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

243# '''

244# return _ALL_MODS.ecef.EcefKarney

245

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

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

248 '''

249 s = start.toNvector()

250 if _isDegrees(end): # bearing

251 gc = s.greatCircle(end)

252 e = None

253 else: # point

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

255 if wrap:

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

257 e = p.toNvector()

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

259 return gc, s, e

260

261 def greatCircle(self, bearing):

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

263 heading on the given bearing from this point.

264

265 Direction of vector is such that initial bearing vector

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

267

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

269

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

271 '''

272 t = Bearing_(bearing)

273 a, b = self.philam

274

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

276

277 sa *= st

278 return Nvector(fdot_(sb, ct, -sa, cb),

279 -fdot_(cb, ct, sa, sb),

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

281

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

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

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

285

286 Direction of vector is such that initial bearing vector

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

288

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

290 this point (compass C{degrees360}).

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

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

293

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

295

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

297

298 @raise Valuerror: Points coincide.

299 '''

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

301 return gc.unit()

302

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

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

305 to an other point.

306

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

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

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

310

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

312

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

314 point or the C{NorthPole} and L{crosserrors

315 <pygeodesy.crosserrors>} is C{True}.

316

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

318 '''

319 n = self.toNvector()

320 p = self.others(other)

321 if wrap:

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

323 p = p.toNvector()

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

325# gc1 = self.greatCircleTo(other)

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

327# gc2 = self.greatCircleTo(NorthPole)

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

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

330

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

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

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

334

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

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

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

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

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

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

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

342

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

344

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

346 '''

347 n = self.toNvector()

348 p = self.others(other)

349 if wrap:

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

351

352 f = Scalar(fraction=fraction)

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

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

355

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

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

358

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

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

361 other point.

362

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

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

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

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

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

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

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

370

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

372

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

374

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

376

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

378 '''

379 n = self.toNvector()

380 p = self.others(other)

381 if wrap:

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

383 p = p.toNvector()

384 f = Scalar(fraction=fraction)

385

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

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

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

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

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

391

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

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

394

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

396 '''Locate an intersection point of two lines each defined by two

397 points or by a point and an (initial) bearing.

398

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

400

401 @see: Method L{intersection2<sphericalNvector.LatLon.intersection2>}

402 for further details.

403 '''

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

405 wrap=wrap, LatLon=self.classof)

406

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

408 '''Locate both intersections of two (great circle) lines each defined

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

410

411 @arg end1: End point of the line starting at this point (L{LatLon})

412 or the bearing at this point (compass C{degrees360}).

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

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

415 at B{C{start2}} (compass C{degrees360}).

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

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

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

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

420

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

422

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

424 point is not L{LatLon}.

425

426 @raise ValueError: Intersection is ambiguous or infinite or

427 the lines are parallel, coincident or null.

428

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

430 '''

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

432 wrap=wrap, LatLon=self.classof)

433

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

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

436

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

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

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

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

441

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

443 C{False} otherwise.

444

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

446

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

448

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

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

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

452 '''

453 if _MODS.booleans.isBoolean(points):

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

455

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

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

458 def _subt(Ps, n0, w):

459 p1 = Ps[0]

460 vs1 = n0.minus(p1.toNvector())

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

462 if w and not Ps.looped:

463 p2 = _unrollon(p1, p2)

464 p1 = p2

465 vs2 = n0.minus(p2.toNvector())

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

467 vs1 = vs2

468

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

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

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

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

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

474 s = fsum(_subt(self.PointsIter(points, loop=1, wrap=wrap),

475 self.toNvector(), wrap)) # normal vector

476 # XXX are winding number optimisations equally applicable to

477 # spherical surface?

478 return fabs(s) > PI

479

480 @deprecated_method

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

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

483 return self.isenclosedBy(points)

484

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

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

487

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

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

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

491 the same hemispere).

492

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

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

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

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

497

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

499 arc, C{False} otherwise.

500

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

502 L{LatLon}.

503 '''

504 p1 = self.others(point1=point1)

505 p2 = self.others(point2=point2)

506 if wrap:

507 p1 = _Wrap.point(p1)

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

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

510

511 # corner case, null arc

512 if n1.isequalTo(n2):

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

514

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

516 return False # PYCHOK returns

517

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

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

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

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

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

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

524

525 @deprecated_method

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

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

528 return self.iswithin(point1, point2)

529

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

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

532

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

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

535 the mean height (C{meter}).

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

537 may be negative or greater than 1.0.

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

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

540

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

542

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

544

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

546 '''

547 if fraction is _0_5:

548 p = self.others(other)

549 if wrap:

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

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

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

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

554 else:

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

556 return r

557

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

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

560 closest to this point.

561

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

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

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

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

566 to interpolate the height.

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

568 given points, otherwise the closest point

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

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

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

572

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

574

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

576 not supported.

577

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

579 '''

580 p1 = self.others(point1=point1)

581 p2 = self.others(point2=point2)

582 if wrap:

583 p1 = _Wrap.point(p1)

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

585 p0 = self

586

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

588 # closer to arc than to its endpoints,

589 # find the closest point on the arc

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

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

592 n = gc1.cross(gc2)

593

594 elif within: # for backward compatibility, XXX unwrapped

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

596 self.distanceTo(point2)) else point2

597

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

599 n1 = p1.toNvector()

600 n2 = p2.toNvector()

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

602 if n is n1:

603 return p1 # is point1

604 elif n is n2:

605 return p2 # is point2 if not wrap

606

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

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

609 d = p1.distanceTo(p2)

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

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

612 return p

613

614 # @deprecated_method

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

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

617

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

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

620 to that point from this point ...

621 '''

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

623 return r.closest, r.distance

624

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

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

627 joining consecutive points) closest to this point.

628

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

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

631

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

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

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

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

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

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

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

639

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

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

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

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

644 C{B{radius} is None}) and the C{angle} from this to

645 the C{closest} point in compass C{degrees360}.

646

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

648

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

650 '''

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

652 _d = self.distanceTo

653 _n = self.nearestOn

654

655 c = p1 = Ps[0]

656 r = _d(c, radius=None) # radians

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

658 if wrap and not Ps.looped:

659 p2 = _unrollon(p1, p2)

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

661 d = _d(p, radius=None) # radians

662 if d < r:

663 c, r = p, d

664 p1 = p2

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

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

667

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

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

670

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

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

673 to override this L{Cartesian} class or specify

674 C{B{Cartesian}=None}.

675

676 @return: A L{Cartesian} or if C{B{Cartesian} is None}, an L{Ecef9Tuple}C{(x, y,

677 z, lat, lon, height, C, M, datum)} with C{C} and C{M} if available.

678

679 @raise TypeError: Invalid L{Cartesian} or other B{C{Cartesian_and_kwds}} item.

680 '''

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

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

683

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

685 '''Convert this point to C{Nvector} components, I{including height}.

686

687 @kwarg Nvector_and_kwds: Optional C{Nvector} and C{Nvector} keyword arguments.

688 Specify C{B{Nvector}=...} to override this C{Nvector}

689 class or use C{B{Nvector}=None}.

690

691 @return: An C{Nvector} or if C{B{Nvector} is None}, a L{Vector4Tuple}C{(x, y, z, h)}.

692

693 @raise TypeError: Invalid C{Nvector} or other B{C{Nvector_and_kwds}} item.

694 '''

695 return LatLonNvectorBase.toNvector(self, **_xkwds(Nvector_and_kwds, Nvector=Nvector))

696

697

698 class Nvector(NvectorBase):

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

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

701 n-vectors have no singularities or discontinuities.

702

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

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

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

706

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

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

709

710 Note commonality with L{pygeodesy.ellipsoidalNvector.Nvector}.

711 '''

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

713

714 @property_RO

715 def sphericalNvector(self):

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

717 '''

718 return type(self)

719

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

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

722 (ECEF) coordinates.

723

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

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

726 to override this L{Cartesian} class or specify

727 C{B{Cartesian}=None}.

728

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

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

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

732

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

734 '''

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

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

737

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

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

740

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

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

743 to override this L{LatLon} class or specify

744 C{B{LatLon}=None}.

745

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

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

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

749

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

751

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

753 '''

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

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

756

757 def greatCircle(self, bearing):

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

759 heading on given (initial) bearing from this point as its

760 n-vector.

761

762 Direction of vector is such that initial bearing vector

763 b = c × p.

764

765 @arg bearing: Initial bearing (compass C{degrees360}).

766

767 @return: N-vector representing great circle (C{Nvector}).

768

769 @raise Valuerror: Polar coincidence.

770 '''

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

772

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

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

775

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

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

778 return n.minus(e)

779

780

781_Nv00 = LatLon(_0_0, _0_0, name=_Nv00_) # reference instance (L{LatLon})

782

783

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

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

786 great circle arcs joining consecutive points).

787

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

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

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

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

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

793

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

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

796

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

798

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

800

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

802 and L{ellipsoidalKarney.areaOf}.

803 '''

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

805 Ps = _Nv00.PointsIter(ps, loop=2, wrap=w)

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

807 n0 = Ps[0].toNvector()

808

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

810 p1 = Ps[1]

811 v1 = p1._N_vector

812 gc = v2.cross(v1)

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

814 if w and not Ps.looped:

815 p2 = _unrollon(p1, p2)

816 p1 = p2

817 v2 = p2._N_vector

818 gc1 = v1.cross(v2)

819 v1 = v2

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

821 gc = gc1

822

823 if _MODS.booleans.isBoolean(points):

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

825 else:

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

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

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

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

830 s = fsum(_interangles(points, wrap))

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

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

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

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

835

836

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

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

839 two points or as a point and bearing.

840

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

842 @arg circle: Radius of the circle (C{meter}, same units as the earth

843 B{C{radius}}) or a point on the circle (L{LatLon}).

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

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

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

847 @kwarg radius_exact_height_wrap: Optional keyword arguments, see method

848 L{intersecant2<pygeodesy.sphericalBase.LatLonSphericalBase.

849 intersecant2>} for further details.

850

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

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

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

854

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

856

857 @raise TypeError: If B{C{center}}, B{C{point}}, B{C{circle}} or B{C{other}}

858 not L{LatLon}.

859

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

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

862 '''

863 c = _Nv00.others(center=center)

864 p = _Nv00.others(point=point)

865 try:

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

867 except (TypeError, ValueError) as x:

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

869 **radius_exact_height_wrap)

870

871

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

873 **LatLon_and_kwds):

874 '''Locate an intersection point of two (great circle) lines each defined

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

876

877 @return: The intersection point (L{LatLon}) or if C{B{LatLon}=None},

878 a cartesian L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,

879 datum)} with C{C} and C{M} if available.

880

881 @see: Function L{intersection2<sphericalNvector.intersection2>}

882 for further details.

883 '''

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

885 kwds = _xkwds(LatLon_and_kwds, height=h, LatLon=LatLon)

886 return i.toLatLon(**kwds)

887

888

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

890 **LatLon_and_kwds):

891 '''Locate both intersections of two (great circle) lines each defined

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

893

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

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

896 B{C{start1}} (compass C{degrees360}).

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

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

899 B{C{start2}} (compass C{degrees360}).

900 @kwarg height: Optional height at the intersection and antipodal point,

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

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

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

904 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=}L{LatLon} to return

905 the intersection points and optionally, additional B{C{LatLon}}

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

907

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

909 C{B{LatLon}=None}, a cartesian L{Ecef9Tuple}C{(x, y, z, lat, lon,

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

911

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

913

914 @raise ValueError: Intersection is ambiguous or infinite or the lines are

915 parallel, coincident or null.

916

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

918 '''

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

920 kwds = _xkwds(LatLon_and_kwds, height=h, LatLon=LatLon)

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

922

923

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

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

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

927 '''

928 p1 = _Nv00.others(start1=start1)

929 p2 = _Nv00.others(start2=start2)

930 if wrap:

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

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

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

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

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

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

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

938 gc1, s1, e1 = _Nv00._gc3(p1, end1, 'end1', wrap=wrap)

939 gc2, s2, e2 = _Nv00._gc3(p2, end2, 'end2', wrap=wrap)

940

941 hs = start1.height, start2.height

942 # there are two (antipodal) candidate intersection

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

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

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

946

947 # selection of intersection point depends on how

948 # lines are defined (by bearings or endpoints)

949 if e1 and e2: # endpoint+endpoint

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

951 hs += end1.height, end2.height

952 elif e1 and not e2: # endpoint+bearing

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

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

955 hs += end1.height,

956 elif e2 and not e1: # bearing+endpoint

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

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

959 hs += end2.height,

960 else: # bearing+bearing

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

962 # towards i1, otherwise towards antipodal i2

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

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

965 if d1 > 0 and d2 > 0:

966 d = 1 # both point to i1

967 elif d1 < 0 and d2 < 0:

968 d = -1 # both point to i2

969 else: # d1, d2 opposite signs

970 # intersection is at further-away intersection point,

971 # take opposite intersection from mid- point of v1

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

973 # get intersection p1 bearing points to, aka being

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

975 # function .sphericalTrigonometry._intersect and

976 # .ellipsoidalBaseDI._intersect3

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

978

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

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

981

982

983 def meanOf(points, height=None, wrap=False, **LatLon_and_kwds):

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

985

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

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

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

989 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=}L{LatLon} to return

990 the mean point and optionally, additional B{C{LatLon}}

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

992

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

994

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

996 B{C{points}} are not C{LatLon}.

997 '''

998 def _N_vs(ps, w):

999 Ps = _Nv00.PointsIter(ps, wrap=w)

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

1001 yield p._N_vector

1002

1003 try:

1004 # geographic mean

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

1006 except (TypeError, ValueError) as x:

1007 raise PointsError(points=points, wrap=wrap, cause=x, **LatLon_and_kwds)

1008 return n.toLatLon(**_xkwds(LatLon_and_kwds, LatLon=LatLon, height=n.h,

1009 name=meanOf.__name__))

1010

1011

1012@deprecated_function

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

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

1015

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

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

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

1019 '''

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

1021 return r.closest, r.distance

1022

1023

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

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

1026 consecutive points) closest to an other point.

1027

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

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

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

1031

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

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

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

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

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

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

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

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

1040

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

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

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

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

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

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

1047 C{degrees360}.

1048

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

1050

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

1052 '''

1053 _MODS.basics._xinstanceof(LatLon, point=point)

1054

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

1056 height=height, wrap=wrap)

1057

1058

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

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

1061 great circle arcs joining consecutive points).

1062

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

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

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

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

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

1068

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

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

1071

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

1073

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

1075

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

1077 C{B{points}} a composite.

1078

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

1080 and L{sphericalTrigonometry.perimeterOf}.

1081 '''

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

1083 Ps = _Nv00.PointsIter(ps, loop=1, wrap=w)

1084 p1 = Ps[0]

1085 v1 = p1._N_vector

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

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

1088 p2 = _unrollon(p1, p2)

1089 p1 = p2

1090 v2 = p2._N_vector

1091 yield v1.angleTo(v2)

1092 v1 = v2

1093

1094 if _MODS.booleans.isBoolean(points):

1095 if not closed:

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

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

1098 else:

1099 r = fsum(_rads(points, closed, wrap))

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

1101

1102

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

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

1105

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

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

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

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

1110

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

1112

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

1114 '''

1115 try:

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

1117 except (TypeError, ValueError) as x:

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

1119

1120

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

1122 height=None, wrap=False,

1123 LatLon=LatLon, **LatLon_kwds):

1124 '''Locate a point given two known, reference points and the (initial)

1125 bearing from those points.

1126

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

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

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

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

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

1132 the mean height (C{meter}).

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

1134 (C{bool}).

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

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

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

1138

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

1140

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

1142

1143 @raise Valuerror: Points coincide.

1144 '''

1145 return _triangulate(_Nv00.others(point1=point1), bearing1,

1146 _Nv00.others(point2=point2), bearing2,

1147 height=height, wrap=wrap,

1148 LatLon=LatLon, **LatLon_kwds)

1149

1150

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

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

1153 LatLon=LatLon, **LatLon_kwds):

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

1155

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

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

1158 as B{C{radius}}).

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

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

1161 as B{C{radius}}).

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

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

1164 as B{C{radius}}).

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

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

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

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

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

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

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

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

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

1174

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

1176

1177 @raise IntersectionError: No intersection, trilateration failed.

1178

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

1180

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

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

1183

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

1185 '''

1186 return _trilaterate(_Nv00.others(point1=point1), distance1,

1187 _Nv00.others(point2=point2), distance2,

1188 _Nv00.others(point3=point3), distance3,

1189 radius=radius, height=height, useZ=useZ,

1190 wrap=wrap, LatLon=LatLon, **LatLon_kwds)

1191

1192

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

1194 areaOf, # functions

1195 intersecant2, intersection, intersection2, ispolar,

1196 meanOf,

1197 nearestOn2, nearestOn3,

1198 perimeterOf,

1199 sumOf,

1200 triangulate, trilaterate)

1201

1202# **) MIT License

1203#

1204# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1205#

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

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

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

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

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

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

1212#

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

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

1215#

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

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

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

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

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

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

1222# OTHER DEALINGS IN THE SOFTWARE.