1

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

3

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

5

6 Ellipsoidal classes geodetic L{LatLon}, geocentric (ECEF) L{Cartesian}

7 and C{Nvector} and DEPRECATED L{Ned} and functions L{meanOf}, L{sumOf}

8 and DEPRECATED L{toNed}.

9

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

11 methods by I{(C) Chris Veness 2011-2016} published under the same MIT

12 Licence**, see U{Vector-based geodesy

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

14

15 These classes and functions work with: (a) geodetic lat-/longitude points on

16 the earth's surface and (b) 3-D vectors used as n-vectors representing points

17 on the earth's surface or vectors normal to the plane of a great circle.

18

19 See also I{Kenneth Gade} U{'A Non-singular Horizontal Position Representation'

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

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

22'''

23# make sure int/int division yields float quotient, see .basics

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

25

26 from pygeodesy.basics import issubclassof, map2, _xinstanceof

27 from pygeodesy.datums import _earth_ellipsoid, _ellipsoidal_datum, _WGS84

28# from pygeodesy.dms import toDMS # _MODS

29 from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase, \

30 _nearestOn, LatLonEllipsoidalBase, \

31 _TOL_M, _Wrap

32 from pygeodesy.errors import _IsnotError, _xkwds, _xkwds_pop2

33# from pygeodesy.fmath import fdot # from .nvectorBase

34 from pygeodesy.interns import _Nv00_, _COMMASPACE_, _pole_ # PYCHOK used!

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

36# from pygeodesy.ltp import Ltp # _MODS

37 from pygeodesy.ltpTuples import Aer as _Aer, Ned as _Ned, Ned4Tuple, \

38 sincos2d_, _xnamed

39# from pygeodesy.named import _xnamed # from .ltpTuples

40 from pygeodesy.nvectorBase import fabs, fdot, NorthPole, LatLonNvectorBase, \

41 NvectorBase, sumOf as _sumOf

42 from pygeodesy.props import deprecated_class, deprecated_function, \

43 deprecated_method, Property_RO, property_RO

44 from pygeodesy.streprs import Fmt, fstr, _xzipairs

45 from pygeodesy.units import Bearing, Distance, Height, Scalar

46# from pygeodesy.utily import sincos2d_, _Wrap # from .ltpTuples, .ellipsoidalBase

47

48# from math import fabs # from .nvectorBase

49

50__all__ = _ALL_LAZY.ellipsoidalNvector

51__version__ = '24.05.19'

52

53

54 class Ned(_Ned):

55 '''DEPRECATED on 2024.02.04, use class L{pygeodesy.Ned}.'''

56

57 def __init__(self, north, east, down, **name):

58 deprecated_class(self.__class__)

59 _Ned.__init__(self, north, east, down, **name)

60

61 @deprecated_method # PYCHOK expected

62 def toRepr(self, prec=None, fmt=Fmt.SQUARE, sep=_COMMASPACE_, **unused):

63 '''DEPRECATED, use class L{pygeodesy.Ned}.

64

65 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).

66 @kwarg fmt: Enclosing backets format (C{str}).

67 @kwarg sep: Separator between NEDs (C{str}).

68

69 @return: This Ned as "[L:f, B:degrees360, E:degrees90]" (C{str})

70 with length or slantrange C{L}, bearing or azimuth C{B}

71 and elevation C{E}.

72 '''

73 dms = _MODS.dms

74 t = (fstr(self.slantrange, prec=prec),

75 dms.toDMS(self.azimuth, form=dms.F_D, prec=prec, ddd=0),

76 dms.toDMS(self.elevation, form=dms.F_D, prec=prec, ddd=0))

77 return _xzipairs('LBE', t, sep=sep, fmt=fmt)

78

79

80 class Cartesian(CartesianEllipsoidalBase):

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

82 C{Nvector} and n-vector-based, geodetic L{LatLon}.

83 '''

84 @property_RO

85 def Ecef(self):

86 '''Get the ECEF I{class} (L{EcefVeness}), I{once}.

87 '''

88 return _Ecef()

89

90 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon, datum=None

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

92

93 @kwarg LatLon_and_kwds: Optional L{LatLon}, B{C{datum}} and other

94 keyword arguments. Use C{B{LatLon}=...} to

95 override this L{LatLon} class or specify

96 C{B{LatLon} is None}.

97

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

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

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

101

102 @raise TypeError: Invalid B{C{LatLon_and_kwds}}.

103 '''

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

105 return CartesianEllipsoidalBase.toLatLon(self, **kwds)

106

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

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

109

110 @kwarg Nvector_and_kwds: Optional C{Nvector}, B{C{datum}} and other

111 keyword arguments. Use C{B{Nvector}=...} to

112 override this C{Nvector} class or specify

113 C{B{Nvector} is None}.

114

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

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

117

118 @raise TypeError: Invalid B{C{Nvector_and_kwds}}.

119 '''

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

121 return CartesianEllipsoidalBase.toNvector(self, **kwds)

122

123

124 class LatLon(LatLonNvectorBase, LatLonEllipsoidalBase):

125 '''An n-vector-based, ellipsoidal L{LatLon} point.

126 '''

127 _Nv = None # cached toNvector (C{Nvector})

128

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

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

131 '''

132 if updated:

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

134 LatLonEllipsoidalBase._update(self, updated, *attrs, **setters)

135

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

137# '''Return the (signed) distance from this point to the great

138# circle defined by a start point and an end point or bearing.

139#

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

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

142# initial bearing (compass C{degrees360}) at the

143# start point.

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

145#

146# @return: Distance to great circle, negative if to left or

147# positive if to right of line (C{meter}, same units

148# as B{C{radius}}).

149#

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

151# '''

152# self.others(start=start)

153#

154# if _isDegrees(end): # gc from point and bearing

155# gc = start.greatCircle(end)

156# else: # gc by two points

157# gc = start.toNvector().cross(end.toNvector())

158#

159# # (signed) angle between point and gc normal vector

160# v = self.toNvector()

161# a = gc.angleTo(v, vSign=v.cross(gc))

162# a = _copysign(PI_2, a) - a

163# return a * float(radius)

164

165 def deltaTo(self, other, wrap=False, **Ned_and_kwds):

166 '''Calculate the local delta from this to an other point.

167

168 @note: This is a linear delta, I{unrelated} to a geodesic on the

169 ellipsoid.

170

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

172 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{other}}

173 point (C{bool}).

174 @kwarg Ned_and_kwds: Optional C{B{Ned}=L{Ned} class and B{name}=NN}

175 to return delta and other B{C{Ned}} keyword arguments.

176

177 @return: Delta from this to the other point (B{C{Ned}}).

178

179 @raise TypeError: The B{C{other}} point is not L{LatLon} or

180 B{C{Ned}} is not L{pygeodesy.Ned} nor

181 L{pygeodesy.Ned4Tuple} nor DEPRECATED L{Ned}.

182

183 @raise ValueError: If ellipsoids are incompatible.

184 '''

185 self.ellipsoids(other) # throws TypeError and ValueError

186

187 p = self.others(other)

188 if wrap:

189 p = _Wrap.point(p)

190 # get delta in cartesian frame

191 dc = p.toCartesian().minus(self.toCartesian())

192 # rotate dc to get delta in n-vector reference

193 # frame using the rotation matrix row vectors

194 ned_ = map2(dc.dot, self._rotation3)

195

196 N, kwds = _xkwds_pop2(Ned_and_kwds, Ned=Ned)

197 if issubclassof(N, Ned4Tuple):

198 ned_ += _MODS.ltp.Ltp(self, ecef=self.Ecef(self.datum)),

199 elif not issubclassof(N, _Ned):

200 raise _IsnotError(Fmt.sub_class(_Ned, Ned4Tuple), Ned=N)

201 return N(*ned_, **_xkwds(kwds, name=self.name))

202

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

204# '''Return the destination point after traveling from this

205# point the given distance on the given initial bearing.

206#

207# @arg distance: Distance traveled (C{meter}, same units as

208# given earth B{C{radius}}).

209# @arg bearing: Initial bearing (compass C{degrees360}).

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

211# @kwarg height: Optional height at destination point,

212# overriding default (C{meter}, same units

213# as B{C{radius}}).

214#

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

216# '''

217# r = _m2radians(distance, radius) # angular distance in radians

218# # great circle by starting from this point on given bearing

219# gc = self.greatCircle(bearing)

220#

221# v1 = self.toNvector()

222# x = v1.times(cos(r)) # component of v2 parallel to v1

223# y = gc.cross(v1).times(sin(r)) # component of v2 perpendicular to v1

224#

225# v2 = x.plus(y).unit()

226# return v2.toLatLon(height=self.height if height is C{None} else height)

227

228 def destinationNed(self, delta):

229 '''Calculate the destination point using the supplied NED delta

230 from this point.

231

232 @arg delta: Delta from this to the other point in the local

233 tangent plane (LTP) of this point (L{Ned}).

234

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

236

237 @raise TypeError: If B{C{delta}} is not L{pygeodesy.Ned} or

238 DEPRECATED L{Ned}.

239 '''

240 _xinstanceof(_Ned, delta=delta)

241

242 nv, ev, dv = self._rotation3

243 # convert NED delta to standard coordinate frame of n-vector

244 dn = delta.ned[:3] # XXX Ned4Tuple.to3Tuple

245 # rotate dn to get delta in cartesian (ECEF) coordinate

246 # reference frame using the rotation matrix column vectors

247 dc = Cartesian(fdot(dn, nv.x, ev.x, dv.x),

248 fdot(dn, nv.y, ev.y, dv.y),

249 fdot(dn, nv.z, ev.z, dv.z))

250

251 # apply (cartesian) delta to this Cartesian to obtain destination as cartesian

252 v = self.toCartesian().plus(dc)

253 return v.toLatLon(datum=self.datum, LatLon=self.classof) # Cartesian(v.x, v.y, v.z).toLatLon(...)

254

255 def distanceTo(self, other, radius=None, wrap=False):

256 '''I{Approximate} the distance from this to an other point.

257

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

259 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},

260 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or

261 L{a_f2Tuple}), overriding the mean radius C{R1}

262 of this point's datum..

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

264 B{C{other}} and angular distance (C{bool}).

265

266 @return: Distance (C{meter}, same units as B{C{radius}}).

267

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

269

270 @raise ValueError: Invalid B{C{radius}}.

271 '''

272 p = self.others(other)

273 if wrap:

274 p = _Wrap.point(p)

275 a = self._N_vector.angleTo(p._N_vector, wrap=wrap)

276 E = self.datum.ellipsoid if radius is None else _earth_ellipsoid(radius)

277 return fabs(a) * E.R1 # see .utily.radians2m

278

279 @property_RO

280 def Ecef(self):

281 '''Get the ECEF I{class} (L{EcefVeness}), I{once}.

282 '''

283 return _Ecef()

284

285 @deprecated_method

286 def equals(self, other, eps=None): # PYCHOK no cover

287 '''DEPRECATED, use method L{isequalTo}.

288 '''

289 return self.isequalTo(other, eps=eps)

290

291 def isequalTo(self, other, eps=None):

292 '''Compare this point with an other point.

293

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

295 @kwarg eps: Optional margin (C{float}).

296

297 @return: C{True} if points are identical, including

298 datum, I{ignoring height}, C{False} otherwise.

299

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

301

302 @raise ValueError: Invalid B{C{eps}}.

303

304 @see: Method C{isequalTo3} to include I{height}.

305 '''

306 return self.datum == self.others(other).datum and \

307 _MODS.formy._isequalTo(self, other, eps=eps)

308

309# def greatCircle(self, bearing):

310# '''Return the great circle heading on the given bearing

311# from this point.

312#

313# Direction of vector is such that initial bearing vector

314# b = c × p, where p is representing this point.

315#

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

317#

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

319# '''

320# a, b, _ = self.philamheight

321# t = radians(bearing)

322#

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

324# return self._xnamed(Nvector(sb * ct - sa * cb * st,

325# -cb * ct - sa * sb * st,

326# ca * st)

327

328# def initialBearingTo(self, other, wrap=False):

329# '''Return the initial bearing (forward azimuth) from

330# this to an other point.

331#

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

333# @kwarg wrap: If C{True}, wrap or I{normalize}

334# and unroll the B{C{other}} (C{bool}).

335#

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

337#

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

339# '''

340# p = self.others(other)

341# if wrap:

342# p = _Wrap.point(p)

343# v1 = self.toNvector()

344#

345# gc1 = v1.cross(p.toNvector()) # gc through v1 & v2

346# gc2 = v1.cross(_NP3) # gc through v1 & North pole

347#

348# # bearing is (signed) angle between gc1 & gc2

349# return degrees360(gc1.angleTo(gc2, vSign=v1))

350

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

352 '''Return the point at given fraction between this and

353 an other point.

354

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

356 @arg fraction: Fraction between both points (C{scalar},

357 0.0 at this to 1.0 at the other point.

358 @kwarg height: Optional height, overriding the fractional

359 height (C{meter}).

360 @kwarg wrap: If C{True}, wrap or I{normalize} the

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

362

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

364

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

366 '''

367 p = self.others(other)

368 if wrap:

369 p = _Wrap.point(p)

370 f = Scalar(fraction=fraction)

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

372 i = self.toNvector().intermediateTo(p.toNvector(), f)

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

374

375 @Property_RO

376 def _rotation3(self):

377 '''(INTERNAL) Get the rotation matrix from n-vector coordinate frame axes.

378 '''

379 nv = self.toNvector() # local (n-vector) coordinate frame

380

381 dv = nv.negate() # down, opposite to n-vector

382 ev = NorthPole.cross(nv, raiser=_pole_).unit() # east, pointing perpendicular to the plane

383 nv = ev.cross(dv) # north, by right hand rule

384 return nv, ev, dv # matrix rows

385

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

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

388

389 @kwarg Cartesian_and_kwds: Optional L{Cartesian}, B{C{datum}} and other

390 keyword arguments. Use C{B{Cartesian}=...}

391 to override this L{Cartesian} class or specify

392 C{B{Cartesian}=None}.

393

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

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

396 datum)} with C{C} and C{M} if available.

