Coverage for pygeodesy/nvectorBase.py: 97%
233 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-21 13:14 -0400
2# -*- coding: utf-8 -*-
4u'''(INTERNAL) Private elliposiodal and spherical C{Nvector} base classes
5L{LatLonNvectorBase} and L{NvectorBase} and function L{sumOf}.
7Pure Python implementation of C{n-vector}-based geodesy tools for ellipsoidal
8earth models, transcoded from JavaScript originals by I{(C) Chris Veness 2005-2016}
9and published under the same MIT Licence**, see U{Vector-based geodesy
10<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>}.
11'''
13from pygeodesy.basics import len2, map1
14from pygeodesy.constants import EPS, EPS1, EPS_2, R_M, _2_0, _N_2_0
15# from pygeodesy.datums import _spherical_datum # from .formy
16from pygeodesy.errors import IntersectionError, _ValueError, VectorError, \
17 _xkwds, _xkwds_pop
18from pygeodesy.fmath import fdot, fidw, hypot_ # PYCHOK fdot shared
19from pygeodesy.fsums import fsum, fsum_
20from pygeodesy.formy import n_xyz2latlon, n_xyz2philam, _spherical_datum
21from pygeodesy.interns import MISSING, NN, _1_, _2_, _3_, _bearing_, \
22 _coincident_, _COMMASPACE_, _distance_, _h_, \
23 _intersection_, _no_, _NorthPole_, _points_, \
24 _pole_, _SPACE_, _SouthPole_, _UNDER
25from pygeodesy.latlonBase import LatLonBase, _ALL_DOCS, _MODS
26# from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS # from .latlonBase
27from pygeodesy.named import notImplemented, _xother3
28from pygeodesy.namedTuples import Trilaterate5Tuple, Vector3Tuple, \
29 Vector4Tuple
30from pygeodesy.props import deprecated_method, property_doc_, \
31 Property_RO, _update_all
32from pygeodesy.streprs import Fmt, hstr, unstr, _xattrs
33from pygeodesy.units import Bearing, Height, Radius_, Scalar
34# from pygeodesy.utily import sincos2d # from vector3d
35from pygeodesy.vector3d import Vector3d, sumOf as _sumOf, sincos2d, _xyzhdn3
37from math import fabs, sqrt # atan2, cos, sin
39__all__ = (_NorthPole_, _SouthPole_) # constants
40__version__ = '23.03.19'
43class NvectorBase(Vector3d): # XXX kept private
44 '''Base class for ellipsoidal and spherical C{Nvector}s.
45 '''
46 _datum = None # L{Datum}, overriden
47 _h = Height(h=0) # height (C{meter})
48 _H = NN # height prefix (C{str}), '↑' in JS version
50 def __init__(self, x_xyz, y=None, z=None, h=0, ll=None, datum=None, name=NN):
51 '''New n-vector normal to the earth's surface.
53 @arg x_xyz: X component of vector (C{scalar}) or (3-D) vector
54 (C{Nvector}, L{Vector3d}, L{Vector3Tuple} or
55 L{Vector4Tuple}).
56 @kwarg y: Y component of vector (C{scalar}), ignored if B{C{x_xyz}}
57 is not C{scalar}, otherwise same units as B{C{x_xyz}}.
58 @kwarg z: Z component of vector (C{scalar}), ignored if B{C{x_xyz}}
59 is not C{scalar}, otherwise same units as B{C{x_xyz}}.
60 @kwarg h: Optional height above surface (C{meter}).
61 @kwarg ll: Optional, original latlon (C{LatLon}).
62 @kwarg datum: Optional, I{pass-thru} datum (L{Datum}).
63 @kwarg name: Optional name (C{str}).
65 @raise TypeError: Non-scalar B{C{x}}, B{C{y}} or B{C{z}}
66 coordinate or B{C{x}} not an C{Nvector},
67 L{Vector3Tuple} or L{Vector4Tuple} or
68 invalid B{C{datum}}.
70 @example:
72 >>> from pygeodesy.sphericalNvector import Nvector
73 >>> v = Nvector(0.5, 0.5, 0.7071, 1)
74 >>> v.toLatLon() # 45.0°N, 045.0°E, +1.00m
75 '''
76 h, d, n = _xyzhdn3(x_xyz, h, datum, ll)
77 Vector3d.__init__(self, x_xyz, y=y, z=z, ll=ll, name=name or n)
78 if h:
79 self.h = h
80 if d is not None:
81 self._datum = _spherical_datum(d, name=self.name) # pass-thru
83 @Property_RO
84 def datum(self):
85 '''Get the I{pass-thru} datum (C{Datum}) or C{None}.
86 '''
87 return self._datum
89 @Property_RO
90 def Ecef(self):
91 '''Get the ECEF I{class} (L{EcefKarney}), I{lazily}.
92 '''
93 return _MODS.ecef.EcefKarney # default
95 @property_doc_(''' the height above surface (C{meter}).''')
96 def h(self):
97 '''Get the height above surface (C{meter}).
98 '''
99 return self._h
101 @h.setter # PYCHOK setter!
102 def h(self, h):
103 '''Set the height above surface (C{meter}).
105 @raise TypeError: If B{C{h}} invalid.
107 @raise VectorError: If B{C{h}} invalid.
108 '''
109 h = Height(h=h, Error=VectorError)
110 if self._h != h:
111 _update_all(self)
112 self._h = h
114 @property_doc_(''' the height prefix (C{str}).''')
115 def H(self):
116 '''Get the height prefix (C{str}).
117 '''
118 return self._H
120 @H.setter # PYCHOK setter!
121 def H(self, H):
122 '''Set the height prefix (C{str}).
123 '''
124 self._H = str(H) if H else NN
126 def hStr(self, prec=-2, m=NN):
127 '''Return a string for the height B{C{h}}.
129 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
130 @kwarg m: Optional unit of the height (C{str}).
132 @see: Function L{pygeodesy.hstr}.
133 '''
134 return NN(self.H, hstr(self.h, prec=prec, m=m))
136 @Property_RO
137 def isEllipsoidal(self):
138 '''Check whether this n-vector is ellipsoidal (C{bool} or C{None} if unknown).
139 '''
140 return self.datum.isEllipsoidal if self.datum else None
142 @Property_RO
143 def isSpherical(self):
144 '''Check whether this n-vector is spherical (C{bool} or C{None} if unknown).
145 '''
146 return self.datum.isSpherical if self.datum else None
148 @Property_RO
149 def lam(self):
150 '''Get the (geodetic) longitude in C{radians} (C{float}).
151 '''
152 return self.philam.lam
154 @Property_RO
155 def lat(self):
156 '''Get the (geodetic) latitude in C{degrees} (C{float}).
157 '''
158 return self.latlon.lat
160 @Property_RO
161 def latlon(self):
162 '''Get the (geodetic) lat-, longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).
163 '''
164 return n_xyz2latlon(self.x, self.y, self.z, name=self.name)
166 @Property_RO
167 def latlonheight(self):
168 '''Get the (geodetic) lat-, longitude in C{degrees} and height (L{LatLon3Tuple}C{(lat, lon, height)}).
169 '''
170 return self.latlon.to3Tuple(self.h)
172 @Property_RO
173 def latlonheightdatum(self):
174 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}).
175 '''
176 return self.latlonheight.to4Tuple(self.datum)
178 @Property_RO
179 def lon(self):
180 '''Get the (geodetic) longitude in C{degrees} (C{float}).
181 '''
182 return self.latlon.lon
184 @Property_RO
185 def phi(self):
186 '''Get the (geodetic) latitude in C{radians} (C{float}).
187 '''
188 return self.philam.phi
190 @Property_RO
191 def philam(self):
192 '''Get the (geodetic) lat-, longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).
193 '''
194 return n_xyz2philam(self.x, self.y, self.z, name=self.name)
196 @Property_RO
197 def philamheight(self):
198 '''Get the (geodetic) lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
199 '''
200 return self.philam.to3Tuple(self.h)
202 @Property_RO
203 def philamheightdatum(self):
204 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLam4Tuple}C{(phi, lam, height, datum)}).
205 '''
206 return self.philamheight.to4Tuple(self.datum)
208 @deprecated_method
209 def to2ab(self): # PYCHOK no cover
210 '''DEPRECATED, use property L{philam}.
212 @return: A L{PhiLam2Tuple}C{(phi, lam)}.
213 '''
214 return self.philam
216 @deprecated_method
217 def to3abh(self, height=None): # PYCHOK no cover
218 '''DEPRECATED, use property L{philamheight} or C{philam.to3Tuple(B{height})}.
220 @kwarg height: Optional height, overriding this
221 n-vector's height (C{meter}).
223 @return: A L{PhiLam3Tuple}C{(phi, lam, height)}.
225 @raise ValueError: Invalid B{C{height}}.
226 '''
227 return self.philamheight if height in (None, self.h) else \
228 self.philam.to3Tuple(height)
230 def toCartesian(self, h=None, Cartesian=None, datum=None, **Cartesian_kwds):
231 '''Convert this n-vector to C{Nvector}-based cartesian (ECEF) coordinates.
233 @kwarg h: Optional height, overriding this n-vector's height (C{meter}).
234 @kwarg Cartesian: Optional class to return the (ECEF) coordinates
235 (C{Cartesian}).
236 @kwarg datum: Optional datum (C{Datum}), overriding this datum.
237 @kwarg Cartesian_kwds: Optional, additional B{C{Cartesian}} keyword
238 arguments, ignored if C{B{Cartesian} is None}.
240 @return: The cartesian (ECEF) coordinates (B{C{Cartesian}}) or
241 if C{B{Cartesian} is None}, an L{Ecef9Tuple}C{(x, y, z,
242 lat, lon, height, C, M, datum)} with C{C} and C{M} if
243 available.
245 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}}
246 argument.
248 @raise ValueError: Invalid B{C{h}}.
250 @example:
252 >>> v = Nvector(0.5, 0.5, 0.7071)
253 >>> c = v.toCartesian() # [3194434, 3194434, 4487327]
254 >>> p = c.toLatLon() # 45.0°N, 45.0°E
255 '''
256 d = _spherical_datum(datum or self.datum, name=self.name)
257 E = d.ellipsoid
258 h = self.h if h is None else Height(h=h)
260 x, y, z = self.x, self.y, self.z
261 # Kenneth Gade eqn 22
262 n = E.b / hypot_(x * E.a_b, y * E.a_b, z)
263 r = h + n * E.a2_b2
265 x *= r
266 y *= r
267 z *= n + h
269 if Cartesian is None:
270 r = self.Ecef(d).reverse(x, y, z, M=True)
271 else:
272 kwds = _xkwds(Cartesian_kwds, datum=d) # h=0
273 r = Cartesian(x, y, z, **kwds)
274 return self._xnamed(r)
276 @deprecated_method
277 def to2ll(self): # PYCHOK no cover
278 '''DEPRECATED, use property L{latlon}.
280 @return: A L{LatLon2Tuple}C{(lat, lon)}.
281 '''
282 return self.latlon
284 @deprecated_method
285 def to3llh(self, height=None): # PYCHOK no cover
286 '''DEPRECATED, use property C{latlonheight} or C{latlon.to3Tuple(B{height})}.
288 @kwarg height: Optional height, overriding this
289 n-vector's height (C{meter}).
291 @return: A L{LatLon3Tuple}C{(lat, lon, height)}.
293 @raise ValueError: Invalid B{C{height}}.
294 '''
295 return self.latlonheight if height in (None, self.h) else \
296 self.latlon.to3Tuple(height)
298 def toLatLon(self, height=None, LatLon=None, datum=None, **LatLon_kwds):
299 '''Convert this n-vector to an C{Nvector}-based geodetic point.
301 @kwarg height: Optional height, overriding this n-vector's
302 height (C{meter}).
303 @kwarg LatLon: Optional class to return the geodetic point
304 (C{LatLon}) or C{None}.
305 @kwarg datum: Optional, spherical datum (C{Datum}).
306 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
307 arguments, ignored if C{B{LatLon} is None}.
309 @return: The geodetic point (C{LatLon}) or if C{B{LatLon} is None},
310 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,
311 datum)} with C{C} and C{M} if available.
313 @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}}
314 argument.
316 @raise ValueError: Invalid B{C{height}}.
318 @example:
320 >>> v = Nvector(0.5, 0.5, 0.7071)
321 >>> p = v.toLatLon() # 45.0°N, 45.0°E
322 '''
323 d = _spherical_datum(datum or self.datum, name=self.name)
324 h = self.h if height is None else Height(height)
325 # use self.Cartesian(Cartesian=None) for better accuracy of the height
326 # than self.Ecef(d).forward(self.lat, self.lon, height=h, M=True)
327 if LatLon is None:
328 r = self.toCartesian(h=h, Cartesian=None, datum=d)
329 else:
330 kwds = _xkwds(LatLon_kwds, height=h, datum=d)
331 r = self._xnamed(LatLon(self.lat, self.lon, **kwds))
332 return r
334 def toStr(self, prec=5, fmt=Fmt.PAREN, sep=_COMMASPACE_): # PYCHOK expected
335 '''Return a string representation of this n-vector.
337 Height component is only included if non-zero.
339 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
340 @kwarg fmt: Enclosing backets format (C{str}).
341 @kwarg sep: Optional separator between components (C{str}).
343 @return: Comma-separated C{"(x, y, z [, h])"} enclosed in
344 B{C{fmt}} brackets (C{str}).
346 @example:
348 >>> Nvector(0.5, 0.5, 0.7071).toStr() # (0.5, 0.5, 0.7071)
349 >>> Nvector(0.5, 0.5, 0.7071, 1).toStr(-3) # (0.500, 0.500, 0.707, +1.00)
350 '''
351 t = Vector3d.toStr(self, prec=prec, fmt=NN, sep=sep)
352 if self.h:
353 t = sep.join((t, self.hStr()))
354 return (fmt % (t,)) if fmt else t
356 def toVector3d(self, norm=True):
357 '''Convert this n-vector to a 3-D vector, I{ignoring
358 the height}.
360 @kwarg norm: Normalize the 3-D vector (C{bool}).
362 @return: The (normalized) vector (L{Vector3d}).
363 '''
364 v = Vector3d.unit(self) if norm else self
365 return Vector3d(v.x, v.y, v.z, name=self.name)
367 @deprecated_method
368 def to4xyzh(self, h=None): # PYCHOK no cover
369 '''DEPRECATED, use property L{xyzh} or C{xyz.to4Tuple(B{h})}.
370 '''
371 return self.xyzh if h in (None, self.h) else Vector4Tuple(
372 self.x, self.y, self.z, h, name=self.name)
374 def unit(self, ll=None):
375 '''Normalize this n-vector to unit length.
377 @kwarg ll: Optional, original latlon (C{LatLon}).
379 @return: Normalized vector (C{Nvector}).
380 '''
381 return _xattrs(Vector3d.unit(self, ll=ll), _UNDER(_h_))
383 @Property_RO
384 def xyzh(self):
385 '''Get this n-vector's components (L{Vector4Tuple}C{(x, y, z, h)})
386 '''
387 return self.xyz.to4Tuple(self.h)
390NorthPole = NvectorBase(0, 0, +1, name=_NorthPole_) # North pole (C{Nvector})
391SouthPole = NvectorBase(0, 0, -1, name=_SouthPole_) # South pole (C{Nvector})
394class _N_vector_(NvectorBase):
395 '''(INTERNAL) Minimal, low-overhead C{n-vector}.
396 '''
397 def __init__(self, x, y, z, h=0, name=NN):
398 self._x, self._y, self._z = x, y, z
399 if h:
400 self._h = h
401 if name:
402 self.name = name
405class LatLonNvectorBase(LatLonBase):
406 '''(INTERNAL) Base class for n-vector-based ellipsoidal
407 and spherical C{LatLon} classes.
408 '''
410 def _update(self, updated, *attrs, **setters): # PYCHOK _Nv=None
411 '''(INTERNAL) Zap cached attributes if updated.
413 @see: C{ellipsoidalNvector.LatLon} and C{sphericalNvector.LatLon}
414 for the special case of B{C{_Nv}}.
415 '''
416 if updated:
417 _Nv = _xkwds_pop(setters, _Nv=None)
418 if _Nv is not None:
419 if _Nv._fromll is not None:
420 _Nv._fromll = None
421 self._Nv = None
422 LatLonBase._update(self, updated, *attrs, **setters)
424# def distanceTo(self, other, **kwds): # PYCHOK no cover
425# '''(INTERNAL) I{Must be overloaded}, see function C{notOverloaded}.
426# '''
427# _MODS.named.notOverloaded(self, other, **kwds)
429 def intersections2(self, radius1, other, radius2, **kwds): # PYCHOK expected
430 '''B{Not implemented}, throws a C{NotImplementedError} always.
431 '''
432 notImplemented(self, radius1, other, radius2, **kwds)
434 def others(self, *other, **name_other_up):
435 '''Refined class comparison.
437 @arg other: The other instance (C{LatLonNvectorBase}).
438 @kwarg name_other_up: Overriding C{name=other} and C{up=1}
439 keyword arguments.
441 @return: The B{C{other}} if compatible.
443 @raise TypeError: Incompatible B{C{other}} C{type}.
444 '''
445 if other:
446 other0 = other[0]
447 if isinstance(other0, (self.__class__, LatLonNvectorBase)): # XXX NvectorBase?
448 return other0
450 other, name, up = _xother3(self, other, **name_other_up)
451 if not isinstance(other, (self.__class__, LatLonNvectorBase)): # XXX NvectorBase?
452 LatLonBase.others(self, other, name=name, up=up + 1)
453 return other
455 def toNvector(self, Nvector=NvectorBase, **Nvector_kwds): # PYCHOK signature
456 '''Convert this point to C{Nvector} components, I{including
457 height}.
459 @kwarg Nvector_kwds: Optional, additional B{C{Nvector}} keyword
460 arguments, ignored if C{B{Nvector} is None}.
462 @return: An B{C{Nvector}} or a L{Vector4Tuple}C{(x, y, z, h)} if
463 B{C{Nvector}} is C{None}.
465 @raise TypeError: Invalid B{C{Nvector}} or B{C{Nvector_kwds}}
466 argument.
467 '''
468 return LatLonBase.toNvector(self, Nvector=Nvector, **Nvector_kwds)
470 def triangulate(self, bearing1, other, bearing2, height=None):
471 '''Locate a point given this and an other point and a bearing
472 at this and the other point.
474 @arg bearing1: Bearing at this point (compass C{degrees360}).
475 @arg other: The other point (C{LatLon}).
476 @arg bearing2: Bearing at the other point (compass C{degrees360}).
477 @kwarg height: Optional height at the triangulated point,
478 overriding the mean height (C{meter}).
480 @return: Triangulated point (C{LatLon}).
482 @raise TypeError: Invalid B{C{other}} point.
484 @raise Valuerror: Points coincide.
486 @example:
488 >>> p = LatLon("47°18.228'N","002°34.326'W") # Basse Castouillet
489 >>> q = LatLon("47°18.664'N","002°31.717'W") # Basse Hergo
490 >>> t = p.triangulate(7, q, 295) # 47.323667°N, 002.568501°W'
491 '''
492 return _triangulate(self, bearing1, self.others(other), bearing2,
493 height=height, LatLon=self.classof)
495 def trilaterate(self, distance1, point2, distance2, point3, distance3,
496 radius=R_M, height=None, useZ=False):
497 '''Locate a point at given distances from this and two other points.
499 @arg distance1: Distance to this point (C{meter}, same units
500 as B{C{radius}}).
501 @arg point2: Second reference point (C{LatLon}).
502 @arg distance2: Distance to point2 (C{meter}, same units as
503 B{C{radius}}).
504 @arg point3: Third reference point (C{LatLon}).
505 @arg distance3: Distance to point3 (C{meter}, same units as
506 B{C{radius}}).
507 @kwarg radius: Mean earth radius (C{meter}).
508 @kwarg height: Optional height at trilaterated point, overriding
509 the mean height (C{meter}, same units as B{C{radius}}).
510 @kwarg useZ: Include Z component iff non-NaN, non-zero (C{bool}).
512 @return: Trilaterated point (C{LatLon}).
514 @raise IntersectionError: No intersection, trilateration failed.
516 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
518 @raise ValueError: Some B{C{points}} coincide or invalid B{C{distance1}},
519 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.
521 @see: U{Trilateration<https://WikiPedia.org/wiki/Trilateration>},
522 Veness' JavaScript U{Trilateration<https://www.Movable-Type.co.UK/
523 scripts/latlong-vectors.html>} and method C{LatLon.trilaterate2}
524 of other, non-C{Nvector LatLon} classes.
525 '''
526 return _trilaterate(self, distance1,
527 self.others(point2=point2), distance2,
528 self.others(point3=point3), distance3,
529 radius=radius, height=height, useZ=useZ,
530 LatLon=self.classof)
532 def trilaterate5(self, distance1, point2, distance2, point3, distance3, # PYCHOK signature
533 area=False, eps=EPS1, radius=R_M, wrap=False):
534 '''B{Not implemented} for C{B{area}=True} or C{B{wrap}=True}
535 and falls back to method C{trilaterate} otherwise.
537 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
538 with a single trilaterated intersection C{minPoint I{is}
539 maxPoint}, C{min I{is} max} the nearest intersection
540 margin and count C{n = 1}.
542 @raise IntersectionError: No intersection, trilateration failed.
544 @raise NotImplementedError: Keyword argument C{B{area}=True} or
545 C{B{wrap}=True} not (yet) supported.
547 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
549 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}},
550 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.
551 '''
552 if area or wrap:
553 notImplemented(self, area=area, wrap=wrap)
555 t = _trilaterate(self, distance1, self.others(point2=point2), distance2,
556 self.others(point3=point3), distance3,
557 radius=radius, height=None, useZ=True,
558 LatLon=self.classof)
559 # ... and handle B{C{eps}} and C{IntersectionError}
560 # like function C{.latlonBase._trilaterate5}
561 d = self.distanceTo(t, radius=radius, wrap=wrap) # PYCHOK distanceTo
562 d = min(fabs(distance1 - d), fabs(distance2 - d), fabs(distance3 - d))
563 if d < eps: # min is max, minPoint is maxPoint
564 return Trilaterate5Tuple(d, t, d, t, 1) # n = 1
565 t = _SPACE_(_no_(_intersection_), Fmt.PAREN(min.__name__, Fmt.f(d, prec=3)))
566 raise IntersectionError(area=area, eps=eps, radius=radius, wrap=wrap, txt=t)
569def sumOf(nvectors, Vector=None, h=None, **Vector_kwds):
570 '''Return the vectorial sum of two or more n-vectors.
572 @arg nvectors: Vectors to be added (C{Nvector}[]).
573 @kwarg Vector: Optional class for the vectorial sum (C{Nvector})
574 or C{None}.
575 @kwarg h: Optional height, overriding the mean height (C{meter}).
576 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword
577 arguments, ignored if C{B{Vector} is None}.
579 @return: Vectorial sum (B{C{Vector}}) or a L{Vector4Tuple}C{(x, y,
580 z, h)} if B{C{Vector}} is C{None}.
582 @raise VectorError: No B{C{nvectors}}.
583 '''
584 n, nvectors = len2(nvectors)
585 if n < 1:
586 raise VectorError(nvectors=n, txt=MISSING)
588 if h is None:
589 h = fsum(v.h for v in nvectors) / float(n)
591 if Vector is None:
592 r = _sumOf(nvectors, Vector=Vector3Tuple).to4Tuple(h)
593 else:
594 r = _sumOf(nvectors, Vector=Vector, **_xkwds(Vector_kwds, h=h))
595 return r
598def _triangulate(point1, bearing1, point2, bearing2, height=None,
599 **LatLon_LatLon_kwds):
600 # (INTERNAL) Locate a point given two known points and initial
601 # bearings from those points, see C{LatLon.triangulate} above
603 def _gc(p, b, _i_):
604 n = p.toNvector()
605 de = NorthPole.cross(n, raiser=_pole_).unit() # east vector @ n
606 dn = n.cross(de) # north vector @ n
607 s, c = sincos2d(Bearing(b, name=_bearing_ + _i_))
608 dest = de.times(s)
609 dnct = dn.times(c)
610 d = dnct.plus(dest) # direction vector @ n
611 return n.cross(d) # great circle point + bearing
613 if point1.isequalTo(point2, EPS):
614 raise _ValueError(points=point2, txt=_coincident_)
616 gc1 = _gc(point1, bearing1, _1_) # great circle p1 + b1
617 gc2 = _gc(point2, bearing2, _2_) # great circle p2 + b2
619 n = gc1.cross(gc2, raiser=_points_) # n-vector of intersection point
621 h = point1._havg(point2) if height is None else Height(height)
622 kwds = _xkwds(LatLon_LatLon_kwds, height=h)
623 return n.toLatLon(**kwds) # Nvector(n.x, n.y, n.z).toLatLon(...)
626def _trilaterate(point1, distance1, point2, distance2, point3, distance3,
627 radius=R_M, height=None, useZ=False,
628 **LatLon_LatLon_kwds):
629 # (INTERNAL) Locate a point at given distances from
630 # three other points, see LatLon.triangulate above
632 def _nd2(p, d, r, _i_, *qs): # .toNvector and angular distance squared
633 for q in qs:
634 if p.isequalTo(q, EPS):
635 raise _ValueError(points=p, txt=_coincident_)
636 return p.toNvector(), (Scalar(d, name=_distance_ + _i_) / r)**2
638 r = Radius_(radius)
640 n1, r12 = _nd2(point1, distance1, r, _1_)
641 n2, r22 = _nd2(point2, distance2, r, _2_, point1)
642 n3, r32 = _nd2(point3, distance3, r, _3_, point1, point2)
644 # the following uses x,y coordinate system with origin at n1, x axis n1->n2
645 y = n3.minus(n1)
646 x = n2.minus(n1)
647 z = None
649 d = x.length # distance n1->n2
650 if d > EPS_2: # and y.length > EPS_2:
651 X = x.unit() # unit vector in x direction n1->n2
652 i = X.dot(y) # signed magnitude of x component of n1->n3
653 Y = y.minus(X.times(i)).unit() # unit vector in y direction
654 j = Y.dot(y) # signed magnitude of y component of n1->n3
655 if fabs(j) > EPS_2:
656 # courtesy of U{Carlos Freitas<https://GitHub.com/mrJean1/PyGeodesy/issues/33>}
657 x = fsum_(r12, -r22, d**2) / (d * _2_0) # n1->intersection x- and ...
658 y = fsum_(r12, -r32, i**2, j**2, x * i * _N_2_0) / (j * _2_0) # ... y-component
659 # courtesy of U{AleixDev<https://GitHub.com/mrJean1/PyGeodesy/issues/43>}
660 z = fsum_(max(r12, r22, r32), -(x**2), -(y**2)) # XXX not just r12!
661 if z > EPS:
662 n = n1.plus(X.times(x)).plus(Y.times(y))
663 if useZ: # include Z component
664 Z = X.cross(Y) # unit vector perpendicular to plane
665 n = n.plus(Z.times(sqrt(z)))
666 if height is None:
667 h = fidw((point1.height, point2.height, point3.height),
668 map1(fabs, distance1, distance2, distance3))
669 else:
670 h = Height(height)
671 kwds = _xkwds(LatLon_LatLon_kwds, height=h)
672 return n.toLatLon(**kwds) # Nvector(n.x, n.y, n.z).toLatLon(...)
674 # no intersection, d < EPS_2 or fabs(j) < EPS_2 or z < EPS
675 t = _SPACE_(_no_, _intersection_, NN)
676 raise IntersectionError(point1=point1, distance1=distance1,
677 point2=point2, distance2=distance2,
678 point3=point3, distance3=distance3,
679 txt=unstr(t, z=z, useZ=useZ))
682__all__ += _ALL_DOCS(LatLonNvectorBase, NvectorBase, sumOf) # classes
684# **) MIT License
685#
686# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
687#
688# Permission is hereby granted, free of charge, to any person obtaining a
689# copy of this software and associated documentation files (the "Software"),
690# to deal in the Software without restriction, including without limitation
691# the rights to use, copy, modify, merge, publish, distribute, sublicense,
692# and/or sell copies of the Software, and to permit persons to whom the
693# Software is furnished to do so, subject to the following conditions:
694#
695# The above copyright notice and this permission notice shall be included
696# in all copies or substantial portions of the Software.
697#
698# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
699# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
700# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
701# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
702# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
703# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
704# OTHER DEALINGS IN THE SOFTWARE.