1

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

3

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

5

6 Ellipsoidal classes geodetic (lat-/longitude) L{LatLon}, geocentric

7(ECEF) L{Cartesian}, DEPRECATED L{Ned} and L{Nvector} and functions

8 L{meanOf}, L{sumOf} 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) geodesic (polar) lat-/longitude

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

17 representing points on the earth's surface or vectors normal to the plane

18 of a great circle.

19

20 See also Kenneth Gade U{'A Non-singular Horizontal Position Representation'

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

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

23'''

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

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

26

27 from pygeodesy.basics import issubclassof, map2, _xinstanceof

28 from pygeodesy.datums import _ellipsoidal_datum, _spherical_datum, _WGS84

29# from pygeodesy.dms import toDMS # _MODS

30 from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase, _TOL_M, \

31 LatLonEllipsoidalBase, _nearestOn

32 from pygeodesy.errors import _IsnotError, _xkwds

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

34 from pygeodesy.interns import NN, _Nv00_, _COMMASPACE_

35 from pygeodesy.interns import _down_, _east_, _north_, _pole_ # PYCHOK used!

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

37# from pygeodesy.ltp import Ltp # _MODS

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

39 sincos2d_, _xnamed

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

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

42 NvectorBase, sumOf as _sumOf

43 from pygeodesy.props import deprecated_class, deprecated_function, \

44 deprecated_method, Property_RO

45 from pygeodesy.streprs import Fmt, fstr, _xzipairs

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

47# from pygeodesy.utily import sincos2d_ # from .ltpTuples

48

49# from math import fabs # from .nvectorBase

50

51__all__ = _ALL_LAZY.ellipsoidalNvector

52__version__ = '23.04.11'

53

54

55 class Ned(_Ned):

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

57

58 def __init__(self, north, east, down, name=NN):

59 deprecated_class(self.__class__)

60 _Ned.__init__(self, north, east, down, name=name)

61

62 @deprecated_method # PYCHOK expected

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

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

65

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

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

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

69

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

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

72 and elevation C{E}.

73 '''

74 dms = _MODS.dms

75 t = (fstr(self.slantrange, prec=3 if prec is None else prec),

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

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

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

79

80

81 class Cartesian(CartesianEllipsoidalBase):

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

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

84 '''

85 @Property_RO

86 def Ecef(self):

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

88 '''

89 return _MODS.ecef.EcefVeness

90

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

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

93

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

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

96 override this L{LatLon} class or specify

97 C{B{LatLon} is None}.

98

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

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

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

102

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

104 '''

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

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

107

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

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

110

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

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

113 override this L{Nvector} class or specify

114 C{B{Nvector} is None}.

115

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

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

118

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

120

121 @example:

122

123 >>> from ellipsoidalNvector import LatLon

124 >>> c = Cartesian(3980581, 97, 4966825)

125 >>> n = c.toNvector() # (0.62282, 0.000002, 0.78237, +0.24)

126 '''

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

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

129

130

131 class LatLon(LatLonNvectorBase, LatLonEllipsoidalBase):

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

133

134 @example:

135

136 >>> from ellipsoidalNvector import LatLon

137 >>> p = LatLon(52.205, 0.119) # height=0, datum=Datums.WGS84

138 '''

139 _Nv = None # cached toNvector (L{Nvector})

140

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

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

143 '''

144 if updated:

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

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

147

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

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

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

151#

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

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

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

155# start point.

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

157#

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

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

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

161#

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

163#

164# @example:

165#

166# >>> p = LatLon(53.2611, -0.7972)

167#

168# >>> s = LatLon(53.3206, -1.7297)

169# >>> b = 96.0

170# >>> d = p.crossTrackDistanceTo(s, b) # -305.7

171#

172# >>> e = LatLon(53.1887, 0.1334)

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

174# '''

175# self.others(start=start)

176#

177# if isscalar(end): # gc from point and bearing

178# gc = start.greatCircle(end)

179# else: # gc by two points

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

181#

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

183# v = self.toNvector()

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

185# a = (-PI_2 - a) if a < 0 else (PI_2 - a)

186# return a * float(radius)

187

188 def deltaTo(self, other, Ned=Ned):

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

190

191 @note: This is a linear delta, I{unrelated} to a geodesic

192 on the ellipsoid.

193

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

195 @kwarg Ned: Class to use (L{pygeodesy.Ned} or

196 L{pygeodesy.Ned4Tuple}), overriding the

197 default DEPRECATED L{Ned}.

198

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

200

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

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

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

204

205 @raise ValueError: If ellipsoids are incompatible.

206

207 @example:

208

209 >>> a = LatLon(49.66618, 3.45063)

210 >>> b = LatLon(48.88667, 2.37472)

211 >>> delta = a.deltaTo(b) # [N:-86126, E:-78900, D:1069]

212 >>> d = delta.length # 116807.681 m

213 >>> b = delta.bearing # 222.493°

214 >>> e = delta.elevation # -0.5245°

215 '''

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

217

218 # get delta in cartesian frame

219 dc = other.toCartesian().minus(self.toCartesian())

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

221 # frame using the rotation matrix row vectors

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

223 if issubclassof(Ned, Ned4Tuple):

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

225 elif not issubclassof(Ned, _Ned):

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

227 return Ned(*ned_, name=self.name)

228

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

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

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

232#

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

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

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

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

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

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

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

240#

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

242#

243# @example:

244#

245# >>> p = LatLon(51.4778, -0.0015)

246# >>> q = p.destination(7794, 300.7)

247# >>> q.toStr() # '51.5135°N, 000.0983°W' ?

248# '''

249# r = _angular(distance, radius) # angular distance in radians

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

251# gc = self.greatCircle(bearing)

252#

253# v1 = self.toNvector()

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

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

256#

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

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

259

260 def destinationNed(self, delta):

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

262 from this point.

263

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

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

266

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

268

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

270 DEPRECATED L{Ned}.

271

272 @example:

273

274 >>> a = LatLon(49.66618, 3.45063)

275 >>> delta = Ned(-86126, -78900, 1069) # from Aer(222.493, -0.5245, 116807.681)

276 >>> b = a.destinationNed(delta) # 48.886669°N, 002.374721°E or 48°53′12.01″N, 002°22′29.0″E +0.20m

277 '''

278 _xinstanceof(_Ned, delta=delta)

279

280 nv, ev, dv = self._rotation3

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

282 dn = delta.ned

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

284 # reference frame using the rotation matrix column vectors

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

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

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

288

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

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

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

292

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

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

295

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

297 @kwarg radius: Mean earth radius, ellipsoid or datum

298 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},

299 L{Datum} or L{a_f2Tuple}), overriding the

300 mean radius C{R1} of this point's datum..

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

302

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

304

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

306

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

308

309 @example:

310

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

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

313 >>> d = p.distanceTo(q) # 404300

314 '''

315 self.others(other)

316

317 a = self._N_vector.angleTo(other._N_vector, wrap=wrap)

318 d = self.datum if radius is None else _spherical_datum(radius)

319 return fabs(a) * d.ellipsoid.R1 # see .utily.radians2m

320

321 @Property_RO

322 def Ecef(self):

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

324 '''

325 return _MODS.ecef.EcefVeness

326

327 @deprecated_method

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

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

330 '''

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

332

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

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

335

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

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

338

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

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

341

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

343

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

345

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

347

348 @example:

349

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

351 >>> q = LatLon(52.205, 0.119)

352 >>> e = p.isequalTo(q) # True

353 '''

354 return LatLonEllipsoidalBase.isequalTo(self, other, eps=eps) \

355 if self.datum == other.datum else False

356

357# def greatCircle(self, bearing):

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

359# from this point.

360#

361# Direction of vector is such that initial bearing vector

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

363#

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

365#

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

367#

368# @example:

369#

370# >>> p = LatLon(53.3206, -1.7297)

371# >>> g = p.greatCircle(96.0)

372# >>> g.toStr() # '(-0.794, 0.129, 0.594)'

373# '''

374# a, b, _ = self.philamheight

375# t = radians(bearing)

376#

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

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

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

380# ca * st)

381

382# def initialBearingTo(self, other):

383# '''Return the initial bearing (forward azimuth) from this

384# to an other point.

385#

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

387#

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

389#

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

391#

392# @example:

393#

394# >>> p1 = LatLon(52.205, 0.119)

395# >>> p2 = LatLon(48.857, 2.351)

396# >>> b = p1.bearingTo(p2) # 156.2

397# '''

398# self.others(other)

399#

400# v1 = self.toNvector()

401# v2 = other.toNvector()

402#

403# gc1 = v1.cross(v2) # gc through v1 & v2

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

405#

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

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

408

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

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

411 an other point.

412

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

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

415 0.0 at this to 1.0 at the other point.

416 @kwarg height: Optional height, overriding the fractional

417 height (C{meter}).

418

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

420

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

422

423 @example:

424

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

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

427 >>> p = p.intermediateTo(q, 0.25) # 51.3721°N, 000.7073°E

428 '''

429 self.others(other)

430

431 f = Scalar(fraction=fraction)

432 i = self.toNvector().intermediateTo(other.toNvector(), f)

433

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

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

436

437 @Property_RO

438 def _rotation3(self):

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

440 '''

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

442

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

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

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

446 return nv, ev, dv # matrix rows

447

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

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

450

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

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

453 to override this L{Cartesian} class or specify

454 C{B{Cartesian} is None}.

455

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

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

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

459

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

461 '''

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

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

464

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

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

467

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

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

470 to override this L{Nvector} class or specify

471 C{B{Nvector} is None}.

472

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

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

475

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

477

478 @example:

479

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

481 >>> n = p.toNvector()

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

483 '''

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

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

486

487

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

489

490

491 class Nvector(NvectorBase):

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

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

494 n-vectors have no singularities or discontinuities.

495

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

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

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

499

500 Note commonality with L{sphericalNvector.Nvector}.

501 '''

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

503

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

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

506

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

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

509 L{Vector4Tuple}).

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

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

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

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

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

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

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

517 L{a_f2Tuple}).

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

519 @kwarg name: Optional name (C{str}).

520

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

522

523 @example:

524

525 >>> from ellipsoidalNvector import Nvector

526 >>> v = Nvector(0.5, 0.5, 0.7071, 1)

527 >>> v.toLatLon() # 45.0°N, 045.0°E, +1.00m

528 '''

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

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

531 self._datum = _ellipsoidal_datum(datum, name=name)

532

533 @Property_RO

534 def datum(self):

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

536 '''

537 return self._datum

538

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

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

541

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

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

544 to override this L{Cartesian} class or specify

545 C{B{Cartesian} is None}.

546

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

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

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

550

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

552

553 @example:

554

555 >>> v = Nvector(0.5, 0.5, 0.7071)

556 >>> c = v.toCartesian() # [3194434, 3194434, 4487327]

557 >>> p = c.toLatLon() # 45.0°N, 45.0°E

558 '''

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

560 datum=self.datum)

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

562

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

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

565

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

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

568 to override this L{LatLon} class or specify

569 C{B{LatLon} is None}.

570

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

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

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

574

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

576

577 @example:

578

579 >>> v = Nvector(0.5, 0.5, 0.7071)

580 >>> p = v.toLatLon() # 45.0°N, 45.0°E

581 '''

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

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

584

585 def unit(self, ll=None):

586 '''Normalize this vector to unit length.

587

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

589

590 @return: Normalised vector (L{Nvector}).

591 '''

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

593 if u.datum != self.datum:

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

595 return u

596

597

598 def meanOf(points, datum=_WGS84, height=None, LatLon=LatLon,

599 **LatLon_kwds):

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

601

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

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

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

605 the mean height (C{meter}).

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

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

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

609 keyword arguments, ignored if

610 C{B{LatLon} is None}.

611

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

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

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

615 C{M} if available.

616

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

618 '''

619 Ps = _Nvll.PointsIter(points)

620 # geographic mean

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

622 kwds = _xkwds(LatLon_kwds, height=height, datum=datum,

623 LatLon=LatLon, name=meanOf.__name__)

624 return m.toLatLon(**kwds)

625

626

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

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

629 '''Iteratively locate the closest point on the geodesic between

630 two other (ellipsoidal) points.

631

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

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

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

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

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

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

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

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

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

641 takes the heights of the points into account.

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

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

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

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

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

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

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

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

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

651

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

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

654

655 @raise ImportError: Package U{geographiclib

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

657 not installed or not found.

658

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

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

661

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

663

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

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

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

667 BOUNDARIES} for more details about the iteration algorithm.

668 '''

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

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

671

672

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

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

675

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

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

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

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

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

681

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

683

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

685 '''

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

687

688

689@deprecated_function

690 def toNed(distance, bearing, elevation, Ned=Ned, name=NN):

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

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

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

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

695

696 Create an NED vector from distance, bearing and elevation

697 (in local coordinate system).

698

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

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

701 @arg elevation: NED vector elevation from local coordinate

702 frame horizontal (C{degrees}).

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

704 C{None}.

705 @kwarg name: Optional name (C{str}).

706

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

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

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

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

711

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

713 or B{C{elevation}}.

714 '''

715 if True: # use new Aer class

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

717 else: # DEPRECATED

718 d = Distance(distance)

719

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

721 Height(elevation=elevation))

722 n = cb * d * ce

723 e = sb * d * ce

724 d *= se

725

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

727 Ned(n, e, -d)

728 return _xnamed(r, name)

729

730

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

732 meanOf, sumOf, toNed)

733

734# **) MIT License

735#

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

737#

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

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

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

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

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

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

744#

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

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

747#

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

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

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

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

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

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

754# OTHER DEALINGS IN THE SOFTWARE.