397

398 @raise TypeError: Invalid B{C{Cartesian}} or other B{C{Cartesian_and_kwds}}.

399 '''

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

401 return LatLonEllipsoidalBase.toCartesian(self, **kwds)

402

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

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

405

406 @kwarg Nvector_and_kwds: Optional C{Nvector}, B{C{datum}} and other

407 keyword arguments. Use C{B{Nvector}=...}

408 to override this C{Nvector} class or specify

409 C{B{Nvector}=None}.

410

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

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

413

414 @raise TypeError: Invalid B{C{Nvector}} or other B{C{Nvector_and_kwds}}.

415 '''

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

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

418

419

420_Nvll = LatLon(0, 0, name=_Nv00_) # reference instance (L{LatLon})

421

422

423 class Nvector(NvectorBase):

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

425 normal to the earth ellipsoid. Unlike lat-/longitude points,

426 n-vectors have no singularities or discontinuities.

427

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

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

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

431

432 Note commonality with L{pygeodesy.sphericalNvector.Nvector}.

433 '''

434 _datum = _WGS84 # default datum (L{Datum})

435

436 def __init__(self, x_xyz, y=None, z=None, h=0, datum=None, ll=None, **name):

437 '''New n-vector normal to the earth's surface.

438

439 @arg x_xyz: X component of vector (C{scalar}) or (3-D) vector

