Coverage for pygeodesy/sphericalNvector.py: 97%

259 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-03-31 10:52 -0400

1 

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

3 

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

5 

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

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

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

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

10 

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

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

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

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

15U{Module latlon-nvector-spherical 

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

17 

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

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

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

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

22 

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

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

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

26 

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

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

29 

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

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

32''' 

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

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

35 

36from pygeodesy.basics import isscalar, _xinstanceof 

37from 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 

40from pygeodesy.errors import _ValueError, _xError, _xkwds 

41from pygeodesy.fmath import fmean, fsum 

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

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

44 _points_, _pole_ 

45from 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 

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

49 sumOf as _sumOf, _triangulate, _trilaterate, \ 

50 notImplemented 

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

52from pygeodesy.props import deprecated_function, deprecated_method 

53from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, \ 

54 Datums, _intersecant2, LatLonSphericalBase 

55from pygeodesy.units import Bearing, Bearing_, Height, Radius, Scalar 

56from pygeodesy.utily import atan2, degrees360, fabs, sincos2, sincos2_, sincos2d 

57 

58# from math import atan2, fabs # from utily 

59 

60__all__ = _ALL_LAZY.sphericalNvector 

61__version__ = '23.03.30' 

62 

63_paths_ = 'paths' 

64 

65 

66class Cartesian(CartesianSphericalBase): 

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

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

69 ''' 

70 

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

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

73 

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

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

76 to override this L{LatLon} class or specify 

77 C{B{LatLon}=None}. 

78 

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

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

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

82 

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

84 ''' 

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

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

87 

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

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

90 

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

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

93 to override this L{Nvector} class or specify 

94 C{B{Nvector}=None}. 

95 

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

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

98 

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

100 ''' 

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

102 # datum=datum or self.datum) 

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

104 # return ll.toNvector(**kwds) 

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

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

107 

108 

109class LatLon(LatLonNvectorBase, LatLonSphericalBase): 

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

111 

112 Tools for working with points and paths on (a spherical 

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

114 

115 @example: 

116 

117 >>> from sphericalNvector import LatLon 

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

119 ''' 

120 _Nv = None # cached_toNvector L{Nvector}) 

121 

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

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

124 ''' 

125 if updated: # reset caches 

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

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

128 

129 def alongTrackDistanceTo(self, start, end, radius=R_M): 

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

131 point on the great circle path defined by a start and an 

132 end point. 

133 

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

135 great circle path, the along-track distance is the distance 

136 from the start point to the point where the perpendicular 

137 crosses the path. 

138 

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

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

141 initial bearing from start point (compass 

142 C{degrees360}). 

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

144 

145 @return: Distance along the great circle path (positive if 

146 after the start toward the end point of the path 

147 or negative if 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 @example: 

154 

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

156 

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

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

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

160 ''' 

161 self.others(start=start) 

162 gc, _, _ = self._gc3(start, end, _end_) 

163 

164 p = self.toNvector() 

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

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

167 

168 @deprecated_method 

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

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

171 ''' 

172 return self.initialBearingTo(other) 

173 

174 def crossTrackDistanceTo(self, start, end, radius=R_M): 

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

176 defined by a start and end point. 

177 

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

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

180 initial bearing from start point (compass 

181 C{degrees360}). 

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

183 

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

185 left or positive if to the right of the path). 

186 

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

188 

189 @raise Valuerror: Some points coincide. 

190 

191 @example: 

192 

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

194 

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

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

197 

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

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

200 ''' 

201 self.others(start=start) 

202 gc, _, _ = self._gc3(start, end, _end_) 

203 

204 p = self.toNvector() 

205 return (gc.angleTo(p) - PI_2) * radius 

206 

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

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

209 travelled the given distance on the given bearing. 

210 

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

212 as B{C{radius}}). 

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

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

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

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

217 

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

219 

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

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

222 

223 @example: 

224 

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

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

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

228 ''' 

229 a = _angular(distance, radius) 

230 sa, ca, sb, cb = sincos2_(a, Bearing_(bearing)) 

231 

232 p = self.toNvector() 

233 e = NorthPole.cross(p, raiser=_pole_).unit() # east vector at p 

234 n = p.cross(e) # north vector at p 

235 q = n.times(cb).plus(e.times(sb)) # direction vector @ p 

236 n = p.times(ca).plus(q.times(sa)) 

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

238 

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

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

241 

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

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

244 @kwarg wrap: Wrap/unroll the angular distance (C{bool}). 

245 

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

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

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

249 

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

251 

252 @example: 

253 

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

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

256 >>> d = p.distanceTo(q) # 404.3 km 

257 ''' 

258 self.others(other) 

259 

260 r = fabs(self.toNvector().angleTo(other.toNvector(), wrap=wrap)) 

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

262 

263# @Property_RO 

264# def Ecef(self): 

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

266# ''' 

267# return _ALL_MODS.ecef.EcefKarney 

268 

269 def _gc3(self, start, end, namend, raiser=_points_): 

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

271 ''' 

272 s = start.toNvector() 

273 if isscalar(end): # bearing 

274 gc = s.greatCircle(end) 

275 e = None 

276 else: 

277 self.others(end, name=namend) 

278 e = end.toNvector() 

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

280 return gc, s, e 

281 

282 def greatCircle(self, bearing): 

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

284 heading on the given bearing from this point. 

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 bearing: Bearing from this point (compass C{degrees360}). 

290 

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

292 ''' 

293 a, b = self.philam 

294 t = Bearing_(bearing) 

295 

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

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

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

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

300 

301 def greatCircleTo(self, other): 

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

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

304 

305 Direction of vector is such that initial bearing vector 

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

307 

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

309 this point (compass C{degrees360}). 

310 

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

312 

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

314 

315 @raise Valuerror: Points coincide. 

316 

317 @example: 

318 

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

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

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

322 

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

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

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

326 ''' 

327 gc, _, _ = self._gc3(self, other, _other_) 

328 return gc.unit() 

329 

330 def initialBearingTo(self, other, **unused): 

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

332 to an other point. 

333 

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

335 @arg unused: Optional keyword argument B{C{wrap}} ignored. 

336 

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

338 

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

340 point or the C{NorthPole}, provided 

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

342 

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

344 

345 @example: 

346 

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

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

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

350 ''' 

351 self.others(other) 

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

353 n = self.toNvector() 

354# gc1 = self.greatCircleTo(other) 

355 gc1 = n.cross(other.toNvector(), raiser=_points_) # .unit() 

356# gc2 = self.greatCircleTo(NorthPole) 

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

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

359 

360 def intermediateChordTo(self, other, fraction, height=None): 

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

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

363 

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

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

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

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

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

369 

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

371 

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

373 

374 @example: 

375 

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

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

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

379 ''' 

380 self.others(other) 

381 

382 f = Scalar(fraction=fraction) 

383 i = other.toNvector().times(f).plus( 

384 self.toNvector().times(1 - f)) 

385# i = other.toNvector() * f + \ 

386# self.toNvector() * (1 - f)) 

387 

388 h = self._havg(other, f=f) if height is None else Height(height) 

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

390 

391 def intermediateTo(self, other, fraction, height=None, **unused): # wrap=False 

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

393 other point. 

394 

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

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

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

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

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

400 

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

402 

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

404 

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

406 

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

408 

409 @example: 

410 

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

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

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

414 ''' 

415 q = self.others(other).toNvector() 

416 p = self.toNvector() 

417 f = Scalar(fraction=fraction) 

418 

419 x = p.cross(q, raiser=_points_) 

420 d = x.unit().cross(p) # unit(p × q) × p 

421 # angular distance α, tan(α) = |p × q| / p ⋅ q 

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

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

424 

425 h = self._havg(other, f=f) if height is None else Height(height) 

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

427 

428 def intersection(self, end1, start2, end2, height=None): 

429 '''Locate the intersection point of two paths each defined 

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

431 

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

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

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

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

436 initial bearing at the second point (compass 

437 C{degrees}). 

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

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

440 

441 @return: The intersection point (L{LatLon}) or C{None} 

442 if no unique intersection exists. 

443 

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

445 is not L{LatLon}. 

446 

447 @raise ValueError: Intersection is ambiguous or infinite or 

448 the paths are parallel, coincident or null. 

449 

450 @example: 

451 

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

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

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

455 ''' 

456 return intersection(self, end1, start2, end2, 

457 height=height, LatLon=self.classof) 

458 

459 def isenclosedBy(self, points): 

460 '''Check whether this point is enclosed by a (convex) polygon. 

461 

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

463 

464 @return: C{True} if the polygon encloses this point, 

465 C{False} otherwise. 

466 

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

468 

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

470 

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

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

473 enclose a pole or wrap around the earth longitudinally. 

474 

475 @example: 

476 

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

478 >>> p = LatLon(45.1, 1.1) 

479 >>> inside = p.isenclosedBy(b) # True 

480 ''' 

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

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

483 def _subtangles(Ps, n0): 

484 vs1 = n0.minus(Ps[0].toNvector()) 

485 for p in Ps.iterate(closed=True): 

486 vs2 = n0.minus(p.toNvector()) 

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

488 vs1 = vs2 

489 

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

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

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

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

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

495 s = fsum(_subtangles(self.PointsIter(points, loop=1), 

496 self.toNvector()), floats=True) # normal vector 

497 # XXX are winding number optimisations equally applicable to 

498 # spherical surface? 

499 return fabs(s) > PI 

500 

501 @deprecated_method 

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

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

504 return self.isenclosedBy(points) 

505 

506 def iswithin(self, point1, point2): 

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

508 

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

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

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

512 the same hemispere). 

513 

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

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

516 

517 @return: C{True} if this point is within the arc, 

518 C{False} otherwise. 

519 

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

521 ''' 

522 n0 = self.toNvector() 

523 n1 = self.others(point1=point1).toNvector() 

524 n2 = self.others(point2=point2).toNvector() 

525 

526 # corner case, null arc 

527 if n1.isequalTo(n2): 

528 return n0.isequalTo(n1) or n0.isequalTo(n2) # PYCHOK returns 

529 

530 if n0.dot(n1) < 0 or n0.dot(n2) < 0: # different hemisphere 

531 return False # PYCHOK returns 

532 

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

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

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

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

537 return n0.minus(n1).dot(n2.minus(n1)) >= 0 and \ 

538 n0.minus(n2).dot(n1.minus(n2)) >= 0 

539 

540 @deprecated_method 

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

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

543 return self.iswithin(point1, point2) 

544 

545 def midpointTo(self, other, height=None, fraction=_0_5): 

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

547 

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

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

550 the mean height (C{meter}). 

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

552 may be negative or greater than 1.0. 

553 

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

555 

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

557 

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

559 

560 @example: 

561 

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

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

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

565 ''' 

566 if fraction is _0_5: 

567 self.others(other) 

568 

569 m = self.toNvector().plus(other.toNvector()) 

570 h = self._havg(other) if height is None else height 

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

572 else: 

573 r = self.intermediateTo(other, fraction, height=height) 

574 return r 

575 

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

577 '''Locate the point on the great circle arc between two 

578 points closest to this point. 

579 

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

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

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

583 for the point within the arc (C{meter}), or 

584 C{None} to interpolate the height. 

585 @kwarg within: If C{True} return the closest point between 

586 both given points, otherwise the closest 

587 point elsewhere on the arc (C{bool}). 

588 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

589 

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

591 

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

593 not supported. 

594 

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

596 

597 @example: 

598 

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

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

601 

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

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

604 

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

606 

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

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

609 ''' 

610 if wrap: # wrap=True throws C{NotImplementedError} always. 

611 notImplemented(self, wrap=wrap) 

612 

613 if self.iswithin(point1, point2) and not point1.isequalTo(point2, EPS): 

614 # closer to arc than to its endpoints, 

615 # find the closest point on the arc 

616 gc1 = point1.toNvector().cross(point2.toNvector()) 

617 gc2 = self.toNvector().cross(gc1) 

618 n = gc1.cross(gc2) 

619 

620 elif within: # for backward compatibility 

621 return point1 if self.distanceTo(point1) < self.distanceTo(point2) else point2 

622 

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

624 n1 = point1.toNvector() 

625 n2 = point2.toNvector() 

626 n = self.toNvector().nearestOn(n1, n2, within=False) 

627 if n is n1: 

628 return point1 

629 elif n is n2: 

630 return point2 

631 

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

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

634 d = point1.distanceTo(point2) 

635 f = (point1.distanceTo(p) / d) if d > EPS0 else _0_5 

636 p.height = point1._havg(point2, f=max(_0_0, min(f, _1_0))) 

637 return p 

638 

639 # @deprecated_method 

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

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

642 

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

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

645 to that point from this point ... 

646 ''' 

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

648 return r.closest, r.distance 

649 

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

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

652 arcs joining consecutive points) closest to this point. 

653 

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

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

656 

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

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

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

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

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

662 

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

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

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

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

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

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

669 

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

671 

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

673 ''' 

674 Ps = self.PointsIter(points, loop=1) 

675 

676 R = self.distanceTo 

677 N = self.nearestOn 

678 

679 c = p1 = Ps[0] 

680 r = R(c, radius=None) # radians 

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

682 p = N(p1, p2, height=height) 

683 d = R(p, radius=None) # radians 

684 if d < r: 

685 c, r = p, d 

686 p1 = p2 

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

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

689 

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

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

692 

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

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

695 to override this L{Cartesian} class or specify 

696 C{B{Cartesian}=None}. 

697 

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

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

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

701 

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

703 ''' 

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

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

706 

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

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

709 

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

711 arguments. Use C{B{Nvector}=...} to override 

712 this L{Nvector} class or specify 

713 C{B{Nvector}=None}. 

714 

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

716 set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)} if B{C{Nvector}}. 

717 

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

719 

720 @example: 

721 

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

723 >>> n = p.toNvector() 

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

725 ''' 

726 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector) 

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

