Coverage for pygeodesy/sphericalNvector.py: 96%

268 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-04-23 16:38 -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 _ValueError, _xError, _xkwds, _xkwds_pop 

41from pygeodesy.fmath import fmean, fsum 

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

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

44 _points_, _pole_ 

45from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER 

46# from pygeodesy.named import notImplemented # from .nvectorBase 

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

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

49 sumOf as _sumOf, _triangulate, _trilaterate, \ 

50 notImplemented 

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

52from pygeodesy.props import deprecated_function, deprecated_method 

53from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, \ 

54 Datums, _intersecant2, LatLonSphericalBase 

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

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

57 

58# from math import atan2, fabs # from utily 

59 

60__all__ = _ALL_LAZY.sphericalNvector 

61__version__ = '23.04.23' 

62 

63_lines_ = 'lines' 

64 

65 

66class Cartesian(CartesianSphericalBase): 

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

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

69 ''' 

70 

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

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

73 

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

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

76 to override this L{LatLon} class or specify 

77 C{B{LatLon}=None}. 

78 

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

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

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

82 

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

84 ''' 

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

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

87 

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

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

90 

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

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

93 to override this L{Nvector} class or specify 

94 C{B{Nvector}=None}. 

95 

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

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

98 

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

100 ''' 

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

102 # datum=datum or self.datum) 

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

104 # return ll.toNvector(**kwds) 

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

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

107 

108 

109class LatLon(LatLonNvectorBase, LatLonSphericalBase): 

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

111 

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

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

114 

115 @example: 

116 

117 >>> from sphericalNvector import LatLon 

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

119 ''' 

120 _Nv = None # cached_toNvector L{Nvector}) 

121 

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

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

124 ''' 

125 if updated: # reset caches 

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

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

128 

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

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

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

132 end point. 

133 

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

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

136 from the start point to the point where the perpendicular 

137 crosses the line. 

138 

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

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

141 initial bearing from start point (compass 

142 C{degrees360}). 

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

144 

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

146 after the start toward the end point of the line 

147 or negative if before the start point). 

148 

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

150 

151 @raise Valuerror: Some points coincide. 

152 

153 @example: 

154 

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

156 

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

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

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

160 ''' 

161 self.others(start=start) 

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

163 

164 p = self.toNvector() 

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

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

167 

168 @deprecated_method 

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

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

171 ''' 

172 return self.initialBearingTo(other) 

173 

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

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

176 defined by a start and end point. 

177 

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

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

180 initial bearing from start point (compass 

181 C{degrees360}). 

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

183 

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

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

186 

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

188 

189 @raise Valuerror: Some points coincide. 

190 

191 @example: 

192 

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

194 

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

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

197 

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

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

200 ''' 

201 self.others(start=start) 

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

203 

204 p = self.toNvector() 

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

206 

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

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

209 travelled the given distance on the given bearing. 

210 

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

212 as B{C{radius}}). 

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

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

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

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

217 

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

219 

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

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

222 

223 @example: 

224 

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

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

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

228 ''' 

229 a = _angular(distance, radius) 

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

231 

232 p = self.toNvector() 

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

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

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

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

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

238 

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

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

241 

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

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

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

245 

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

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

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

249 

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

251 

252 @example: 

253 

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

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

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

257 ''' 

258 self.others(other) 

259 

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

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

262 

263# @Property_RO 

264# def Ecef(self): 

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

266# ''' 

267# return _ALL_MODS.ecef.EcefKarney 

268 

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

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

271 ''' 

272 s = start.toNvector() 

273 if isscalar(end): # bearing 

274 gc = s.greatCircle(end) 

275 e = None 

276 else: 

277 e = self.others(end, name=namend).toNvector() 

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

279 return gc, s, e 

280 

281 def greatCircle(self, bearing): 

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

283 heading on the given bearing from this point. 

284 

285 Direction of vector is such that initial bearing vector 

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

287 

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

289 

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

291 ''' 

292 a, b = self.philam 

293 t = Bearing_(bearing) 

294 

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

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

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

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

299 

300 def greatCircleTo(self, other): 

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

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

303 

304 Direction of vector is such that initial bearing vector 

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

306 

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

308 this point (compass C{degrees360}). 

309 

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

311 

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

313 

314 @raise Valuerror: Points coincide. 

315 

316 @example: 

317 

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

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

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

321 

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

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

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