440 (C{Nvector}, L{Vector3d}, L{Vector3Tuple} or

441 L{Vector4Tuple}).

442 @kwarg y: Y component of vector (C{scalar}), ignored if B{C{x_xyz}}

443 is not C{scalar}, otherwise same units as B{C{x_xyz}}.

444 @kwarg z: Z component of vector (C{scalar}), ignored if B{C{x_xyz}}

445 is not C{scalar}, otherwise same units as B{C{x_xyz}}.

446 @kwarg h: Optional height above model surface (C{meter}).

447 @kwarg datum: Optional datum this n-vector is defined in

448 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or

449 L{a_f2Tuple}).

450 @kwarg ll: Optional, original latlon (C{LatLon}).

451 @kwarg name: Optional C{B{name}=NN} (C{str}).

452

453 @raise TypeError: If B{C{datum}} is not a L{Datum}.

454 '''

455 NvectorBase.__init__(self, x_xyz, y=y, z=z, h=h, ll=ll, **name)

456 if datum not in (None, self._datum):

457 self._datum = _ellipsoidal_datum(datum, **name)

458

459 @Property_RO

460 def datum(self):

461 '''Get this n-vector's datum (L{Datum}).

462 '''

463 return self._datum

464

465 @property_RO

466 def ellipsoidalNvector(self):

467 '''Get this C{Nvector}'s ellipsoidal class.

