Coverage for pygeodesy/sphericalNvector.py: 92%

320 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-28 15: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 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, _r2m 

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.30' 

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}) or C{None}. 

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 (C{radians} 

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

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

152 toward the end point of the line or negative if 

153 "before" the start point. 

154 

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

156 

157 @raise Valuerror: Some points coincide. 

158 

159 @example: 

160 

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

162 

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

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

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

166 ''' 

167 p = self.others(start=start) 

168 n = self.toNvector() 

169 

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

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

172 return _r2m(start.toNvector().angleTo(a, vSign=gc), radius) 

173 

174 @deprecated_method 

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

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

177 ''' 

178 return self.initialBearingTo(other) 

179 

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

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

182 defined by a start and end point. 

183 

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

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

186 initial bearing from start point (compass 

187 C{degrees360}). 

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

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

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

191 

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

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

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

195 of the line . 

196 

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

198 

199 @raise Valuerror: Some points coincide. 

200 

201 @example: 

202 

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

204 

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

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

207 

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

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

210 ''' 

211 p = self.others(start=start) 

212 n = self.toNvector() 

213 

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

215 return _r2m(gc.angleTo(n) - PI_2, radius) 

216 

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

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

219 travelled the given distance on the given bearing. 

220 

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

222 as B{C{radius}}). 

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

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

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

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

227 

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

229 

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

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

232 

233 @example: 

234 

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

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

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

238 ''' 

239 b = Bearing_(bearing) 

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

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

242 

243 n = self.toNvector() 

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

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

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

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

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

249 

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

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

252 

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

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

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

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

257 

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

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

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

261 

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

263 

264 @example: 

265 

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

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

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

269 ''' 

270 p = self.others(other) 

271 if wrap: 

272 p = _unrollon(self, p) 

273 n = p.toNvector() 

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

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

276 

277# @Property_RO 

278# def Ecef(self): 

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

280# ''' 

281# return _ALL_MODS.ecef.EcefKarney 

282 

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

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

285 ''' 

286 s = start.toNvector() 

287 if isscalar(end): # bearing 

288 gc = s.greatCircle(end) 

289 e = None 

290 else: # point 

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

292 if wrap: 

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

294 e = p.toNvector() 

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

296 return gc, s, e 

297 

298 def greatCircle(self, bearing): 

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

300 heading on the given bearing from this point. 

301 

302 Direction of vector is such that initial bearing vector 

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

304 

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

306 

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

308 ''' 

309 t = Bearing_(bearing) 

310 a, b = self.philam 

311 

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

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

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

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

316 

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

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

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

320 

321 Direction of vector is such that initial bearing vector 

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

323 

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

325 this point (compass C{degrees360}). 

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

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

328 

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

330 

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

332 

333 @raise Valuerror: Points coincide. 

334 

335 @example: 

336 

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

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

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

340 

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

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

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

344 ''' 

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

346 return gc.unit() 

347 

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

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

350 to an other point. 

351 

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

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

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

355 

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

357 

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

359 point or the C{NorthPole}, provided 

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

361 

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

363 

364 @example: 

365 

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

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

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

369 ''' 

370 n = self.toNvector() 

371 p = self.others(other) 

372 if wrap: 

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

374 p = p.toNvector() 

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

376# gc1 = self.greatCircleTo(other) 

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

378# gc2 = self.greatCircleTo(NorthPole) 

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

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

381 

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

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

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

385 

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

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

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

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

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

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

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

393 

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

395 

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

397 

398 @example: 

399 

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

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

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

403 ''' 

404 n = self.toNvector() 

405 p = self.others(other) 

406 if wrap: 

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

408 

409 f = Scalar(fraction=fraction) 

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

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

412 

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

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

415 

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

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

418 other point. 

419 

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

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

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

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

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

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

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

427 

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

429 

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

431 

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

433 

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

435 

436 @example: 

437 

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

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

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

441 ''' 

442 n = self.toNvector() 

443 p = self.others(other) 

444 if wrap: 

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

446 p = p.toNvector() 

447 f = Scalar(fraction=fraction) 

448 

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

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

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

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

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

454 

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

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

457 

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

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

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

461 

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

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

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

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

466 initial bearing at the second point (compass 

467 C{degrees}). 

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

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

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

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

472 

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

474 

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

476 point is not L{LatLon}. 

477 

478 @raise ValueError: Intersection is ambiguous or infinite or 

479 the lines are parallel, coincident or null. 

480 

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

482 L{intersection2}. 

483 

484 @example: 

485 

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

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

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

489 ''' 

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

491 wrap=wrap, LatLon=self.classof) 

492 

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

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

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

496 

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

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

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

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

501 initial bearing at the second start point (compass 

502 C{degrees360}). 

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

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

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

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

507 

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

509 

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

511 point is not L{LatLon}. 

512 

513 @raise ValueError: Intersection is ambiguous or infinite or 

514 the lines are parallel, coincident or null. 

515 

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

517 L{intersection}. 

518 ''' 

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

520 wrap=wrap, LatLon=self.classof) 