325 ''' 

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

327 return gc.unit() 

328 

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

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

331 to an other point. 

332 

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

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

335 

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

337 

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

339 point or the C{NorthPole}, provided 

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

341 

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

343 

344 @example: 

345 

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

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

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

349 ''' 

350 self.others(other) 

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

352 n = self.toNvector() 

353# gc1 = self.greatCircleTo(other) 

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

355# gc2 = self.greatCircleTo(NorthPole) 

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

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

358 

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

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

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

362 

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

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

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

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

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

368 

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

370 

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

372 

373 @example: 

374 

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

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

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

378 ''' 

379 self.others(other) 

380 

381 f = Scalar(fraction=fraction) 

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

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

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

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

386 

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

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

389 

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

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

392 other point. 

393 

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

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

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

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

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

399 

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

401 

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

403 

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

405 

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

407 

408 @example: 

409 

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

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

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

413 ''' 

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

415 p = self.toNvector() 

416 f = Scalar(fraction=fraction) 

417 

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

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

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

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

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

423 

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

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

426 

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

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

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

430 

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

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

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

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

435 initial bearing at the second point (compass 

436 C{degrees}). 

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

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

439 

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

441 if no unique intersection exists. 

442 

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

444 is not L{LatLon}. 

445 

446 @raise ValueError: Intersection is ambiguous or infinite or 

447 the lines are parallel, coincident or null. 

448 

449 @example: 

450 

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

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

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

454 ''' 

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

456 height=height, LatLon=self.classof) 

457 

458 def isenclosedBy(self, points): 

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

460 

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

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

463 

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

465 C{False} otherwise. 

466 

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

468 

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

470 

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

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

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

474 ''' 

475 if _MODS.booleans.isBoolean(points): 

476 return points._encloses(self.lat, self.lon) 

477 

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

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

480 def _subtangles(Ps, n0): 

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

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

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

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

485 vs1 = vs2 

486 

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

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

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

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

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

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

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

494 # XXX are winding number optimisations equally applicable to 

495 # spherical surface? 

496 return fabs(s) > PI 

497 

498 @deprecated_method 

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

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

501 return self.isenclosedBy(points) 

502 

503 def iswithin(self, point1, point2): 

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

505 

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

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

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

509 the same hemispere). 

510 

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

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

513 

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

515 C{False} otherwise. 

516 

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

518 ''' 

519 n0 = self.toNvector() 

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

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

522 

523 # corner case, null arc 

524 if n1.isequalTo(n2): 

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

526 

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

528 return False # PYCHOK returns 

529 

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

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

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

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

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

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

536 

537 @deprecated_method 

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

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

540 return self.iswithin(point1, point2) 

541 

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

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

544 

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

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

547 the mean height (C{meter}). 

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

549 may be negative or greater than 1.0. 

550 

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

552 

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

554 

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

556 

557 @example: 

558 

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

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

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

562 ''' 

563 if fraction is _0_5: 

564 self.others(other) 

565 

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

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

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

569 else: 

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

571 return r 

572 

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

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

575 points closest to this point. 

576 

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

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

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

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

581 C{None} to interpolate the height. 

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

583 both given points, otherwise the closest 

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

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

586 

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

588 

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

590 not supported. 

591 

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

593 

594 @example: 

595 

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

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

598 

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

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

601 

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

603 

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

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

606 ''' 

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

608 notImplemented(self, wrap=wrap) 

609 

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

611 # closer to arc than to its endpoints, 

612 # find the closest point on the arc 

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

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

615 n = gc1.cross(gc2) 

616 

617 elif within: # for backward compatibility 

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

619 

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

621 n1 = point1.toNvector() 

622 n2 = point2.toNvector() 

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

624 if n is n1: 

625 return point1 

626 elif n is n2: 

627 return point2 

628 

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

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

631 d = point1.distanceTo(point2) 

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

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

634 return p 

635 

636 # @deprecated_method 

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

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

639 

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

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

642 to that point from this point ... 

643 ''' 

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

645 return r.closest, r.distance 

646 

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

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

649 arcs joining consecutive points) closest to this point. 

650 

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

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

653 

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

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

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

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

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

659 

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

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

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

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

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

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

666 

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

668 

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

670 ''' 

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

672 

673 R = self.distanceTo 

674 N = self.nearestOn 

675 

676 c = p1 = Ps[0] 

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

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

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

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

681 if d < r: 

682 c, r = p, d 

683 p1 = p2 

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

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

686 

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

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

689 

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

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

692 to override this L{Cartesian} class or specify 

693 C{B{Cartesian}=None}. 

694 

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

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

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

698 

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

700 ''' 

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

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