728 

729 

730class Nvector(NvectorBase): 

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

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

733 n-vectors have no singularities or discontinuities. 

734 

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

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

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

738 

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

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

741 

742 Note commonality with L{ellipsoidalNvector.Nvector}. 

743 ''' 

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

745 

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

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

748 (ECEF) coordinates. 

749 

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

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

752 to override this L{Cartesian} class or specify 

753 C{B{Cartesian}=None}. 

754 

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

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

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

758 

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

760 ''' 

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

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

763 

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

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

766 

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

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

769 to override this L{LatLon} class or specify 

770 C{B{LatLon}=None}. 

771 

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

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

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

775 

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

777 

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

779 ''' 

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

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

782 

783 def greatCircle(self, bearing): 

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

785 heading on given compass bearing from this point as its 

786 n-vector. 

787 

788 Direction of vector is such that initial bearing vector 

789 b = c × p. 

790 

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

792 

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

794 

795 @raise Valuerror: Polar coincidence. 

796 

797 @example: 

798 

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

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

801 ''' 

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

803 

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

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

806 

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

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

809 return n.minus(e) 

810 

811 

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

813 

814 

815def areaOf(points, radius=R_M): 

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

817 (with great circle arcs joining consecutive points). 

818 

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

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

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

822 

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

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