521 

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

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

524 

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

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

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

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

529 

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

531 C{False} otherwise. 

532 

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

534 

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

536 

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

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

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

540 ''' 

541 if _MODS.booleans.isBoolean(points): 

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

543 

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

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

546 def _subtangles(ps, w): 

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

548 n0 = self.toNvector() 

549 _m0 = n0.minus 

550 p1 = Ps[0] 

551 vs1 = _m0(p1.toNvector()) 

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

553 if w and not Ps.looped: 

554 p2 = _unrollon(p1, p2) 

555 p1 = p2 

556 vs2 = _m0(p2.toNvector()) 

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

558 vs1 = vs2 

559 

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

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

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

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

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

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

566 # XXX are winding number optimisations equally applicable to 

567 # spherical surface? 

568 return fabs(s) > PI 

569 

570 @deprecated_method 

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

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

573 return self.isenclosedBy(points) 

574 

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

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

577 

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

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

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

581 the same hemispere). 

582 

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

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

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

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

587 

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

589 arc, C{False} otherwise. 

590 

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

592 L{LatLon}. 

593 ''' 

594 p1 = self.others(point1=point1) 

595 p2 = self.others(point2=point2) 

596 if wrap: 

597 p1 = _Wrap.point(p1) 

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

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

600 

601 # corner case, null arc 

602 if n1.isequalTo(n2): 

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

604 

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

606 return False # PYCHOK returns 

607 

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

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

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

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

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

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

614 

615 @deprecated_method 

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

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

618 return self.iswithin(point1, point2) 

619 

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

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

622 

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

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

625 the mean height (C{meter}). 

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

627 may be negative or greater than 1.0. 

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

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

630 

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

632 

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

634 

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

636 

637 @example: 

638 

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

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

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

642 ''' 

643 if fraction is _0_5: 

644 p = self.others(other) 

645 if wrap: 

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

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

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

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

650 else: 

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

652 return r 

653 

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

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

656 closest to this point. 

657 

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

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

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

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

662 to interpolate the height. 

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

664 given points, otherwise the closest point 

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

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

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

668 

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

670 

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

672 not supported. 

673 

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

675 

676 @example: 

677 

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

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

680 

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

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

683 

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

685 

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

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

688 ''' 

689 p1 = self.others(point1=point1) 

690 p2 = self.others(point2=point2) 

691 if wrap: 

692 p1 = _Wrap.point(p1) 

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

694 p0 = self 

695 

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

697 # closer to arc than to its endpoints, 

698 # find the closest point on the arc 

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

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

701 n = gc1.cross(gc2) 

702 

703 elif within: # for backward compatibility, XXX unwrapped 

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

705 self.distanceTo(point2)) else point2 

706 

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

708 n1 = p1.toNvector() 

709 n2 = p2.toNvector() 

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

711 if n is n1: 

712 return p1 # is point1 

713 elif n is n2: 

714 return p2 # is point2 if not wrap 

715 

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

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

718 d = p1.distanceTo(p2) 

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

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

721 return p 

722 

723 # @deprecated_method 

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

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

726 

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

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

729 to that point from this point ... 

730 ''' 

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

732 return r.closest, r.distance 

733 

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

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

736 arcs joining consecutive points) closest to this point. 

737 

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

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

740 

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

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

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

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

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

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

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

748 

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

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

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

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

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

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

755 

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

757 

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

759 ''' 

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

761 _r = self.distanceTo 

762 _n = self.nearestOn 

763 

764 c = p1 = Ps[0] 

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

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

767 if wrap and not Ps.looped: 

768 p2 = _unrollon(p1, p2) 

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

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

771 if d < r: 

772 c, r = p, d 

773 p1 = p2 

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

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

776 

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

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

779 

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

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

782 to override this L{Cartesian} class or specify 

783 C{B{Cartesian}=None}. 

784 

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

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

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

788 

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