468 '''

469 return type(self)

470

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

472 '''Convert this n-vector to C{Nvector}-based cartesian (ECEF) coordinates.

473

474 @kwarg Cartesian_and_kwds: Optional L{Cartesian}, B{C{h}}, B{C{datum}} and

475 other keyword arguments. Use C{B{Cartesian}=...}

476 to override this L{Cartesian} class or specify

477 C{B{Cartesian} is None}.

478

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

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

481 datum)} with C{C} and C{M} if available.

482

483 @raise TypeError: Invalid B{C{Cartesian_and_kwds}}.

484 '''

485 kwds = _xkwds(Cartesian_and_kwds, h=self.h, Cartesian=Cartesian,

486 datum=self.datum)

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

488

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

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

491

492 @kwarg LatLon_and_kwds: Optional L{LatLon}, B{C{height}}, B{C{datum}}

493 and other keyword arguments. Use C{B{LatLon}=...}

494 to override this L{LatLon} class or specify

495 C{B{LatLon} is None}.

496

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

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

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

500

501 @raise TypeError: Invalid B{C{LatLon_and_kwds}}.

502 '''

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

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

505

506 def unit(self, ll=None):

507 '''Normalize this vector to unit length.

508

509 @kwarg ll: Optional, original latlon (C{LatLon}).

510

511 @return: Normalised vector (C{Nvector}).

512 '''