825 

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

827 

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

829 

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

831 and L{ellipsoidalKarney.areaOf}. 

832 

833 @example: 

834 

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

836 >>> areaOf(b) # 8666058750.718977 

837 ''' 

838 def _interangles(Ps): 

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

840 n0 = Ps[0].toNvector() 

841 

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

843 v1 = Ps[1]._N_vector 

844 gc = v2.cross(v1) 

845 for p in Ps.iterate(closed=True): 

846 v2 = p._N_vector 

847 gc1 = v1.cross(v2) 

848 v1 = v2 

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

850 gc = gc1 

851 

852 if _MODS.booleans.isBoolean(points): 

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

854 else: 

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

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

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

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

859 s = fsum(_interangles(_Nvll.PointsIter(points, loop=2)), floats=True) 

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

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

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

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

864 

865 

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

867 height=None, wrap=True): 

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

869 

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

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

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

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

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

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

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

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

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

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

880 circle} methods. 

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

882 conventionally) or C{None}. 

883 @kwarg wrap: Wrap and unroll longitudes (C{bool}). 

884 

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

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

887 point C{is} this very instance. 

888 

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

890 

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

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

893 

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

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

896 ''' 

897 c = _Nvll.others(center=center) 

898 p = _Nvll.others(point=point) 

899 try: 

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