790 ''' 

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

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

793 

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

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

796 

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

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

799 override this L{Nvector} class or use 

800 C{B{Nvector}=None}. 

801 

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

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

804 

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

806 

807 @example: 

808 

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

810 >>> n = p.toNvector() 

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

812 ''' 

813 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector) 

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

815 

816 

817class Nvector(NvectorBase): 

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

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

820 n-vectors have no singularities or discontinuities. 

821 

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

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

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

825 

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

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

828 

829 Note commonality with L{ellipsoidalNvector.Nvector}. 

830 ''' 

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

832 

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

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

835 (ECEF) coordinates. 

836 

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

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

839 to override this L{Cartesian} class or specify 

840 C{B{Cartesian}=None}. 

841 

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

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

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

845 

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

847 ''' 

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

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

850 

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

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

853 

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

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

856 to override this L{LatLon} class or specify 

857 C{B{LatLon}=None}. 

858 

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

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

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

862 

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

864 

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

866 ''' 

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

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

869 

870 def greatCircle(self, bearing): 

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

872 heading on given compass bearing from this point as its 

873 n-vector. 

874 

875 Direction of vector is such that initial bearing vector 

876 b = c × p. 

877 

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

879 

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

881 

882 @raise Valuerror: Polar coincidence. 

883 

884 @example: 

885 

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

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

888 ''' 

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

890 

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

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

893 

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

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

896 return n.minus(e) 

897 

898 

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

900 

901 

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

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

904 great circle arcs joining consecutive points). 

905 

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

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

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

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

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

911 

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

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

914 

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

916 

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

918 

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

920 and L{ellipsoidalKarney.areaOf}. 

921 

922 @example: 

923 

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

925 >>> areaOf(b) # 8666058750.718977 

926 ''' 

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

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

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

930 n0 = Ps[0].toNvector() 

931 

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

933 p1 = Ps[1] 

934 v1 = p1._N_vector 

935 gc = v2.cross(v1) 

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

937 if w and not Ps.looped: 

938 p2 = _unrollon(p1, p2) 

939 p1 = p2 

940 v2 = p2._N_vector 

941 gc1 = v1.cross(v2) 

942 v1 = v2 

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

944 gc = gc1 

945 

946 if _MODS.booleans.isBoolean(points): 

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

948 else: 

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

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

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

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

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

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

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

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

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

958 

959 

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

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

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

963 

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

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

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

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

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

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

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

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

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

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

974 circle} methods. 

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

976 conventionally) or C{None}. 

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

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

979 

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

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

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

983 

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

985 

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

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

988 

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

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

991 ''' 

992 c = _Nvll.others(center=center) 

993 p = _Nvll.others(point=point) 

994 try: 

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

996 height=height, wrap=wrap) 

997 except (TypeError, ValueError) as x: 

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

999 exact=exact, wrap=wrap) 

1000 

1001 

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

1003 LatLon=LatLon, **LatLon_kwds): 

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

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

1006 

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

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

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

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

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

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

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

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

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

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

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

1018 (L{LatLon}). 

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

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

1021 

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

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

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

1025 

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

1027 

1028 @raise ValueError: Intersection is ambiguous or infinite or 

1029 the lines are parallel, coincident or null. 

1030 

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

1032 

1033 @example: 

1034 

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

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

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

1038 ''' 

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

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

1041 return i.toLatLon(**kwds) 

1042 

1043 

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

1045 LatLon=LatLon, **LatLon_kwds): 

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

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

1048 

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

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

1051 initial bearing at the first start point 

1052 (compass C{degrees360}). 

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

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

1055 initial bearing at the second start point 

1056 (compass C{degrees360}). 

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

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

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

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

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

1062 antipodal points (L{LatLon}). 

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

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

1065 

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

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

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

1069 if available. 

1070 

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

1072 

1073 @raise ValueError: Intersection is ambiguous or infinite or 

1074 the lines are parallel, coincident or null. 

1075 

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

1077 ''' 

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

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

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

1081 

1082 

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

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

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

