Coverage for pygeodesy/sphericalNvector.py: 92%

323 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-12-02 13:46 -0500

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 C{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, \ 

53 property_RO 

54from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \ 

55 _intersecant2, LatLonSphericalBase, \ 

56 _radians2m, Datums 

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

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

59 sincos2d, _unrollon, _Wrap 

60 

61# from math import atan2, fabs # from utily 

62 

63__all__ = _ALL_LAZY.sphericalNvector 

64__version__ = '23.11.27' 

65 

66_lines_ = 'lines' 

67 

68 

69class Cartesian(CartesianSphericalBase): 

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

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

72 ''' 

73 

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

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

76 

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

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

79 to override this L{LatLon} class or specify 

80 C{B{LatLon}=None}. 

81 

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

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

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

85 

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

87 ''' 

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

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

90 

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

92 '''Convert this cartesian to C{Nvector} components, I{including height}. 

93 

94 @kwarg Nvector_and_kwds: Optional C{Nvector} and C{Nvector} keyword 

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

96 to override this C{Nvector} class or specify 

97 C{B{Nvector}=None}. 

98 

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

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

101 

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

103 ''' 

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

105 # datum=datum or self.datum) 

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

107 # return ll.toNvector(**kwds) 

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

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

110 

111 

112class LatLon(LatLonNvectorBase, LatLonSphericalBase): 

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

114 

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

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

117 

118 @example: 

119 

120 >>> from sphericalNvector import LatLon 

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

122 ''' 

123 _Nv = None # cached_toNvector C{Nvector}) 

124 

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

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

127 ''' 

128 if updated: # reset caches 

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

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

131 

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

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

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

135 end point. 

136 

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

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

139 from the start point to the point where the perpendicular 

140 crosses the line. 

141 

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

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

144 initial bearing from start point (compass 

145 C{degrees360}). 

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

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

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

149 

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

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

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

153 toward the end point of the line or negative if 

154 "before" the start point. 

155 

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

157 

158 @raise Valuerror: Some points coincide. 

159 

160 @example: 

161 

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

163 

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

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

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

167 ''' 

168 p = self.others(start=start) 

169 n = self.toNvector() 

170 

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

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

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

174 

175 @deprecated_method 

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

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

178 ''' 

179 return self.initialBearingTo(other) 

180 

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

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

183 defined by a start and end point. 

184 

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

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

187 initial bearing from start point (compass 

188 C{degrees360}). 

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

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

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

192 

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

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

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

196 of the line . 

197 

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

199 

200 @raise Valuerror: Some points coincide. 

201 

202 @example: 

203 

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

205 

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

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

208 

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

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

211 ''' 

212 p = self.others(start=start) 

213 n = self.toNvector() 

214 

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

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

217 

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

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

220 travelled the given distance on the given bearing. 

221 

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

223 as B{C{radius}}). 

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

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

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

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

228 

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

230 

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

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

233 

234 @example: 

235 

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

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

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

239 ''' 

240 b = Bearing_(bearing) 

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

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

243 

244 n = self.toNvector() 

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

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

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

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

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

250 

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

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

253 

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

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

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

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

258 

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

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

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

262 

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

264 

265 @example: 

266 

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

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

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

270 ''' 

271 p = self.others(other) 

272 if wrap: 

273 p = _unrollon(self, p) 

274 n = p.toNvector() 

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

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

277 

278# @Property_RO 

279# def Ecef(self): 

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

281# ''' 

282# return _ALL_MODS.ecef.EcefKarney 

283 

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

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

286 ''' 

287 s = start.toNvector() 

288 if isscalar(end): # bearing 

289 gc = s.greatCircle(end) 

290 e = None 

291 else: # point 

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

293 if wrap: 

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

295 e = p.toNvector() 

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

297 return gc, s, e 

298 

299 def greatCircle(self, bearing): 

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

301 heading on the given bearing from this point. 

302 

303 Direction of vector is such that initial bearing vector 

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

305 

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

307 

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

309 ''' 

310 t = Bearing_(bearing) 

311 a, b = self.philam 

312 

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

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

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

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

317 

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

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

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

321 

322 Direction of vector is such that initial bearing vector 

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

324 

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

326 this point (compass C{degrees360}). 

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

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

329 

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

331 

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

333 

334 @raise Valuerror: Points coincide. 

335 

336 @example: 

337 

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

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

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

341 

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

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

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

345 ''' 

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

347 return gc.unit() 

348 

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

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

351 to an other point. 

352 

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

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

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

356 

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

358 

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

360 point or the C{NorthPole}, provided 

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

362 

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

364 

365 @example: 

366 

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

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

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

370 ''' 

371 n = self.toNvector() 

372 p = self.others(other) 

373 if wrap: 

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

375 p = p.toNvector() 

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

377# gc1 = self.greatCircleTo(other) 

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

379# gc2 = self.greatCircleTo(NorthPole) 

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

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

382 

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

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

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

386 

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

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

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

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

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

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

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

394 

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

396 

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

398 

399 @example: 

400 

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

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

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

404 ''' 

405 n = self.toNvector() 

406 p = self.others(other) 

407 if wrap: 

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

409 

410 f = Scalar(fraction=fraction) 

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

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

413 

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

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

416 

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

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

419 other point. 

420 

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

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

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

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

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

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

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

428 

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

430 

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

432 

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

434 

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

436 

437 @example: 

438 

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

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

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

442 ''' 

443 n = self.toNvector() 

444 p = self.others(other) 

445 if wrap: 

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

447 p = p.toNvector() 

448 f = Scalar(fraction=fraction) 

449 

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

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

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

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

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

455 

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

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

458 

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

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

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

462 

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

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

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

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

467 initial bearing at the second point (compass 

468 C{degrees}). 

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

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

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

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

473 

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

475 

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

477 point is not L{LatLon}. 

478 

479 @raise ValueError: Intersection is ambiguous or infinite or 

480 the lines are parallel, coincident or null. 

481 

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

483 L{intersection2}. 

484 

485 @example: 

486 

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

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

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

490 ''' 

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

492 wrap=wrap, LatLon=self.classof) 

493 

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

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

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

497 

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

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

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

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

502 initial bearing at the second start point (compass 

503 C{degrees360}). 

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

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

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

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

508 

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

510 

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

512 point is not L{LatLon}. 

513 

514 @raise ValueError: Intersection is ambiguous or infinite or 

515 the lines are parallel, coincident or null. 

516 

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

518 L{intersection}. 

519 ''' 

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

521 wrap=wrap, LatLon=self.classof) 