901 height=height, wrap=wrap) 

902 except (TypeError, ValueError) as x: 

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

904 

905 

906def intersection(start1, end1, start2, end2, 

907 height=None, LatLon=LatLon, **LatLon_kwds): 

908 '''Locate the intersection of two paths each defined by two 

909 points or by a start point and an initial bearing. 

910 

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

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

913 initial bearing at the first start point 

914 (compass C{degrees360}). 

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

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

917 initial bearing at the second start point 

918 (compass C{degrees360}). 

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

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

921 @kwarg LatLon: Optional class to return the intersection 

922 point (L{LatLon}). 

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

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

925 

926 @return: The intersection point (B{C{LatLon}}) or 3-tuple 

927 (C{degrees90}, C{degrees180}, height) if B{C{LatLon}} 

928 is C{None} or C{None} if no unique intersection 

929 exists. 

930 

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

932 

933 @raise ValueError: Intersection is ambiguous or infinite or 

934 the paths are parallel, coincident or null. 

935 

936 @example: 

937 

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

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

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

941 ''' 

942 _Nvll.others(start1=start1) 

943 _Nvll.others(start2=start2) 

944 

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

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

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

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

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

950 # paths are defined by start/end points, take closer intersection. 

951 gc1, s1, e1 = _Nvll._gc3(start1, end1, 'end1') 