1086 ''' 

1087 p1 = _Nvll.others(start1=start1) 

1088 p2 = _Nvll.others(start2=start2) 

1089 if wrap: 

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

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

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

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

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

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

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

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

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

1099 

1100 hs = start1.height, start2.height 

1101 # there are two (antipodal) candidate intersection 

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

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

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

1105 

1106 # selection of intersection point depends on how 

1107 # lines are defined (by bearings or endpoints) 

1108 if e1 and e2: # endpoint+endpoint 

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

1110 hs += end1.height, end2.height 

1111 elif e1 and not e2: # endpoint+bearing 

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

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

1114 hs += end1.height, 

1115 elif e2 and not e1: # bearing+endpoint 

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

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

1118 hs += end2.height, 

1119 else: # bearing+bearing 

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

1121 # towards i1, otherwise towards antipodal i2 

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

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

1124 if d1 > 0 and d2 > 0: 

1125 d = 1 # both point to i1 

1126 elif d1 < 0 and d2 < 0: 

1127 d = -1 # both point to i2 

1128 else: # d1, d2 opposite signs 

1129 # intersection is at further-away intersection point, 

1130 # take opposite intersection from mid- point of v1 

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

1132 # get intersection p1 bearing points to, aka being 

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

1134 # function .sphericalTrigonometry._intersect and 

1135 # .ellipsoidalBaseDI._intersect3 

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

1137 

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

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

1140 

1141 

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

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

1144 

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

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

1147 (C{meter}). 

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

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

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

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

1152 

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

1154 

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

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

1157 ''' 

1158 def _N_vs(ps, w): 

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

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

1161 yield p._N_vector 

1162 

1163 try: 

1164 # geographic mean 

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

1166 except (TypeError, ValueError) as x: 

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

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

1169 name=meanOf.__name__)) 

1170 

1171 

1172@deprecated_function 

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

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

1175 

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

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

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

1179 ''' 

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

1181 return r.closest, r.distance 

1182 

1183 

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

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

1186 consecutive points) closest to an other point. 

1187 

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

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

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

1191 

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

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

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

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

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

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

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

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

1200 

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

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

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

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

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

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

1207 C{degrees360}. 

1208 

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

1210 

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

1212 ''' 

1213 _xinstanceof(LatLon, point=point) 

1214 

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

1216 height=height, wrap=wrap) 

1217 

1218 

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

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

1221 great circle arcs joining consecutive points). 

1222 

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

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

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

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

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

1228 

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

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

1231 

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

1233 

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

1235 

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

1237 C{B{points}} a composite. 

1238 

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

1240 and L{sphericalTrigonometry.perimeterOf}. 

1241 ''' 

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

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

1244 p1 = Ps[0] 

1245 v1 = p1._N_vector 

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

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

1248 p2 = _unrollon(p1, p2) 

1249 p1 = p2 

1250 v2 = p2._N_vector 

1251 yield v1.angleTo(v2) 

1252 v1 = v2 

1253 

1254 if _MODS.booleans.isBoolean(points): 

1255 if not closed: 

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

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

1258 else: 

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

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

1261 

1262 

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

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

1265 

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

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

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

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

1270 

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

1272 

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

1274 ''' 

1275 try: 

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

1277 except (TypeError, ValueError) as x: 

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

1279 

1280 

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

1282 height=None, wrap=False, 

1283 LatLon=LatLon, **LatLon_kwds): 

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

1285 from those points. 

1286 

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

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

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

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

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

1292 the mean height (C{meter}). 

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

1294 (C{bool}). 

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

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

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

1298 

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

1300 

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

1302 

1303 @raise Valuerror: Points coincide. 

1304 

1305 @example: 

1306 

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

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

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

1310 ''' 

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

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

1313 height=height, wrap=wrap, 

1314 LatLon=LatLon, **LatLon_kwds) 

1315 

1316 

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

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

1319 LatLon=LatLon, **LatLon_kwds): 

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

1321 

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

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

1324 as B{C{radius}}). 

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

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

1327 as B{C{radius}}). 

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

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

1330 as B{C{radius}}). 

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

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

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

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

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

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

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

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

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

1340 

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

1342 

1343 @raise IntersectionError: No intersection, trilateration failed. 

1344 

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

1346 

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

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

1349 

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

1351 ''' 

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

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

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

1355 radius=radius, height=height, useZ=useZ, 

1356 wrap=wrap, LatLon=LatLon, **LatLon_kwds) 

1357 

1358 

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

1360 areaOf, # functions 

1361 intersecant2, intersection, intersection2, ispolar, 

1362 meanOf, 

1363 nearestOn2, nearestOn3, 

1364 perimeterOf, 

1365 sumOf, 

1366 triangulate, trilaterate) 

1367 

1368# **) MIT License 

1369# 

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

1371# 

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

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

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

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

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

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

1378# 

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

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

1381# 

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

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

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

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

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

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

1388# OTHER DEALINGS IN THE SOFTWARE.