Coverage for pygeodesy/sphericalNvector.py: 92%

320 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-07-12 13:40 -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 lines 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 PointsError, VectorError, _xError, _xkwds 

41from pygeodesy.fmath import fmean, fsum 

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

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

44 _point_, _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 LatLonNvectorBase, NorthPole, \ 

49 notImplemented, NvectorBase, _nsumOf, \ 

50 _triangulate, _trilaterate 

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

52from pygeodesy.props import deprecated_function, deprecated_method 

53from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, \ 

54 _intersecant2, LatLonSphericalBase, \ 

55 Datums 

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

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

58 sincos2d, _unrollon, _Wrap 

59 

60# from math import atan2, fabs # from utily 

61 

62__all__ = _ALL_LAZY.sphericalNvector 

63__version__ = '23.06.12' 

64 

65_lines_ = 'lines' 

66 

67 

68class Cartesian(CartesianSphericalBase): 

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

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

71 ''' 

72 

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

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

75 

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

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

78 to override this L{LatLon} class or specify 

79 C{B{LatLon}=None}. 

80 

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

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

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

84 

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

86 ''' 

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

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

89 

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

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

92 

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

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

95 to override this L{Nvector} class or specify 

96 C{B{Nvector}=None}. 

97 

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

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

100 

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

102 ''' 

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

104 # datum=datum or self.datum) 

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

106 # return ll.toNvector(**kwds) 

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

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

109 

110 

111class LatLon(LatLonNvectorBase, LatLonSphericalBase): 

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

113 

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

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

116 

117 @example: 

118 

119 >>> from sphericalNvector import LatLon 

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

121 ''' 

122 _Nv = None # cached_toNvector L{Nvector}) 

123 

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

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

126 ''' 

127 if updated: # reset caches 

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

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

130 

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

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

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

134 end point. 

135 

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

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

138 from the start point to the point where the perpendicular 

139 crosses the line. 

140 

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

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

143 initial bearing from start point (compass 

144 C{degrees360}). 

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

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

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

148 

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

150 after the start toward the end point of the line 

151 or negative if before the start point). 

152 

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

154 

155 @raise Valuerror: Some points coincide. 

156 

157 @example: 

158 

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

160 

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

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

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

164 ''' 

165 p = self.others(start=start) 

166 n = self.toNvector() 

167 

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

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

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

171 

172 @deprecated_method 

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

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

175 ''' 

176 return self.initialBearingTo(other) 

177 

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

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

180 defined by a start and end point. 

181 

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

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

184 initial bearing from start point (compass 

185 C{degrees360}). 

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

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

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

189 

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

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

192 

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

194 

195 @raise Valuerror: Some points coincide. 

196 

197 @example: 

198 

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

200 

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

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

203 

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

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

206 ''' 

207 p = self.others(start=start) 

208 n = self.toNvector() 

209 

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

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

212 

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

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

215 travelled the given distance on the given bearing. 

216 

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

218 as B{C{radius}}). 

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

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

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

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

223 

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

225 

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

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

228 

229 @example: 

230 

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

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

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

234 ''' 

235 b = Bearing_(bearing) 

236 a = _angular(distance, radius, low=None) 

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

238 

239 n = self.toNvector() 

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

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

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

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

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

245 

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

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

248 

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

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

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

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

253 

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

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

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

257 

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

259 

260 @example: 

261 

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

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

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

265 ''' 

266 p = self.others(other) 

267 if wrap: 

268 p = _unrollon(self, p) 

269 n = p.toNvector() 

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

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

272 

273# @Property_RO 

274# def Ecef(self): 

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

276# ''' 

277# return _ALL_MODS.ecef.EcefKarney 

278 

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

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

281 ''' 

282 s = start.toNvector() 

283 if isscalar(end): # bearing 

284 gc = s.greatCircle(end) 

285 e = None 

286 else: # point 

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

288 if wrap: 

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

290 e = p.toNvector() 

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

292 return gc, s, e 

293 

294 def greatCircle(self, bearing): 

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

296 heading on the given bearing from this point. 

297 

298 Direction of vector is such that initial bearing vector 

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

300 

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

302 

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

304 ''' 