522 

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

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

525 

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

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

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

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

530 

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

532 C{False} otherwise. 

533 

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

535 

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

537 

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

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

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

541 ''' 

542 if _MODS.booleans.isBoolean(points): 

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

544 

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

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

547 def _subtangles(ps, w): 

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

549 n0 = self.toNvector() 

550 _m0 = n0.minus 

551 p1 = Ps[0] 

552 vs1 = _m0(p1.toNvector()) 

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

554 if w and not Ps.looped: 

555 p2 = _unrollon(p1, p2) 

556 p1 = p2 

557 vs2 = _m0(p2.toNvector()) 

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

559 vs1 = vs2 

560 

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

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

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

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

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

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

567 # XXX are winding number optimisations equally applicable to 

568 # spherical surface? 

569 return fabs(s) > PI 

570 

571 @deprecated_method 

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

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

574 return self.isenclosedBy(points) 

575 

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

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

578 

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

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

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

582 the same hemispere). 

583 

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

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

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

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

588 

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

590 arc, C{False} otherwise. 

591 

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

593 L{LatLon}. 

594 ''' 

595 p1 = self.others(point1=point1) 

596 p2 = self.others(point2=point2) 

597 if wrap: 

598 p1 = _Wrap.point(p1) 

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

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

601 

602 # corner case, null arc 