952 gc2, s2, e2 = _Nvll._gc3(start2, end2, 'end2') 

953 

954 hs = start1.height, start2.height 

955 # there are two (antipodal) candidate intersection 

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

957 i1 = gc1.cross(gc2, raiser=_paths_) 

958 # postpone computing i2 until needed 

959 # i2 = gc2.cross(gc1, raiser=_paths_) 

960 

961 # selection of intersection point depends on how 

962 # paths are defined (by bearings or endpoints) 

963 if e1 and e2: # endpoint+endpoint 

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

965 hs += end1.height, end2.height 

966 elif e1 and not e2: # endpoint+bearing 

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

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

969 hs += end1.height, 

970 elif e2 and not e1: # bearing+endpoint 

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

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

973 hs += end2.height, 

974 else: # bearing+bearing 

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

976 # towards i1, otherwise towards antipodal i2 

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

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

979 if d1 > 0 and d2 > 0: 

980 d = 1 # both point to i1 

981 elif d1 < 0 and d2 < 0: 

982 d = -1 # both point to i2 

983 else: # d1, d2 opposite signs 

984 # intersection is at further-away intersection point, 

985 # take opposite intersection from mid- point of v1 

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

987 # get intersection p1 bearing points to, aka being 

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

989 # function .sphericalTrigonometry._intersect and 

990 # .ellipsoidalBaseDI._intersect3 

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

992 

993 i = i1 if d > 0 else gc2.cross(gc1, raiser=_paths_) 

994 

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

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

997 return i.toLatLon(**kwds) # Nvector(i.x, i.y, i.z).toLatLon(...) 

998 

999 

1000def meanOf(points, height=None, LatLon=LatLon, **LatLon_kwds): 

1001 '''Compute the geographic mean of the supplied points. 

1002 

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

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

1005 (C{meter}). 

1006 @kwarg LatLon: Optional class to return the mean point 

1007 (L{LatLon}). 

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

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

1010 

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

1012 

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

1014 

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

1016 ''' 