305 t = Bearing_(bearing) 

306 a, b = self.philam 

307 

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

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

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

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

312 

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

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

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

316 

317 Direction of vector is such that initial bearing vector 

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

319 

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

321 this point (compass C{degrees360}). 

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

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

324 

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

326 

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

328 

329 @raise Valuerror: Points coincide. 

330 

331 @example: 

332 

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

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

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

336 

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

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

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

340 ''' 

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

342 return gc.unit() 

343 

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

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

346 to an other point. 

347 

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

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

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

351 

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

353 

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

355 point or the C{NorthPole}, provided 

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

357 

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

359 

360 @example: 

361 

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

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

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

365 ''' 

366 n = self.toNvector() 

367 p = self.others(other) 

368 if wrap: 

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

370 p = p.toNvector() 

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

372# gc1 = self.greatCircleTo(other) 

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

374# gc2 = self.greatCircleTo(NorthPole) 

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

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

377 

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

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

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

381 

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

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

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

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

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

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

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

389 

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

391 

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

393 

394 @example: 

395 

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

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

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

399 ''' 

400 n = self.toNvector() 

401 p = self.others(other) 

402 if wrap: 

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

404 

405 f = Scalar(fraction=fraction) 

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

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

408 

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

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

411 

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

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

414 other point. 

415 

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

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

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

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

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

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

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

423 

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

425 

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

427 

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

429 

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

431 

432 @example: 

433 

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

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

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

437 ''' 

438 n = self.toNvector() 

439 p = self.others(other) 

440 if wrap: 

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

442 p = p.toNvector() 

443 f = Scalar(fraction=fraction) 

444 

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

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

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

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

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

450 

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

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

453 

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

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

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

457 

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

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

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

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

462 initial bearing at the second point (compass 

463 C{degrees}). 

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

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

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

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

468 

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

470 

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

472 point is not L{LatLon}. 

473 

474 @raise ValueError: Intersection is ambiguous or infinite or 

475 the lines are parallel, coincident or null. 

476 

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

478 L{intersection2}. 

479 

480 @example: 

481 

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

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

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

485 ''' 

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

487 wrap=wrap, LatLon=self.classof) 

488 

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

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

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

492 

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

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

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

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

497 initial bearing at the second start point (compass 

498 C{degrees360}). 

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

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

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

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

503 

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

505 

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

507 point is not L{LatLon}. 

508 

509 @raise ValueError: Intersection is ambiguous or infinite or 

510 the lines are parallel, coincident or null. 

511 

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

513 L{intersection}. 

514 ''' 

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

516 wrap=wrap, LatLon=self.classof) 

517 

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

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

520 

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

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

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

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

525 

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

527 C{False} otherwise. 

528 

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

530 

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

532 

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

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

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

536 ''' 

537 if _MODS.booleans.isBoolean(points): 

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

539 

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

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

542 def _subtangles(ps, w): 

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

544 n0 = self.toNvector() 

545 _m0 = n0.minus 

546 p1 = Ps[0] 

547 vs1 = _m0(p1.toNvector()) 

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

549 if w and not Ps.looped: 

550 p2 = _unrollon(p1, p2) 

551 p1 = p2 

552 vs2 = _m0(p2.toNvector()) 

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

554 vs1 = vs2 

555 

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

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

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

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

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

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

562 # XXX are winding number optimisations equally applicable to 

563 # spherical surface? 

564 return fabs(s) > PI 

565 

566 @deprecated_method 

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

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

569 return self.isenclosedBy(points) 

570 

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

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

573 

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

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

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

577 the same hemispere). 

578 

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

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

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

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

583 

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

585 arc, C{False} otherwise. 

586 

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

588 L{LatLon}. 

589 ''' 

