Coverage for pygeodesy/ellipsoidalNvector.py: 97%
130 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-05-20 11:54 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-05-20 11:54 -0400
2# -*- coding: utf-8 -*-
4u'''Ellipsoidal, C{N-vector}-based geodesy.
6Ellipsoidal classes geodetic (lat-/longitude) L{LatLon}, geocentric
7(ECEF) L{Cartesian}, DEPRECATED L{Ned} and L{Nvector} and functions
8L{meanOf}, L{sumOf} and DEPRECATED L{toNed}.
10Pure Python implementation of n-vector-based geodetic (lat-/longitude)
11methods by I{(C) Chris Veness 2011-2016} published under the same MIT
12Licence**, see U{Vector-based geodesy
13<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>}.
15These classes and functions work with: (a) geodesic (polar) lat-/longitude
16points on the earth's surface and (b) 3-D vectors used as n-vectors
17representing points on the earth's surface or vectors normal to the plane
18of a great circle.
20See also Kenneth Gade U{'A Non-singular Horizontal Position Representation'
21<https://www.NavLab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf>},
22The Journal of Navigation (2010), vol 63, nr 3, pp 395-417.
23'''
24# make sure int/int division yields float quotient, see .basics
25from __future__ import division as _; del _ # PYCHOK semicolon
27from pygeodesy.basics import issubclassof, map2, _xinstanceof
28from pygeodesy.datums import _ellipsoidal_datum, _spherical_datum, _WGS84
29# from pygeodesy.dms import toDMS # _MODS
30from pygeodesy.ellipsoidalBase import CartesianEllipsoidalBase, \
31 _TOL_M, LatLonEllipsoidalBase, \
32 _nearestOn, _Wrap
33from pygeodesy.errors import _IsnotError, _xkwds
34# from pygeodesy.fmath import fdot # from .nvectorBase
35from pygeodesy.interns import NN, _Nv00_, _COMMASPACE_
36from pygeodesy.interns import _down_, _east_, _north_, _pole_ # PYCHOK used!
37from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER
38# from pygeodesy.ltp import Ltp # _MODS
39from pygeodesy.ltpTuples import Aer as _Aer, Ned as _Ned, Ned4Tuple, \
40 sincos2d_, _xnamed
41# from pygeodesy.named import _xnamed # from .ltpTuples
42from pygeodesy.nvectorBase import fabs, fdot, NorthPole, LatLonNvectorBase, \
43 NvectorBase, sumOf as _sumOf
44from pygeodesy.props import deprecated_class, deprecated_function, \
45 deprecated_method, Property_RO
46from pygeodesy.streprs import Fmt, fstr, _xzipairs
47from pygeodesy.units import Bearing, Distance, Height, Scalar
48# from pygeodesy.utily import sincos2d_, _Wrap # from .ltpTuples, .ellipsoidalBase
50# from math import fabs # from .nvectorBase
52__all__ = _ALL_LAZY.ellipsoidalNvector
53__version__ = '23.05.04'
56class Ned(_Ned):
57 '''DEPRECATED, use class L{pygeodesy.Ned}.'''
59 def __init__(self, north, east, down, name=NN):
60 deprecated_class(self.__class__)
61 _Ned.__init__(self, north, east, down, name=name)
63 @deprecated_method # PYCHOK expected
64 def toRepr(self, prec=None, fmt=Fmt.SQUARE, sep=_COMMASPACE_, **unused):
65 '''DEPRECATED, use class L{pygeodesy.Ned}.
67 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
68 @kwarg fmt: Enclosing backets format (C{str}).
69 @kwarg sep: Separator between NEDs (C{str}).
71 @return: This Ned as "[L:f, B:degrees360, E:degrees90]" (C{str})
72 with length or slantrange C{L}, bearing or azimuth C{B}
73 and elevation C{E}.
74 '''
75 dms = _MODS.dms
76 t = (fstr(self.slantrange, prec=3 if prec is None else prec),
77 dms.toDMS(self.azimuth, form=dms.F_D, prec=prec, ddd=0),
78 dms.toDMS(self.elevation, form=dms.F_D, prec=prec, ddd=0))
79 return _xzipairs('LBE', t, sep=sep, fmt=fmt)
82class Cartesian(CartesianEllipsoidalBase):
83 '''Extended to convert geocentric, L{Cartesian} points to
84 L{Nvector} and n-vector-based, geodetic L{LatLon}.
85 '''
86 @Property_RO
87 def Ecef(self):
88 '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
89 '''
90 return _MODS.ecef.EcefVeness
92 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon, datum=None
93 '''Convert this cartesian to an C{Nvector}-based geodetic point.
95 @kwarg LatLon_and_kwds: Optional L{LatLon}, B{C{datum}} and other
96 keyword arguments. Use C{B{LatLon}=...} to
97 override this L{LatLon} class or specify
98 C{B{LatLon} is None}.
100 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set
101 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
102 C, M, datum)} with C{C} and C{M} if available.
104 @raise TypeError: Invalid B{C{LatLon_and_kwds}}.
105 '''
106 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
107 return CartesianEllipsoidalBase.toLatLon(self, **kwds)
109 def toNvector(self, **Nvector_and_kwds): # PYCHOK Datums.WGS84
110 '''Convert this cartesian to L{Nvector} components, I{including height}.
112 @kwarg Nvector_and_kwds: Optional L{Nvector}, B{C{datum}} and other
113 keyword arguments. Use C{B{Nvector}=...} to
114 override this L{Nvector} class or specify
115 C{B{Nvector} is None}.
117 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}}
118 is set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}
120 @raise TypeError: Invalid B{C{Nvector_and_kwds}}.
122 @example:
124 >>> from ellipsoidalNvector import LatLon
125 >>> c = Cartesian(3980581, 97, 4966825)
126 >>> n = c.toNvector() # (0.62282, 0.000002, 0.78237, +0.24)
127 '''
128 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector, datum=self.datum)
129 return CartesianEllipsoidalBase.toNvector(self, **kwds)
132class LatLon(LatLonNvectorBase, LatLonEllipsoidalBase):
133 '''An n-vector-based, ellipsoidal L{LatLon} point.
135 @example:
137 >>> from ellipsoidalNvector import LatLon
138 >>> p = LatLon(52.205, 0.119) # height=0, datum=Datums.WGS84
139 '''
140 _Nv = None # cached toNvector (L{Nvector})
142 def _update(self, updated, *attrs, **setters): # PYCHOK args
143 '''(INTERNAL) Zap cached attributes if updated.
144 '''
145 if updated:
146 LatLonNvectorBase._update(self, updated, _Nv=self._Nv) # special case
147 LatLonEllipsoidalBase._update(self, updated, *attrs, **setters)
149# def crossTrackDistanceTo(self, start, end, radius=R_M):
150# '''Return the (signed) distance from this point to the great
151# circle defined by a start point and an end point or bearing.
152#
153# @arg start: Start point of great circle line (L{LatLon}).
154# @arg end: End point of great circle line (L{LatLon}) or
155# initial bearing (compass C{degrees360}) at the
156# start point.
157# @kwarg radius: Mean earth radius (C{meter}).
158#
159# @return: Distance to great circle, negative if to left or
160# positive if to right of line (C{meter}, same units
161# as B{C{radius}}).
162#
163# @raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.
164#
165# @example:
166#
167# >>> p = LatLon(53.2611, -0.7972)
168#
169# >>> s = LatLon(53.3206, -1.7297)
170# >>> b = 96.0
171# >>> d = p.crossTrackDistanceTo(s, b) # -305.7
172#
173# >>> e = LatLon(53.1887, 0.1334)
174# >>> d = p.crossTrackDistanceTo(s, e) # -307.5
175# '''
176# self.others(start=start)
177#
178# if isscalar(end): # gc from point and bearing
179# gc = start.greatCircle(end)
180# else: # gc by two points
181# gc = start.toNvector().cross(end.toNvector())
182#
183# # (signed) angle between point and gc normal vector
184# v = self.toNvector()
185# a = gc.angleTo(v, vSign=v.cross(gc))
186# a = (-PI_2 - a) if a < 0 else (PI_2 - a)
187# return a * float(radius)
189 def deltaTo(self, other, Ned=Ned, wrap=False):
190 '''Calculate the local delta from this to an other point.
192 @note: This is a linear delta, I{unrelated} to a geodesic
193 on the ellipsoid.
195 @arg other: The other point (L{LatLon}).
196 @kwarg Ned: Class to use (L{pygeodesy.Ned} or
197 L{pygeodesy.Ned4Tuple}), overriding the
198 default DEPRECATED L{Ned}.
199 @kwarg wrap: If C{True}, wrap or I{normalize} the
200 B{C{other}} point (C{bool}).
202 @return: Delta from this to the other point (B{C{Ned}}).
204 @raise TypeError: The B{C{other}} point is not L{LatLon} or
205 B{C{Ned}} is not L{pygeodesy.Ned} nor
206 L{pygeodesy.Ned4Tuple} nor DEPRECATED L{Ned}.
208 @raise ValueError: If ellipsoids are incompatible.
210 @example:
212 >>> a = LatLon(49.66618, 3.45063)
213 >>> b = LatLon(48.88667, 2.37472)
214 >>> delta = a.deltaTo(b) # [N:-86126, E:-78900, D:1069]
215 >>> d = delta.length # 116807.681 m
216 >>> b = delta.bearing # 222.493°
217 >>> e = delta.elevation # -0.5245°
218 '''
219 self.ellipsoids(other) # throws TypeError and ValueError
221 p = self.others(other)
222 if wrap:
223 p = _Wrap.point(p)
224 # get delta in cartesian frame
225 dc = p.toCartesian().minus(self.toCartesian())
226 # rotate dc to get delta in n-vector reference
227 # frame using the rotation matrix row vectors
228 ned_ = map2(dc.dot, self._rotation3)
229 if issubclassof(Ned, Ned4Tuple):
230 ned_ += (_MODS.ltp.Ltp(self, ecef=self.Ecef(self.datum)),)
231 elif not issubclassof(Ned, _Ned):
232 raise _IsnotError(Fmt.sub_class(_Ned, Ned4Tuple), Ned=Ned)
233 return Ned(*ned_, name=self.name)
235# def destination(self, distance, bearing, radius=R_M, height=None):
236# '''Return the destination point after traveling from this
237# point the given distance on the given initial bearing.
238#
239# @arg distance: Distance traveled (C{meter}, same units as
240# given earth B{C{radius}}).
241# @arg bearing: Initial bearing (compass C{degrees360}).
242# @kwarg radius: Mean earth radius (C{meter}).
243# @kwarg height: Optional height at destination point,
244# overriding default (C{meter}, same units
245# as B{C{radius}}).
246#
247# @return: Destination point (L{LatLon}).
248#
249# @example:
250#
251# >>> p = LatLon(51.4778, -0.0015)
252# >>> q = p.destination(7794, 300.7)
253# >>> q.toStr() # '51.5135°N, 000.0983°W' ?
254# '''
255# r = _angular(distance, radius) # angular distance in radians
256# # great circle by starting from this point on given bearing
257# gc = self.greatCircle(bearing)
258#
259# v1 = self.toNvector()
260# x = v1.times(cos(r)) # component of v2 parallel to v1
261# y = gc.cross(v1).times(sin(r)) # component of v2 perpendicular to v1
262#
263# v2 = x.plus(y).unit()
264# return v2.toLatLon(height=self.height if height is C{None} else height)
266 def destinationNed(self, delta):
267 '''Calculate the destination point using the supplied NED delta
268 from this point.
270 @arg delta: Delta from this to the other point in the local
271 tangent plane (LTP) of this point (L{Ned}).
273 @return: Destination point (L{LatLon}).
275 @raise TypeError: If B{C{delta}} is not L{pygeodesy.Ned} or
276 DEPRECATED L{Ned}.
278 @example:
280 >>> a = LatLon(49.66618, 3.45063)
281 >>> delta = Ned(-86126, -78900, 1069) # from Aer(222.493, -0.5245, 116807.681)
282 >>> b = a.destinationNed(delta) # 48.886669°N, 002.374721°E or 48°53′12.01″N, 002°22′29.0″E +0.20m
283 '''
284 _xinstanceof(_Ned, delta=delta)
286 nv, ev, dv = self._rotation3
287 # convert NED delta to standard coordinate frame of n-vector
288 dn = delta.ned
289 # rotate dn to get delta in cartesian (ECEF) coordinate
290 # reference frame using the rotation matrix column vectors
291 dc = Cartesian(fdot(dn, nv.x, ev.x, dv.x),
292 fdot(dn, nv.y, ev.y, dv.y),
293 fdot(dn, nv.z, ev.z, dv.z))
295 # apply (cartesian) delta to this Cartesian to obtain destination as cartesian
296 v = self.toCartesian().plus(dc)
297 return v.toLatLon(datum=self.datum, LatLon=self.classof) # Cartesian(v.x, v.y, v.z).toLatLon(...)
299 def distanceTo(self, other, radius=None, wrap=False):
300 '''I{Approximate} the distance from this to an other point.
302 @arg other: The other point (L{LatLon}).
303 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
304 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or
305 L{a_f2Tuple}), overriding the mean radius C{R1}
306 of this point's datum..
307 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
308 B{C{other}} and angular distance (C{bool}).
310 @return: Distance (C{meter}, same units as B{C{radius}}).
312 @raise TypeError: The B{C{other}} point is not L{LatLon}.
314 @raise ValueError: Invalid B{C{radius}}.
316 @example:
318 >>> p = LatLon(52.205, 0.119)
319 >>> q = LatLon(48.857, 2.351);
320 >>> d = p.distanceTo(q) # 404300
321 '''
322 p = self.others(other)
323 if wrap:
324 p = _Wrap.point(p)
325 a = self._N_vector.angleTo(p._N_vector, wrap=wrap)
326 d = self.datum if radius is None else _spherical_datum(radius)
327 return fabs(a) * d.ellipsoid.R1 # see .utily.radians2m
329 @Property_RO
330 def Ecef(self):
331 '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
332 '''
333 return _MODS.ecef.EcefVeness
335 @deprecated_method
336 def equals(self, other, eps=None): # PYCHOK no cover
337 '''DEPRECATED, use method L{isequalTo}.
338 '''
339 return self.isequalTo(other, eps=eps)
341 def isequalTo(self, other, eps=None):
342 '''Compare this point with an other point.
344 @arg other: The other point (L{LatLon}).
345 @kwarg eps: Optional margin (C{float}).
347 @return: C{True} if points are identical, including
348 datum, I{ignoring height}, C{False} otherwise.
350 @raise TypeError: The B{C{other}} point is not L{LatLon}.
352 @raise ValueError: Invalid B{C{eps}}.
354 @see: Method C{isequalTo3} to include I{height}.
356 @example:
358 >>> p = LatLon(52.205, 0.119)
359 >>> q = LatLon(52.205, 0.119)
360 >>> e = p.isequalTo(q) # True
361 '''
362 return self.datum == self.others(other).datum and \
363 _MODS.formy._isequalTo(self, other, eps=eps)
365# def greatCircle(self, bearing):
366# '''Return the great circle heading on the given bearing
367# from this point.
368#
369# Direction of vector is such that initial bearing vector
370# b = c × p, where p is representing this point.
371#
372# @arg bearing: Bearing from this point (compass C{degrees360}).
373#
374# @return: N-vector representing great circle (L{Nvector}).
375#
376# @example:
377#
378# >>> p = LatLon(53.3206, -1.7297)
379# >>> g = p.greatCircle(96.0)
380# >>> g.toStr() # '(-0.794, 0.129, 0.594)'
381# '''
382# a, b, _ = self.philamheight
383# t = radians(bearing)
384#
385# sa, ca, sb, cb, st, ct = sincos2_(a, b, t)
386# return self._xnamed(Nvector(sb * ct - sa * cb * st,
387# -cb * ct - sa * sb * st,
388# ca * st)
390# def initialBearingTo(self, other, wrap=False):
391# '''Return the initial bearing (forward azimuth) from
392# this to an other point.
393#
394# @arg other: The other point (L{LatLon}).
395# @kwarg wrap: If C{True}, wrap or I{normalize}
396# and unroll the B{C{other}} (C{bool}).
397#
398# @return: Initial bearing (compass C{degrees360}).
399#
400# @raise TypeError: The B{C{other}} point is not L{LatLon}.
401#
402# @example:
403#
404# >>> p1 = LatLon(52.205, 0.119)
405# >>> p2 = LatLon(48.857, 2.351)
406# >>> b = p1.bearingTo(p2) # 156.2
407# '''
408# p = self.others(other)
409# if wrap:
410# p = _Wrap.point(p)
411# v1 = self.toNvector()
412#
413# gc1 = v1.cross(p.toNvector()) # gc through v1 & v2
414# gc2 = v1.cross(_NP3) # gc through v1 & North pole
415#
416# # bearing is (signed) angle between gc1 & gc2
417# return degrees360(gc1.angleTo(gc2, vSign=v1))
419 def intermediateTo(self, other, fraction, height=None, wrap=False):
420 '''Return the point at given fraction between this and
421 an other point.
423 @arg other: The other point (L{LatLon}).
424 @arg fraction: Fraction between both points (C{scalar},
425 0.0 at this to 1.0 at the other point.
426 @kwarg height: Optional height, overriding the fractional
427 height (C{meter}).
428 @kwarg wrap: If C{True}, wrap or I{normalize} the
429 B{C{other}} point (C{bool}).
431 @return: Intermediate point (L{LatLon}).
433 @raise TypeError: The B{C{other}} point is not L{LatLon}.
435 @example:
437 >>> p = LatLon(52.205, 0.119)
438 >>> q = LatLon(48.857, 2.351)
439 >>> p = p.intermediateTo(q, 0.25) # 51.3721°N, 000.7073°E
440 '''
441 p = self.others(other)
442 if wrap:
443 p = _Wrap.point(p)
444 f = Scalar(fraction=fraction)
445 h = self._havg(other, f=f, h=height)
446 i = self.toNvector().intermediateTo(p.toNvector(), f)
447 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)
449 @Property_RO
450 def _rotation3(self):
451 '''(INTERNAL) Get the rotation matrix from n-vector coordinate frame axes.
452 '''
453 nv = self.toNvector() # local (n-vector) coordinate frame
455 dv = nv.negate() # down, opposite to n-vector
456 ev = NorthPole.cross(nv, raiser=_pole_).unit() # east, pointing perpendicular to the plane
457 nv = ev.cross(dv) # north, by right hand rule
458 return nv, ev, dv # matrix rows
460 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian, datum=None
461 '''Convert this point to an C{Nvector}-based geodetic point.
463 @kwarg Cartesian_and_kwds: Optional L{Cartesian}, B{C{datum}} and other
464 keyword arguments. Use C{B{Cartesian}=...}
465 to override this L{Cartesian} class or specify
466 C{B{Cartesian}=None}.
468 @return: The geodetic point (L{Cartesian}) or if B{C{Cartesian}} is set
469 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,
470 datum)} with C{C} and C{M} if available.
472 @raise TypeError: Invalid B{C{Cartesian}} or other B{C{Cartesian_and_kwds}}.
473 '''
474 kwds = _xkwds(Cartesian_and_kwds, Cartesian=Cartesian, datum=self.datum)
475 return LatLonEllipsoidalBase.toCartesian(self, **kwds)
477 def toNvector(self, **Nvector_and_kwds): # PYCHOK signature
478 '''Convert this point to L{Nvector} components, I{including height}.
480 @kwarg Nvector_and_kwds: Optional L{Nvector}, B{C{datum}} and other
481 keyword arguments. Use C{B{Nvector}=...}
482 to override this L{Nvector} class or specify
483 C{B{Nvector}=None}.
485 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}}
486 is set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}.
488 @raise TypeError: Invalid B{C{Nvector}} or other B{C{Nvector_and_kwds}}.
490 @example:
492 >>> p = LatLon(45, 45)
493 >>> n = p.toNvector()
494 >>> n.toStr() # [0.50, 0.50, 0.70710]
495 '''
496 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector, datum=self.datum)
497 return LatLonNvectorBase.toNvector(self, **kwds)
500_Nvll = LatLon(0, 0, name=_Nv00_) # reference instance (L{LatLon})
503class Nvector(NvectorBase):
504 '''An n-vector is a position representation using a (unit) vector
505 normal to the earth ellipsoid. Unlike lat-/longitude points,
506 n-vectors have no singularities or discontinuities.
508 For many applications, n-vectors are more convenient to work
509 with than other position representations like lat-/longitude,
510 earth-centred earth-fixed (ECEF) vectors, UTM coordinates, etc.
512 Note commonality with L{sphericalNvector.Nvector}.
513 '''
514 _datum = _WGS84 # default datum (L{Datum})
516 def __init__(self, x_xyz, y=None, z=None, h=0, datum=None, ll=None, name=NN):
517 '''New n-vector normal to the earth's surface.
519 @arg x_xyz: X component of vector (C{scalar}) or (3-D) vector
520 (C{Nvector}, L{Vector3d}, L{Vector3Tuple} or
521 L{Vector4Tuple}).
522 @kwarg y: Y component of vector (C{scalar}), ignored if B{C{x_xyz}}
523 is not C{scalar}, otherwise same units as B{C{x_xyz}}.
524 @kwarg z: Z component of vector (C{scalar}), ignored if B{C{x_xyz}}
525 is not C{scalar}, otherwise same units as B{C{x_xyz}}.
526 @kwarg h: Optional height above model surface (C{meter}).
527 @kwarg datum: Optional datum this n-vector is defined in
528 (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or
529 L{a_f2Tuple}).
530 @kwarg ll: Optional, original latlon (C{LatLon}).
531 @kwarg name: Optional name (C{str}).
533 @raise TypeError: If B{C{datum}} is not a L{Datum}.
535 @example:
537 >>> from ellipsoidalNvector import Nvector
538 >>> v = Nvector(0.5, 0.5, 0.7071, 1)
539 >>> v.toLatLon() # 45.0°N, 045.0°E, +1.00m
540 '''
541 NvectorBase.__init__(self, x_xyz, y=y, z=z, h=h, ll=ll, name=name)
542 if datum not in (None, self._datum):
543 self._datum = _ellipsoidal_datum(datum, name=name)
545 @Property_RO
546 def datum(self):
547 '''Get this n-vector's datum (L{Datum}).
548 '''
549 return self._datum
551 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian
552 '''Convert this n-vector to C{Nvector}-based cartesian (ECEF) coordinates.
554 @kwarg Cartesian_and_kwds: Optional L{Cartesian}, B{C{h}}, B{C{datum}} and
555 other keyword arguments. Use C{B{Cartesian}=...}
556 to override this L{Cartesian} class or specify
557 C{B{Cartesian} is None}.
559 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} is set
560 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,
561 datum)} with C{C} and C{M} if available.
563 @raise TypeError: Invalid B{C{Cartesian_and_kwds}}.
565 @example:
567 >>> v = Nvector(0.5, 0.5, 0.7071)
568 >>> c = v.toCartesian() # [3194434, 3194434, 4487327]
569 >>> p = c.toLatLon() # 45.0°N, 45.0°E
570 '''
571 kwds = _xkwds(Cartesian_and_kwds, h=self.h, Cartesian=Cartesian,
572 datum=self.datum)
573 return NvectorBase.toCartesian(self, **kwds) # class or .classof
575 def toLatLon(self, **LatLon_and_kwds): # PYCHOK height=None, LatLon=LatLon
576 '''Convert this n-vector to an C{Nvector}-based geodetic point.
578 @kwarg LatLon_and_kwds: Optional L{LatLon}, B{C{height}}, B{C{datum}}
579 and other keyword arguments. Use C{B{LatLon}=...}
580 to override this L{LatLon} class or specify
581 C{B{LatLon} is None}.
583 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set
584 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
585 C, M, datum)} with C{C} and C{M} if available.
587 @raise TypeError: Invalid B{C{LatLon_and_kwds}}.
589 @example:
591 >>> v = Nvector(0.5, 0.5, 0.7071)
592 >>> p = v.toLatLon() # 45.0°N, 45.0°E
593 '''
594 kwds = _xkwds(LatLon_and_kwds, height=self.h, datum=self.datum, LatLon=LatLon)
595 return NvectorBase.toLatLon(self, **kwds) # class or .classof
597 def unit(self, ll=None):
598 '''Normalize this vector to unit length.
600 @kwarg ll: Optional, original latlon (C{LatLon}).
602 @return: Normalised vector (L{Nvector}).
603 '''
604 u = NvectorBase.unit(self, ll=ll)
605 if u.datum != self.datum:
606 u._update(False, datum=self.datum)
607 return u
610def meanOf(points, datum=_WGS84, height=None, wrap=False,
611 **LatLon_and_kwds):
612 '''Compute the geographic mean of several points.
614 @arg points: Points to be averaged (L{LatLon}[]).
615 @kwarg datum: Optional datum to use (L{Datum}).
616 @kwarg height: Optional height at mean point, overriding
617 the mean height (C{meter}).
618 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{points}}
619 (C{bool}).
620 @kwarg LatLon_and_kwds: Optional B{C{LatLon}} class to return
621 the mean points and overriding this L{LatLon}
622 (or C{None}) and additional B{C{LatLon}}
623 keyword arguments, ignored if C{B{LatLon}
624 is None}.
626 @return: Geographic mean point and mean height (B{C{LatLon}})
627 or if B{C{LatLon}} is C{None}, an L{Ecef9Tuple}C{(x,
628 y, z, lat, lon, height, C, M, datum)} with C{C} and
629 C{M} if available.
631 @raise ValueError: Insufficient number of B{C{points}}.
632 '''
633 Ps = _Nvll.PointsIter(points, wrap=wrap)
634 # geographic mean
635 m = sumOf(p._N_vector for p in Ps.iterate(closed=False))
636 kwds = _xkwds(LatLon_and_kwds, height=height, datum=datum,
637 LatLon=LatLon, name=meanOf.__name__)
638 return m.toLatLon(**kwds)
641def nearestOn(point, point1, point2, within=True, height=None, wrap=False,
642 equidistant=None, tol=_TOL_M, LatLon=LatLon, **LatLon_kwds):
643 '''I{Iteratively} locate the closest point on the geodesic between
644 two other (ellipsoidal) points.
646 @arg point: Reference point (C{LatLon}).
647 @arg point1: Start point of the geodesic (C{LatLon}).
648 @arg point2: End point of the geodesic (C{LatLon}).
649 @kwarg within: If C{True} return the closest point I{between}
650 B{C{point1}} and B{C{point2}}, otherwise the
651 closest point elsewhere on the geodesic (C{bool}).
652 @kwarg height: Optional height for the closest point (C{meter},
653 conventionally) or C{None} or C{False} for the
654 interpolated height. If C{False}, the closest
655 takes the heights of the points into account.
656 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll I{only}
657 B{C{point1}} and B{C{point2}} (C{bool}).
658 @kwarg equidistant: An azimuthal equidistant projection (I{class}
659 or function L{pygeodesy.equidistant}) or C{None}
660 for the preferred C{B{point}.Equidistant}.
661 @kwarg tol: Convergence tolerance (C{meter}).
662 @kwarg LatLon: Optional class to return the closest point
663 (L{LatLon}) or C{None}.
664 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
665 arguments, ignored if C{B{LatLon} is None}.
667 @return: Closest point, a B{C{LatLon}} instance or if C{B{LatLon}
668 is None}, a L{LatLon4Tuple}C{(lat, lon, height, datum)}.
670 @raise ImportError: Package U{geographiclib
671 <https://PyPI.org/project/geographiclib>}
672 not installed or not found.
674 @raise TypeError: Invalid or non-ellipsoidal B{C{point}}, B{C{point1}}
675 or B{C{point2}} or invalid B{C{equidistant}}.
677 @raise ValueError: No convergence for the B{C{tol}}.
679 @see: U{The B{ellipsoidal} case<https://GIS.StackExchange.com/questions/48937/
680 calculating-intersection-of-two-circles>} and U{Karney's paper
681 <https://ArXiv.org/pdf/1102.1215.pdf>}, pp 20-21, section B{14. MARITIME
682 BOUNDARIES} for more details about the iteration algorithm.
683 '''
684 return _nearestOn(point, point1, point2, within=within, height=height, wrap=wrap,
685 equidistant=equidistant, tol=tol, LatLon=LatLon, **LatLon_kwds)
688def sumOf(nvectors, Vector=Nvector, h=None, **Vector_kwds):
689 '''Return the vectorial sum of two or more n-vectors.
691 @arg nvectors: Vectors to be added (L{Nvector}[]).
692 @kwarg Vector: Optional class for the vectorial sum (L{Nvector}).
693 @kwarg h: Optional height, overriding the mean height (C{meter}).
694 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword
695 arguments, ignored if C{B{Vector} is None}.
697 @return: Vectorial sum (B{C{Vector}}).
699 @raise VectorError: No B{C{nvectors}}.
700 '''
701 return _sumOf(nvectors, Vector=Vector, h=h, **Vector_kwds)
704@deprecated_function
705def toNed(distance, bearing, elevation, Ned=Ned, name=NN):
706 '''DEPRECATED, use L{pygeodesy.Aer}C{(bearing, elevation,
707 distance).xyzLocal.toNed(B{Ned}, name=B{name})} or
708 L{XyzLocal}C{(pygeodesy.Aer(bearing, elevation,
709 distance)).toNed(B{Ned}, name=B{name})}.
711 Create an NED vector from distance, bearing and elevation
712 (in local coordinate system).
714 @arg distance: NED vector length (C{meter}).
715 @arg bearing: NED vector bearing (compass C{degrees360}).
716 @arg elevation: NED vector elevation from local coordinate
717 frame horizontal (C{degrees}).
718 @kwarg Ned: Optional class to return the NED (C{Ned}) or
719 C{None}.
720 @kwarg name: Optional name (C{str}).
722 @return: An NED vector equivalent to this B{C{distance}},
723 B{C{bearing}} and B{C{elevation}} (DEPRECATED L{Ned})
724 or a DEPRECATED L{Ned3Tuple}C{(north, east, down)}
725 if C{B{Ned} is None}.
727 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}}
728 or B{C{elevation}}.
729 '''
730 if True: # use new Aer class
731 n, e, d, _ = _Aer(bearing, elevation, distance).xyz4
732 else: # DEPRECATED
733 d = Distance(distance)
735 sb, cb, se, ce = sincos2d_(Bearing(bearing),
736 Height(elevation=elevation))
737 n = cb * d * ce
738 e = sb * d * ce
739 d *= se
741 r = _MODS.deprecated.Ned3Tuple(n, e, -d) if Ned is None else \
742 Ned(n, e, -d)
743 return _xnamed(r, name)
746__all__ += _ALL_OTHER(Cartesian, LatLon, Ned, Nvector, # classes
747 meanOf, sumOf, toNed)
749# **) MIT License
750#
751# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
752#
753# Permission is hereby granted, free of charge, to any person obtaining a
754# copy of this software and associated documentation files (the "Software"),
755# to deal in the Software without restriction, including without limitation
756# the rights to use, copy, modify, merge, publish, distribute, sublicense,
757# and/or sell copies of the Software, and to permit persons to whom the
758# Software is furnished to do so, subject to the following conditions:
759#
760# The above copyright notice and this permission notice shall be included
761# in all copies or substantial portions of the Software.
762#
763# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
764# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
765# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
766# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
767# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
768# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
769# OTHER DEALINGS IN THE SOFTWARE.