513 u = NvectorBase.unit(self, ll=ll)

514 if u.datum != self.datum:

515 u._update(False, datum=self.datum)

516 return u

517

518

519 def _Ecef():

520 # return the Ecef class and overwrite property_RO

521 Cartesian.Ecef = LatLon.Ecef = E = _MODS.ecef.EcefVeness

522 return E

523

524

525 def meanOf(points, datum=_WGS84, height=None, wrap=False,

526 **LatLon_and_kwds):

527 '''Compute the geographic mean of several points.

528

529 @arg points: Points to be averaged (L{LatLon}[]).

530 @kwarg datum: Optional datum to use (L{Datum}).

531 @kwarg height: Optional height at mean point, overriding

532 the mean height (C{meter}).

533 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{points}}

534 (C{bool}).

535 @kwarg LatLon_and_kwds: Optional B{C{LatLon}} class to return

536 the mean points and overriding this L{LatLon}

537 (or C{None}) and additional B{C{LatLon}}

538 keyword arguments, ignored if C{B{LatLon}

539 is None}.

540

541 @return: Geographic mean point and mean height (B{C{LatLon}})

542 or if B{C{LatLon}} is C{None}, an L{Ecef9Tuple}C{(x,

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

544 C{M} if available.

545

546 @raise ValueError: Insufficient number of B{C{points}}.

547 '''

548 Ps = _Nvll.PointsIter(points, wrap=wrap)

549 # geographic mean

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

551 kwds = _xkwds(LatLon_and_kwds, height=height, datum=datum,

552 LatLon=LatLon, name__=meanOf)

553 return m.toLatLon(**kwds)

554

555

556 def nearestOn(point, point1, point2, within=True, height=None, wrap=False,

557 equidistant=None, tol=_TOL_M, LatLon=LatLon, **LatLon_kwds):