590 p1 = self.others(point1=point1) 

591 p2 = self.others(point2=point2) 

592 if wrap: 

593 p1 = _Wrap.point(p1) 

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

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

596 

597 # corner case, null arc 

598 if n1.isequalTo(n2): 

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

600 

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

602 return False # PYCHOK returns 

603 

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

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

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

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

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

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

610 

611 @deprecated_method 

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

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

614 return self.iswithin(point1, point2) 

615 

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

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

618 

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

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

621 the mean height (C{meter}). 

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

623 may be negative or greater than 1.0. 

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

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

626 

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

628 

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

630 

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

632 

633 @example: 

634 

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

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

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

638 ''' 

639 if fraction is _0_5: 

640 p = self.others(other) 

641 if wrap: 

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

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

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

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

646 else: 

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

648 return r 

649 

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

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

652 closest to this point. 

653 

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

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

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

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

658 to interpolate the height. 

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

660 given points, otherwise the closest point 

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

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

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

664 

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

666 

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

668 not supported. 

669 

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

671 

672 @example: 

673 

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

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

676 

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

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

679 

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

681 

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

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

684 ''' 

685 p1 = self.others(point1=point1) 

686 p2 = self.others(point2=point2) 

687 if wrap: 

688 p1 = _Wrap.point(p1) 

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

690 p0 = self 

691 

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

693 # closer to arc than to its endpoints, 

694 # find the closest point on the arc 

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

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

697 n = gc1.cross(gc2) 

698 

699 elif within: # for backward compatibility, XXX unwrapped 

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

701 self.distanceTo(point2)) else point2 

702 

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

704 n1 = p1.toNvector() 

705 n2 = p2.toNvector() 

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

707 if n is n1: 

708 return p1 # is point1 

709 elif n is n2: 

710 return p2 # is point2 if not wrap 

711 

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

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

714 d = p1.distanceTo(p2) 

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

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

717 return p 

718 

719 # @deprecated_method 

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

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

722 

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

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

725 to that point from this point ... 

726 ''' 

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

728 return r.closest, r.distance 

729 

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

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

732 arcs joining consecutive points) closest to this point. 

733 

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

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

736 

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

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

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

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

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

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

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

744 

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

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

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

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

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

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

751 

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

753 

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

755 ''' 

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

757 _r = self.distanceTo 

758 _n = self.nearestOn 

759 

760 c = p1 = Ps[0] 

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

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

763 if wrap and not Ps.looped: 

764 p2 = _unrollon(p1, p2) 

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

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

767 if d < r: 

768 c, r = p, d 

769 p1 = p2 

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

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

772 

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

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

775 

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

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

778 to override this L{Cartesian} class or specify 

779 C{B{Cartesian}=None}. 

780 

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

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

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

784 

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

786 ''' 

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

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

789 

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

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

792 

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

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

795 override this L{Nvector} class or use 

796 C{B{Nvector}=None}. 

797 

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

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

800 

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

802 

803 @example: 

804 

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

806 >>> n = p.toNvector() 

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

808 ''' 

809 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector) 

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

811 

812 

813class Nvector(NvectorBase): 

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

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

816 n-vectors have no singularities or discontinuities. 

817 

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

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

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

821 

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

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

824 

825 Note commonality with L{ellipsoidalNvector.Nvector}. 

826 ''' 

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

828 

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

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

831 (ECEF) coordinates. 

832 

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

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

835 to override this L{Cartesian} class or specify 

836 C{B{Cartesian}=None}. 

837 

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

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

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

841 

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

843 ''' 

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

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

846 

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

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

849 

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

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

852 to override this L{LatLon} class or specify 

853 C{B{LatLon}=None}. 

854 

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

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

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

858 

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

860 

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

862 ''' 

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

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

865 

866 def greatCircle(self, bearing): 

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

868 heading on given compass bearing from this point as its 

869 n-vector. 

870 

871 Direction of vector is such that initial bearing vector 

872 b = c × p. 

873 

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

875 

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

877 

878 @raise Valuerror: Polar coincidence. 

879 

880 @example: 

881 

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

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

884 ''' 

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