1017 Ps = _Nvll.PointsIter(points) 

1018 # geographic mean 

1019 m = sumOf(p._N_vector for p in Ps.iterate(closed=False)) 

1020 kwds = _xkwds(LatLon_kwds, height=height, LatLon=LatLon, 

1021 name=meanOf.__name__) 

1022 return m.toLatLon(**kwds) 

1023 

1024 

1025@deprecated_function 

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

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

1028 

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

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

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

1032 ''' 

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

1034 return r.closest, r.distance 

1035 

1036 

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

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

1039 joining consecutive points) closest to an other point. 

1040 

1041 If the given point is within the extent of any great circle 

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

1043 closest is the nearest of the arc's end points. 

1044 

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

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

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

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

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

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

1051 

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

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

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

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

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

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

1058 C{degrees360}. 

1059 

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

1061 

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

1063 ''' 

1064 _xinstanceof(LatLon, point=point) 

1065 

1066 return point.nearestOn3(points, closed=closed, radius=radius, height=height) 

1067 

1068 

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

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

1071 (with great circle arcs joining consecutive points). 

1072 

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

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

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

1076 

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

1078 C{radians} if B{C{radius}} is C{None}). 

1079 

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

1081 

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

1083 

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

1085 C{B{points}} a composite. 

1086 

1087 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalTrigonometry.perimeterOf} 

1088 and L{ellipsoidalKarney.perimeterOf}. 

1089 ''' 

1090 def _rads(Ps, closed): # angular edge lengths in radians 

1091 v1 = Ps[0]._N_vector 

1092 for p in Ps.iterate(closed=closed): 

1093 v2 = p._N_vector 

1094 yield v1.angleTo(v2) 

1095 v1 = v2 

1096 

1097 if _MODS.booleans.isBoolean(points): 

1098 if not closed: 

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

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

1101 else: 

1102 r = fsum(_rads(_Nvll.PointsIter(points, loop=1), closed), floats=True) 

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

1104 

1105 

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

1107 '''Return the vectorial sum of two or more n-vectors. 

1108 

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

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

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

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

1113 

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

1115 

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

1117 ''' 

1118 return _sumOf(nvectors, Vector=Vector, h=h, **Vector_kwds) 

1119 

1120 

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

1122 height=None, LatLon=LatLon, **LatLon_kwds): 

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

1124 from those points. 

1125 

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

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

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

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

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

1131 the mean height (C{meter}). 

1132 @kwarg LatLon: Optional class to return the triangulated point 

1133 (L{LatLon}). 

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

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

1136 

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

1138 

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

1140 

1141 @raise Valuerror: Points coincide. 

1142 

1143 @example: 

1144 

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

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

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

1148 ''' 

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

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

1151 height=height, LatLon=LatLon, **LatLon_kwds) 

1152 

1153 

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

1155 radius=R_M, height=None, useZ=False, 

1156 LatLon=LatLon, **LatLon_kwds): 

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

1158 

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

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

1161 as B{C{radius}}). 

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

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

1164 as B{C{radius}}). 

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

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

1167 as B{C{radius}}). 

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

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

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

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

1172 @kwarg LatLon: Optional class to return the trilaterated 

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

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

1175 

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

1177 

1178 @raise IntersectionError: No intersection, trilateration failed. 

1179 

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

1181 

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

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

1184 

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

1186 ''' 

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

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

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

1190 radius=radius, height=height, useZ=useZ, 

1191 LatLon=LatLon, **LatLon_kwds) 

1192 

1193 

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

1195 areaOf, # functions 

1196 intersecant2, intersection, ispolar, 

1197 meanOf, 

1198 nearestOn2, nearestOn3, 

1199 perimeterOf, 

1200 sumOf, 

1201 triangulate, trilaterate) 

1202 

1203# **) MIT License 

1204# 

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

1206# 

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

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

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

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

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

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

1213# 

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

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

1216# 

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

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

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

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

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

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

1223# OTHER DEALINGS IN THE SOFTWARE.