603 if n1.isequalTo(n2): 

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

605 

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

607 return False # PYCHOK returns 

608 

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

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

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

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

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

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

615 

616 @deprecated_method 

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

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

619 return self.iswithin(point1, point2) 

620 

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

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

623 

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

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

626 the mean height (C{meter}). 

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

628 may be negative or greater than 1.0. 

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

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

631 

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

633 

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

635 

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

637 

638 @example: 

639 

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

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

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

643 ''' 

644 if fraction is _0_5: 

645 p = self.others(other) 

646 if wrap: 

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

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

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

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

651 else: 

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

653 return r 

654 

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

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

657 closest to this point. 

658 

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

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

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

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

663 to interpolate the height. 

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

665 given points, otherwise the closest point 

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

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

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

669 

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

671 

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

673 not supported. 

674 

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

676 

677 @example: 

678 

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

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

681 

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

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

684 

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

686 

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

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

689 ''' 

690 p1 = self.others(point1=point1) 

691 p2 = self.others(point2=point2) 

692 if wrap: 

693 p1 = _Wrap.point(p1) 

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

695 p0 = self 

696 

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

698 # closer to arc than to its endpoints, 

699 # find the closest point on the arc 

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

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

702 n = gc1.cross(gc2) 

703 

704 elif within: # for backward compatibility, XXX unwrapped 

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

706 self.distanceTo(point2)) else point2 

707 

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

709 n1 = p1.toNvector() 

710 n2 = p2.toNvector() 

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

712 if n is n1: 

713 return p1 # is point1 

714 elif n is n2: 

715 return p2 # is point2 if not wrap 

716 

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

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

719 d = p1.distanceTo(p2) 

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

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

722 return p 

723 

724 # @deprecated_method 

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

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

727 

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

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

730 to that point from this point ... 

731 ''' 

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

733 return r.closest, r.distance 

734 

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

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

737 arcs joining consecutive points) closest to this point. 

738 

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

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

741 

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

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

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

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

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

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

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

749 

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

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

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

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

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

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

756 

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

758 

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

760 ''' 

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

762 _r = self.distanceTo 

763 _n = self.nearestOn 

764 

765 c = p1 = Ps[0] 

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

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

768 if wrap and not Ps.looped: 

769 p2 = _unrollon(p1, p2) 

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

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

772 if d < r: 

773 c, r = p, d 

774 p1 = p2 

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

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

777 

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

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

780 

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

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

783 to override this L{Cartesian} class or specify 

784 C{B{Cartesian}=None}. 

785 

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

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

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

789 

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

791 ''' 

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

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

794 

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

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

797 

798 @kwarg Nvector_and_kwds: Optional C{Nvector} and C{Nvector} keyword 

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

800 override this C{Nvector} class or use 

801 C{B{Nvector}=None}. 

802 

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

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

805 

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

807 

808 @example: 

809 

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

811 >>> n = p.toNvector() 

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

813 ''' 

814 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector) 

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

816 

817 

818class Nvector(NvectorBase): 

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

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

821 n-vectors have no singularities or discontinuities. 

822 

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

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

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

826 

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

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

829 

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

831 ''' 

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

833 

834 @property_RO 

835 def sphericalNvector(self): 

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

837 ''' 

838 return type(self) 

839 

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

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

842 (ECEF) coordinates. 

843 

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

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

846 to override this L{Cartesian} class or specify 

847 C{B{Cartesian}=None}. 

848 

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

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

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

852 

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

854 ''' 

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

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

857 

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

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

860 

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

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

863 to override this L{LatLon} class or specify 

864 C{B{LatLon}=None}. 

865 

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

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

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

869 

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

871 

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

873 ''' 

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

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

876 

877 def greatCircle(self, bearing): 

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

879 heading on given compass bearing from this point as its 

880 n-vector. 

881 

882 Direction of vector is such that initial bearing vector 