886 

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

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

889 

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

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

892 return n.minus(e) 

893 

894 

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

896 

897 

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

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

900 great circle arcs joining consecutive points). 

901 

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

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

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

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

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

907 

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

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

910 

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

912 

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

914 

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

916 and L{ellipsoidalKarney.areaOf}. 

917 

918 @example: 

919 

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

921 >>> areaOf(b) # 8666058750.718977 

922 ''' 

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

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

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

926 n0 = Ps[0].toNvector() 

927 

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

929 p1 = Ps[1] 

930 v1 = p1._N_vector 

931 gc = v2.cross(v1) 

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

933 if w and not Ps.looped: 

934 p2 = _unrollon(p1, p2) 

935 p1 = p2 

936 v2 = p2._N_vector 

937 gc1 = v1.cross(v2) 

938 v1 = v2 

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

940 gc = gc1 

941 

942 if _MODS.booleans.isBoolean(points): 

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

944 else: 

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

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

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

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

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

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

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

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

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

954 

955 

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

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

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

959 

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

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

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

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

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

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

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

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

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

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

970 circle} methods. 

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

972 conventionally) or C{None}. 

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

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

975 

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

977 an instance of this class. For a tangent line, both points are 

978 the same instance, the B{C{point}} or wrapped or I{normalized}. 

979 

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

981 

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

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

984 

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

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

987 ''' 

988 c = _Nvll.others(center=center) 

989 p = _Nvll.others(point=point) 

990 try: 

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

992 height=height, wrap=wrap) 

993 except (TypeError, ValueError) as x: 

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

995 exact=exact, wrap=wrap) 

996 

997 

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

999 LatLon=LatLon, **LatLon_kwds): 

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

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

1002 

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

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

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

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

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

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

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

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

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

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

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

1014 (L{LatLon}). 

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

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

1017 

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

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

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

1021 

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

1023 

1024 @raise ValueError: Intersection is ambiguous or infinite or 

1025 the lines are parallel, coincident or null. 

1026 

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

1028 

1029 @example: 

1030 

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

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

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

1034 ''' 

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

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

1037 return i.toLatLon(**kwds) 

1038 

1039 

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

1041 LatLon=LatLon, **LatLon_kwds): 

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

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

1044 

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

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

1047 initial bearing at the first start point 

1048 (compass C{degrees360}). 

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

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

1051 initial bearing at the second start point 

1052 (compass C{degrees360}). 

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

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

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

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

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

1058 antipodal points (L{LatLon}). 

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

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

1061 

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

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

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

1065 if available. 

1066 

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

1068 

1069 @raise ValueError: Intersection is ambiguous or infinite or 

1070 the lines are parallel, coincident or null. 

1071 

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

1073 ''' 

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

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

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

1077 

1078 

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

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

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

1082 ''' 

1083 p1 = _Nvll.others(start1=start1) 

1084 p2 = _Nvll.others(start2=start2) 

1085 if wrap: 

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

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

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

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

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

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

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

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

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

1095 

1096 hs = start1.height, start2.height 

1097 # there are two (antipodal) candidate intersection 

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

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

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

1101 

1102 # selection of intersection point depends on how 

1103 # lines are defined (by bearings or endpoints) 

1104 if e1 and e2: # endpoint+endpoint 

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

1106 hs += end1.height, end2.height 

1107 elif e1 and not e2: # endpoint+bearing 

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

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

1110 hs += end1.height, 

1111 elif e2 and not e1: # bearing+endpoint 

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

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

1114 hs += end2.height, 

1115 else: # bearing+bearing 

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

1117 # towards i1, otherwise towards antipodal i2 

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

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

1120 if d1 > 0 and d2 > 0: 

1121 d = 1 # both point to i1 

1122 elif d1 < 0 and d2 < 0: 

1123 d = -1 # both point to i2 

1124 else: # d1, d2 opposite signs 

1125 # intersection is at further-away intersection point, 

1126 # take opposite intersection from mid- point of v1 

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

1128 # get intersection p1 bearing points to, aka being 

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

1130 # function .sphericalTrigonometry._intersect and 

1131 # .ellipsoidalBaseDI._intersect3 

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

1133 

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

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

1136 

1137 

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

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

1140 

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

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

1143 (C{meter}). 

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

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

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

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

1148 

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

1150 

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

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

1153 ''' 