558 '''I{Iteratively} locate the closest point on the geodesic between

559 two other (ellipsoidal) points.

560

561 @arg point: Reference point (C{LatLon}).

562 @arg point1: Start point of the geodesic (C{LatLon}).

563 @arg point2: End point of the geodesic (C{LatLon}).

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

565 B{C{point1}} and B{C{point2}}, otherwise the

566 closest point elsewhere on the geodesic (C{bool}).

567 @kwarg height: Optional height for the closest point (C{meter},

568 conventionally) or C{None} or C{False} for the

569 interpolated height. If C{False}, the closest

570 takes the heights of the points into account.

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

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

573 @kwarg equidistant: An azimuthal equidistant projection (I{class}

574 or function L{pygeodesy.equidistant}) or C{None}

575 for the preferred C{B{point}.Equidistant}.

576 @kwarg tol: Convergence tolerance (C{meter}).

577 @kwarg LatLon: Optional class to return the closest point

578 (L{LatLon}) or C{None}.

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

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

581

582 @return: Closest point, a B{C{LatLon}} instance or if C{B{LatLon}

583 is None}, a L{LatLon4Tuple}C{(lat, lon, height, datum)}.

584

585 @raise ImportError: Package U{geographiclib

586 <https://PyPI.org/project/geographiclib>}

587 not installed or not found.

588

589 @raise TypeError: Invalid or non-ellipsoidal B{C{point}}, B{C{point1}}

590 or B{C{point2}} or invalid B{C{equidistant}}.

591

592 @raise ValueError: No convergence for the B{C{tol}}.

593

594 @see: U{The B{ellipsoidal} case<https://GIS.StackExchange.com/questions/48937/

595 calculating-intersection-of-two-circles>} and U{Karney's paper

596 <https://ArXiv.org/pdf/1102.1215.pdf>}, pp 20-21, section B{14. MARITIME

597 BOUNDARIES} for more details about the iteration algorithm.

598 '''

599 return _nearestOn(point, point1, point2, within=within, height=height, wrap=wrap,

600 equidistant=equidistant, tol=tol, LatLon=LatLon, **LatLon_kwds)

601

602

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

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

605

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

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

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

609 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword

610 arguments, ignored if C{B{Vector} is None}.

611

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

613

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

615 '''

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

617

618

619@deprecated_function

620 def toNed(distance, bearing, elevation, Ned=Ned, **name):

621 '''DEPRECATED, use L{pygeodesy.Aer}C{(bearing, elevation,

622 distance).xyzLocal.toNed(B{Ned}, name=B{name})} or

623 L{XyzLocal}C{(pygeodesy.Aer(bearing, elevation,

624 distance)).toNed(B{Ned}, name=B{name})}.

625

626 Create an NED vector from distance, bearing and elevation

627 (in local coordinate system).

628

629 @arg distance: NED vector length (C{meter}).

630 @arg bearing: NED vector bearing (compass C{degrees360}).

631 @arg elevation: NED vector elevation from local coordinate

632 frame horizontal (C{degrees}).

633 @kwarg Ned: Optional class to return the NED (C{Ned}) or

634 C{None}.

635 @kwarg name: Optional C{B{name}=NN} (C{str}).

636

637 @return: An NED vector equivalent to this B{C{distance}},

638 B{C{bearing}} and B{C{elevation}} (DEPRECATED L{Ned})

639 or a DEPRECATED L{Ned3Tuple}C{(north, east, down)}

640 if C{B{Ned} is None}.

641

642 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}}

643 or B{C{elevation}}.

644 '''

645 if True: # use new Aer class

646 n, e, d, _ = _Aer(bearing, elevation, distance).xyz4

647 else: # DEPRECATED

648 d = Distance(distance)

649

650 sb, cb, se, ce = sincos2d_(Bearing(bearing),

651 Height(elevation=elevation))

652 n = cb * d * ce

653 e = sb * d * ce

654 d *= se

655

656 r = _MODS.deprecated.classes.Ned3Tuple(n, e, -d) if Ned is None else \

657 Ned(n, e, -d)

658 return _xnamed(r, name)

659

660

661__all__ += _ALL_OTHER(Cartesian, LatLon, Ned, Nvector, # classes

662 meanOf, sumOf, toNed)

663

664# **) MIT License

665#

666# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

667#

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

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

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

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

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

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

674#

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

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

677#

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

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

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

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

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

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

684# OTHER DEALINGS IN THE SOFTWARE.