703 

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

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

706 

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

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

709 this L{Nvector} class or specify 

710 C{B{Nvector}=None}. 

711 

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

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

714 

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

716 

717 @example: 

718 

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

720 >>> n = p.toNvector() 

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

722 ''' 

723 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector) 

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

725 

726 

727class Nvector(NvectorBase): 

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

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

730 n-vectors have no singularities or discontinuities. 

731 

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

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

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

735 

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

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

738 

739 Note commonality with L{ellipsoidalNvector.Nvector}. 

740 ''' 

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

742 

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

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

745 (ECEF) coordinates. 

746 

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

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

749 to override this L{Cartesian} class or specify 

750 C{B{Cartesian}=None}. 

751 

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

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

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

755 

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

757 ''' 

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

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

760 

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

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

763 

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

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

766 to override this L{LatLon} class or specify 

767 C{B{LatLon}=None}. 

768 

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

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

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

772 

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

774 

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

776 ''' 

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

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

779 

780 def greatCircle(self, bearing): 

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

782 heading on given compass bearing from this point as its 

783 n-vector. 

784 

785 Direction of vector is such that initial bearing vector 

786 b = c × p. 

787 

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

789 

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

791 

792 @raise Valuerror: Polar coincidence. 

793 

794 @example: 

795 

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

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

798 ''' 

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

800 

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

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

803 

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

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

806 return n.minus(e) 

807 

808 

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

810 

811 

812def areaOf(points, radius=R_M): 

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

814 (with great circle arcs joining consecutive points). 

815 

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

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

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

819 

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

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

822 

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

824 

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

826 

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

828 and L{ellipsoidalKarney.areaOf}. 

829 

830 @example: 

831 

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

833 >>> areaOf(b) # 8666058750.718977 

834 ''' 

835 def _interangles(Ps): 

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

837 n0 = Ps[0].toNvector() 

838 

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

840 v1 = Ps[1]._N_vector 

841 gc = v2.cross(v1) 

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

843 v2 = p._N_vector 

844 gc1 = v1.cross(v2) 

845 v1 = v2 

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

847 gc = gc1 

848 

849 if _MODS.booleans.isBoolean(points): 

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

851 else: 

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

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

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

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

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

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

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

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

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

861 

862 

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

864 height=None, wrap=True): 

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

866 

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

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

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

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

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

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

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

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

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

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

877 circle} methods. 

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

879 conventionally) or C{None}. 

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

881 

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

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

884 point C{is} this very instance. 

885 

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

887 

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

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

890 

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

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

893 ''' 

894 c = _Nvll.others(center=center) 

895 p = _Nvll.others(point=point) 

896 try: 

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

898 height=height, wrap=wrap) 

899 except (TypeError, ValueError) as x: 

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

901 

902 

903def intersection(start1, end1, start2, end2, height=None, 

904 LatLon=LatLon, **LatLon_kwds): 

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

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

907 

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

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

910 initial bearing at the first start point 

911 (compass C{degrees360}). 

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

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

914 initial bearing at the second start point 

915 (compass C{degrees360}). 

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

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

918 @kwarg LatLon: Optional class to return the intersection 

919 point (L{LatLon}). 

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

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

922 

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

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

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

926 

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

928 

929 @raise ValueError: Intersection is ambiguous or infinite or 

930 the lines are parallel, coincident or null. 

931 

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

933 

934 @example: 

935 

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

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

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

939 ''' 

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

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

942 _ = _xkwds_pop(kwds, wrap=None) # from .formy.intersection2 

943 return i.toLatLon(**kwds) 

944 

945 

946def intersection2(start1, end1, start2, end2, height=None, 

947 LatLon=LatLon, **LatLon_kwds): 

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

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

950 

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

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

953 initial bearing at the first start point 

954 (compass C{degrees360}). 

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

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

957 initial bearing at the second start point 

958 (compass C{degrees360}). 

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

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

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

962 antipodal point (L{LatLon}). 

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

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

965 

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

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

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

969 if available. 

970 

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

972 

973 @raise ValueError: Intersection is ambiguous or infinite or 

974 the lines are parallel, coincident or null. 

975 

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

977 ''' 

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

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

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

981 

982 

983def _intersect3(start1, end1, start2, end2, height): 

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

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

986 ''' 

987 _Nvll.others(start1=start1) 

988 _Nvll.others(start2=start2) 