883 b = c × p. 

884 

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

886 

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

888 

889 @raise Valuerror: Polar coincidence. 

890 

891 @example: 

892 

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

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

895 ''' 

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

897 

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

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

900 

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

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

903 return n.minus(e) 

904 

905 

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

907 

908 

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

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

911 great circle arcs joining consecutive points). 

912 

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

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

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

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

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

918 

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

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

921 

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

923 

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

925 

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

927 and L{ellipsoidalKarney.areaOf}. 

928 

929 @example: 

930 

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

932 >>> areaOf(b) # 8666058750.718977 

933 ''' 

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

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

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

937 n0 = Ps[0].toNvector() 

938 

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

940 p1 = Ps[1] 

941 v1 = p1._N_vector 

942 gc = v2.cross(v1) 

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

944 if w and not Ps.looped: 

945 p2 = _unrollon(p1, p2) 

946 p1 = p2 

947 v2 = p2._N_vector 

948 gc1 = v1.cross(v2) 

949 v1 = v2 

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

951 gc = gc1 

952 

953 if _MODS.booleans.isBoolean(points): 

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

955 else: 

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

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

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

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

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

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

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

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

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

965 

966 

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

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

969 two points or as a point and bearing. 

970 

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

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

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

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

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

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

977 @kwarg radius_exact_height_wrap: Optional keyword arguments, see 

978 method L{LatLon.intersecant2} for further details. 

979 

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

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

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

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{other}} invalid. 

988 

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

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

991 ''' 

992 c = _Nvll.others(center=center) 

993 p = _Nvll.others(point=point) 

994 try: 

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

996 except (TypeError, ValueError) as x: 

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

998 **radius_exact_height_wrap) 

999 

1000 

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

1002 LatLon=LatLon, **LatLon_kwds): 

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

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

1005 

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

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

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

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

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

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

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

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

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

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

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

1017 (L{LatLon}). 

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

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

1020 

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

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

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

1024 

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

1026 

1027 @raise ValueError: Intersection is ambiguous or infinite or 

1028 the lines are parallel, coincident or null. 

1029 

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

1031 

1032 @example: 

1033 

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

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

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

1037 ''' 

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

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

1040 return i.toLatLon(**kwds) 

1041 

1042 

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

1044 LatLon=LatLon, **LatLon_kwds): 

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

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

1047 

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

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

1050 initial bearing at the first start point 

1051 (compass C{degrees360}). 

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

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

1054 initial bearing at the second start point 

1055 (compass C{degrees360}). 

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

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

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

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

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

1061 antipodal points (L{LatLon}). 

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

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

1064 

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

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

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

1068 if available. 

1069 

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

1071 

1072 @raise ValueError: Intersection is ambiguous or infinite or 

1073 the lines are parallel, coincident or null. 

1074 

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

1076 ''' 

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

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

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

1080 

1081 

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

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

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

1085 ''' 

1086 p1 = _Nvll.others(start1=start1) 

1087 p2 = _Nvll.others(start2=start2) 

1088 if wrap: 

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

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

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

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

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

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

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

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

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

1098 

1099 hs = start1.height, start2.height 

1100 # there are two (antipodal) candidate intersection 

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

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

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

1104 

1105 # selection of intersection point depends on how 

1106 # lines are defined (by bearings or endpoints) 

1107 if e1 and e2: # endpoint+endpoint 

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

1109 hs += end1.height, end2.height 

1110 elif e1 and not e2: # endpoint+bearing 

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

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

1113 hs += end1.height, 

1114 elif e2 and not e1: # bearing+endpoint 

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

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

1117 hs += end2.height, 

1118 else: # bearing+bearing 

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

1120 # towards i1, otherwise towards antipodal i2 

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

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

1123 if d1 > 0 and d2 > 0: 

1124 d = 1 # both point to i1 

1125 elif d1 < 0 and d2 < 0: 

1126 d = -1 # both point to i2 