1154 def _N_vs(ps, w): 

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

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

1157 yield p._N_vector 

1158 

1159 try: 

1160 # geographic mean 

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

1162 except (TypeError, ValueError) as x: 

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

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

1165 name=meanOf.__name__)) 

1166 

1167 

1168@deprecated_function 

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

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

1171 

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

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

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

1175 ''' 

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

1177 return r.closest, r.distance 

1178 

1179 

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

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

1182 consecutive points) closest to an other point. 

1183 

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

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

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

1187 

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

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

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

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

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

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

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

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

1196 

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

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

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

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

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

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

1203 C{degrees360}. 

1204 

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

1206 

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

1208 ''' 

1209 _xinstanceof(LatLon, point=point) 

1210 

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

1212 height=height, wrap=wrap) 

1213 

1214 

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

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

1217 great circle arcs joining consecutive points). 

1218 

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

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

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

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

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

1224 

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

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

1227 

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

1229 

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

1231 

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

1233 C{B{points}} a composite. 

1234 

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

1236 and L{sphericalTrigonometry.perimeterOf}. 

1237 ''' 

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

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

1240 p1 = Ps[0] 

1241 v1 = p1._N_vector 

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

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

1244 p2 = _unrollon(p1, p2) 

1245 p1 = p2 

1246 v2 = p2._N_vector 

1247 yield v1.angleTo(v2) 

1248 v1 = v2 

1249 

1250 if _MODS.booleans.isBoolean(points): 

1251 if not closed: 

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

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

1254 else: 

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

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

1257 

1258 

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

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

1261 

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

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

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

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

1266 

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

1268 

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

1270 ''' 

1271 try: 

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

1273 except (TypeError, ValueError) as x: 

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

1275 

1276 

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

1278 height=None, wrap=False, 

1279 LatLon=LatLon, **LatLon_kwds): 

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

1281 from those points. 

1282 

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

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

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

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

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

1288 the mean height (C{meter}). 

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

1290 (C{bool}). 

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

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

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

1294 

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

1296 

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

1298 

1299 @raise Valuerror: Points coincide. 

1300 

1301 @example: 

1302 

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

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

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

1306 ''' 

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

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

1309 height=height, wrap=wrap, 

1310 LatLon=LatLon, **LatLon_kwds) 

1311 

1312 

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

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

1315 LatLon=LatLon, **LatLon_kwds): 

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

1317 

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

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

1320 as B{C{radius}}). 

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

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

1323 as B{C{radius}}). 

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

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

1326 as B{C{radius}}). 

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

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

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

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

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

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

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

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

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

1336 

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

1338 

1339 @raise IntersectionError: No intersection, trilateration failed. 

1340 

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

1342 

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

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

1345 

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

1347 ''' 

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

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

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

1351 radius=radius, height=height, useZ=useZ, 

1352 wrap=wrap, LatLon=LatLon, **LatLon_kwds) 

1353 

1354 

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

1356 areaOf, # functions 

1357 intersecant2, intersection, intersection2, ispolar, 

1358 meanOf, 

1359 nearestOn2, nearestOn3, 

1360 perimeterOf, 

1361 sumOf, 

1362 triangulate, trilaterate) 

1363 

1364# **) MIT License 

1365# 

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

1367# 

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

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

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

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

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

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

1374# 

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

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

1377# 

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

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

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

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

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

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

1384# OTHER DEALINGS IN THE SOFTWARE.