989 

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

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

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

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

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

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

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

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

998 

999 hs = start1.height, start2.height 

1000 # there are two (antipodal) candidate intersection 

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

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

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

1004 

1005 # selection of intersection point depends on how 

1006 # lines are defined (by bearings or endpoints) 

1007 if e1 and e2: # endpoint+endpoint 

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

1009 hs += end1.height, end2.height 

1010 elif e1 and not e2: # endpoint+bearing 

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

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

1013 hs += end1.height, 

1014 elif e2 and not e1: # bearing+endpoint 

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

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

1017 hs += end2.height, 

1018 else: # bearing+bearing 

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

1020 # towards i1, otherwise towards antipodal i2 

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

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

1023 if d1 > 0 and d2 > 0: 

1024 d = 1 # both point to i1 

1025 elif d1 < 0 and d2 < 0: 

1026 d = -1 # both point to i2 

1027 else: # d1, d2 opposite signs 

1028 # intersection is at further-away intersection point, 

1029 # take opposite intersection from mid- point of v1 

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

1031 # get intersection p1 bearing points to, aka being 

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

1033 # function .sphericalTrigonometry._intersect and 

1034 # .ellipsoidalBaseDI._intersect3 

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

1036 

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

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

1039 

1040 

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

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

1043 

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

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

1046 (C{meter}). 

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

1048 (L{LatLon}). 

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

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

1051 

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

1053 

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

1055 

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

1057 ''' 

1058 Ps = _Nvll.PointsIter(points) 

1059 # geographic mean 

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

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

1062 name=meanOf.__name__) 

1063 return m.toLatLon(**kwds) 

1064 

1065 

1066@deprecated_function 

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

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

1069 

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

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

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

1073 ''' 

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

1075 return r.closest, r.distance 

1076 

1077 

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

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

1080 joining consecutive points) closest to an other point. 

1081 

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

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

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

1085 

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

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

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

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

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

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

1092 

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

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

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

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

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

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

1099 C{degrees360}. 

1100 

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

1102 

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

1104 ''' 

1105 _xinstanceof(LatLon, point=point) 

1106 

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

1108 

1109 

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

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

1112 (with great circle arcs joining consecutive points). 

1113 

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

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

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

1117 

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

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

1120 

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

1122 

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

1124 

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

1126 C{B{points}} a composite. 

1127 

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

1129 and L{ellipsoidalKarney.perimeterOf}. 

1130 ''' 

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

1132 v1 = Ps[0]._N_vector 

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

1134 v2 = p._N_vector 

1135 yield v1.angleTo(v2) 

1136 v1 = v2 

1137 

1138 if _MODS.booleans.isBoolean(points): 

1139 if not closed: 

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

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

1142 else: 

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

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

1145 

1146 

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

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

1149 

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

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

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

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

1154 

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

1156 

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

1158 ''' 

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

1160 

1161 

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

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

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

1165 from those points. 

1166 

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

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

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

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

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

1172 the mean height (C{meter}). 

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

1174 (L{LatLon}). 

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

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

1177 

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

1179 

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

1181 

1182 @raise Valuerror: Points coincide. 

1183 

1184 @example: 

1185 

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

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

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

1189 ''' 

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

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

1192 height=height, LatLon=LatLon, **LatLon_kwds) 

1193 

1194 

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

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

1197 LatLon=LatLon, **LatLon_kwds): 

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

1199 

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

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

1202 as B{C{radius}}). 

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

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

1205 as B{C{radius}}). 

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

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

1208 as B{C{radius}}). 

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

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

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

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

1213 @kwarg LatLon: Optional class to return the trilaterated 

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

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

1216 

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

1218 

1219 @raise IntersectionError: No intersection, trilateration failed. 

1220 

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

1222 

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

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

1225 

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

1227 ''' 

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

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

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

1231 radius=radius, height=height, useZ=useZ, 

1232 LatLon=LatLon, **LatLon_kwds) 

1233 

1234 

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

1236 areaOf, # functions 

1237 intersecant2, intersection, intersection2, ispolar, 

1238 meanOf, 

1239 nearestOn2, nearestOn3, 

1240 perimeterOf, 

1241 sumOf, 

1242 triangulate, trilaterate) 

1243 

1244# **) MIT License 

1245# 

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

1247# 

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

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

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

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

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

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

1254# 

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

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

1257# 

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

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

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

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

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

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

1264# OTHER DEALINGS IN THE SOFTWARE.