1127 else: # d1, d2 opposite signs 

1128 # intersection is at further-away intersection point, 

1129 # take opposite intersection from mid- point of v1 

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

1131 # get intersection p1 bearing points to, aka being 

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

1133 # function .sphericalTrigonometry._intersect and 

1134 # .ellipsoidalBaseDI._intersect3 

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

1136 

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

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

1139 

1140 

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

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

1143 

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

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

1146 (C{meter}). 

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

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

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

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

1151 

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

1153 

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

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

1156 ''' 

1157 def _N_vs(ps, w): 

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

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

1160 yield p._N_vector 

1161 

1162 try: 

1163 # geographic mean 

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

1165 except (TypeError, ValueError) as x: 

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

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

1168 name=meanOf.__name__)) 

1169 

1170 

1171@deprecated_function 

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

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

1174 

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

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

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

1178 ''' 

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

1180 return r.closest, r.distance 

1181 

1182 

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

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

1185 consecutive points) closest to an other point. 

1186 

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

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

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

1190 

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

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

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

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

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

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

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

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

1199 

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

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

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

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

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

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

1206 C{degrees360}. 

1207 

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

1209 

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

1211 ''' 

1212 _xinstanceof(LatLon, point=point) 

1213 

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

1215 height=height, wrap=wrap) 

1216 

1217 

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

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

1220 great circle arcs joining consecutive points). 

1221 

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

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

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

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

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

1227 

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

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

1230 

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

1232 

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

1234 

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

1236 C{B{points}} a composite. 

1237 

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

1239 and L{sphericalTrigonometry.perimeterOf}. 

1240 ''' 

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

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

1243 p1 = Ps[0] 

1244 v1 = p1._N_vector 

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

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

1247 p2 = _unrollon(p1, p2) 

1248 p1 = p2 

1249 v2 = p2._N_vector 

1250 yield v1.angleTo(v2) 

1251 v1 = v2 

1252 

1253 if _MODS.booleans.isBoolean(points): 

1254 if not closed: 

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

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

1257 else: 

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

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

1260 

1261 

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

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

1264 

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

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

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

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

1269 

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

1271 

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

1273 ''' 

1274 try: 

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

1276 except (TypeError, ValueError) as x: 

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

1278 

1279 

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

1281 height=None, wrap=False, 

1282 LatLon=LatLon, **LatLon_kwds): 

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

1284 from those points. 

1285 

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

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

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

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

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

1291 the mean height (C{meter}). 

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

1293 (C{bool}). 

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

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

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

1297 

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

1299 

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

1301 

1302 @raise Valuerror: Points coincide. 

1303 

1304 @example: 

1305 

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

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

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

1309 ''' 

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

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

1312 height=height, wrap=wrap, 

1313 LatLon=LatLon, **LatLon_kwds) 

1314 

1315 

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

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

1318 LatLon=LatLon, **LatLon_kwds): 

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

1320 

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

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

1323 as B{C{radius}}). 

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

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

1326 as B{C{radius}}). 

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

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

1329 as B{C{radius}}). 

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

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

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

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

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

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

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

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

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

1339 

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

1341 

1342 @raise IntersectionError: No intersection, trilateration failed. 

1343 

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

1345 

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

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

1348 

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

1350 ''' 

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

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

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

1354 radius=radius, height=height, useZ=useZ, 

1355 wrap=wrap, LatLon=LatLon, **LatLon_kwds) 

1356 

1357 

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

1359 areaOf, # functions 

1360 intersecant2, intersection, intersection2, ispolar, 

1361 meanOf, 

1362 nearestOn2, nearestOn3, 

1363 perimeterOf, 

1364 sumOf, 

1365 triangulate, trilaterate) 

1366 

1367# **) MIT License 

1368# 

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

1370# 

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

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

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

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

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

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

1377# 

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

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

1380# 

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

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

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

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

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

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

1387# OTHER DEALINGS IN THE SOFTWARE.