Coverage for pygeodesy/nvectorBase.py: 96%
249 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-01 11:43 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-01 11:43 -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'''
13# from pygeodesy.basics import map1 # from .namedTuples
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_pop2
18from pygeodesy.fmath import fdot, fidw, hypot_ # PYCHOK fdot shared
19from pygeodesy.fsums import Fsum, fsumf_
20from pygeodesy.formy import _isequalTo, n_xyz2latlon, n_xyz2philam, \
21 _spherical_datum
22# from pygeodesy.internals import _under # from .named
23from pygeodesy.interns import NN, _1_, _2_, _3_, _bearing_, _coincident_, \
24 _COMMASPACE_, _distance_, _h_, _insufficient_, \
25 _intersection_, _no_, _NorthPole_, _point_, \
26 _pole_, _SPACE_, _SouthPole_
27from pygeodesy.latlonBase import LatLonBase, _ALL_DOCS, _ALL_LAZY, _MODS
28# from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS # from .latlonBase
29from pygeodesy.named import _xother3, _under
30from pygeodesy.namedTuples import Trilaterate5Tuple, Vector3Tuple, \
31 Vector4Tuple, map1
32from pygeodesy.props import deprecated_method, Property_RO, property_doc_, \
33 property_RO, _update_all
34from pygeodesy.streprs import Fmt, hstr, unstr, _xattrs
35from pygeodesy.units import Bearing, Height, Radius_, Scalar
36from pygeodesy.utily import sincos2d, _unrollon, _unrollon3
37from pygeodesy.vector3d import Vector3d, _xyzhdn3
39from math import fabs, sqrt
41__all__ = _ALL_LAZY.nvectorBase
42__version__ = '24.05.31'
45class NvectorBase(Vector3d): # XXX kept private
46 '''Base class for ellipsoidal and spherical C{Nvector}s.
47 '''
48 _datum = None # L{Datum}, overriden
49 _h = Height(h=0) # height (C{meter})
50 _H = NN # height prefix (C{str}), '↑' in JS version
52 def __init__(self, x_xyz, y=None, z=None, h=0, ll=None, datum=None, **name):
53 '''New n-vector normal to the earth's surface.
55 @arg x_xyz: X component of vector (C{scalar}) or (3-D) vector
56 (C{Nvector}, L{Vector3d}, L{Vector3Tuple} or
57 L{Vector4Tuple}).
58 @kwarg y: Y 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 z: Z component of vector (C{scalar}), ignored if B{C{x_xyz}}
61 is not C{scalar}, otherwise same units as B{C{x_xyz}}.
62 @kwarg h: Optional height above surface (C{meter}).
63 @kwarg ll: Optional, original latlon (C{LatLon}).
64 @kwarg datum: Optional, I{pass-thru} datum (L{Datum}).
65 @kwarg name: Optional C{B{name}=NN} (C{str}).
67 @raise TypeError: Non-scalar B{C{x}}, B{C{y}} or B{C{z}}
68 coordinate or B{C{x}} not an C{Nvector},
69 L{Vector3Tuple} or L{Vector4Tuple} or
70 invalid B{C{datum}}.
71 '''
72 h, d, n = _xyzhdn3(x_xyz, h, datum, ll, **name)
73 Vector3d.__init__(self, x_xyz, y=y, z=z, ll=ll, name=n)
74 if h:
75 self.h = h
76 if d is not None:
77 self._datum = _spherical_datum(d, name=self.name) # pass-thru
79 @Property_RO
80 def datum(self):
81 '''Get the I{pass-thru} datum (C{Datum}) or C{None}.
82 '''
83 return self._datum
85 @property_RO
86 def Ecef(self):
87 '''Get the ECEF I{class} (L{EcefKarney}), I{once}.
88 '''
89 NvectorBase.Ecef = E = _MODS.ecef.EcefKarney # overwrite property_RO
90 return E
92 @property_RO
93 def ellipsoidalNvector(self):
94 '''Get the C{Nvector type} iff ellipsoidal, overloaded in L{pygeodesy.ellipsoidalNvector.Nvector}.
95 '''
96 return False
98 @property_doc_(''' the height above surface (C{meter}).''')
99 def h(self):
100 '''Get the height above surface (C{meter}).
101 '''
102 return self._h
104 @h.setter # PYCHOK setter!
105 def h(self, h):
106 '''Set the height above surface (C{meter}).
108 @raise TypeError: If B{C{h}} invalid.
110 @raise VectorError: If B{C{h}} invalid.
111 '''
112 h = Height(h=h, Error=VectorError)
113 if self._h != h:
114 _update_all(self)
115 self._h = h
117 @property_doc_(''' the height prefix (C{str}).''')
118 def H(self):
119 '''Get the height prefix (C{str}).
120 '''
121 return self._H
123 @H.setter # PYCHOK setter!
124 def H(self, H):
125 '''Set the height prefix (C{str}).
126 '''
127 self._H = str(H) if H else NN
129 def hStr(self, prec=-2, m=NN):
130 '''Return a string for the height B{C{h}}.
132 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
133 @kwarg m: Optional unit of the height (C{str}).
135 @see: Function L{pygeodesy.hstr}.
136 '''
137 return NN(self.H, hstr(self.h, prec=prec, m=m))
139 @Property_RO
140 def isEllipsoidal(self):
141 '''Check whether this n-vector is ellipsoidal (C{bool} or C{None} if unknown).
142 '''
143 return self.datum.isEllipsoidal if self.datum else None
145 @Property_RO
146 def isSpherical(self):
147 '''Check whether this n-vector is spherical (C{bool} or C{None} if unknown).
148 '''
149 return self.datum.isSpherical if self.datum else None
151 @Property_RO
152 def lam(self):
153 '''Get the (geodetic) longitude in C{radians} (C{float}).
154 '''
155 return self.philam.lam
157 @Property_RO
158 def lat(self):
159 '''Get the (geodetic) latitude in C{degrees} (C{float}).
160 '''
161 return self.latlon.lat
163 @Property_RO
164 def latlon(self):
165 '''Get the (geodetic) lat-, longitude in C{degrees} (L{LatLon2Tuple}C{(lat, lon)}).
166 '''
167 return n_xyz2latlon(self.x, self.y, self.z, name=self.name)
169 @Property_RO
170 def latlonheight(self):
171 '''Get the (geodetic) lat-, longitude in C{degrees} and height (L{LatLon3Tuple}C{(lat, lon, height)}).
172 '''
173 return self.latlon.to3Tuple(self.h)
175 @Property_RO
176 def latlonheightdatum(self):
177 '''Get the lat-, longitude in C{degrees} with height and datum (L{LatLon4Tuple}C{(lat, lon, height, datum)}).
178 '''
179 return self.latlonheight.to4Tuple(self.datum)
181 @Property_RO
182 def lon(self):
183 '''Get the (geodetic) longitude in C{degrees} (C{float}).
184 '''
185 return self.latlon.lon
187 @Property_RO
188 def phi(self):
189 '''Get the (geodetic) latitude in C{radians} (C{float}).
190 '''
191 return self.philam.phi
193 @Property_RO
194 def philam(self):
195 '''Get the (geodetic) lat-, longitude in C{radians} (L{PhiLam2Tuple}C{(phi, lam)}).
196 '''
197 return n_xyz2philam(self.x, self.y, self.z, name=self.name)
199 @Property_RO
200 def philamheight(self):
201 '''Get the (geodetic) lat-, longitude in C{radians} and height (L{PhiLam3Tuple}C{(phi, lam, height)}).
202 '''
203 return self.philam.to3Tuple(self.h)
205 @Property_RO
206 def philamheightdatum(self):
207 '''Get the lat-, longitude in C{radians} with height and datum (L{PhiLam4Tuple}C{(phi, lam, height, datum)}).
208 '''
209 return self.philamheight.to4Tuple(self.datum)
211 @property_RO
212 def sphericalNvector(self):
213 '''Get the C{Nvector type} iff spherical, overloaded in L{pygeodesy.sphericalNvector.Nvector}.
214 '''
215 return False
217 @deprecated_method
218 def to2ab(self): # PYCHOK no cover
219 '''DEPRECATED, use property L{philam}.
221 @return: A L{PhiLam2Tuple}C{(phi, lam)}.
222 '''
223 return self.philam
225 @deprecated_method
226 def to3abh(self, height=None): # PYCHOK no cover
227 '''DEPRECATED, use property L{philamheight} or C{philam.to3Tuple(B{height})}.
229 @kwarg height: Optional height, overriding this
230 n-vector's height (C{meter}).
232 @return: A L{PhiLam3Tuple}C{(phi, lam, height)}.
234 @raise ValueError: Invalid B{C{height}}.
235 '''
236 return self.philamheight if height in (None, self.h) else \
237 self.philam.to3Tuple(height)
239 def toCartesian(self, h=None, Cartesian=None, datum=None, **Cartesian_kwds):
240 '''Convert this n-vector to C{Nvector}-based cartesian (ECEF) coordinates.
242 @kwarg h: Optional height, overriding this n-vector's height (C{meter}).
243 @kwarg Cartesian: Optional class to return the (ECEF) coordinates
244 (C{Cartesian}).
245 @kwarg datum: Optional datum (C{Datum}), overriding this datum.
246 @kwarg Cartesian_kwds: Optional, additional B{C{Cartesian}} keyword
247 arguments, ignored if C{B{Cartesian} is None}.
249 @return: The cartesian (ECEF) coordinates (B{C{Cartesian}}) or
250 if C{B{Cartesian} is None}, an L{Ecef9Tuple}C{(x, y, z,
251 lat, lon, height, C, M, datum)} with C{C} and C{M} if
252 available.
254 @raise TypeError: Invalid B{C{Cartesian}} or B{C{Cartesian_kwds}}
255 argument.
257 @raise ValueError: Invalid B{C{h}}.
258 '''
259 D = _spherical_datum(datum or self.datum, name=self.name)
260 E = D.ellipsoid
261 h = self.h if h is None else Height(h)
263 x, y, z = self.x, self.y, self.z
264 # Kenneth Gade eqn 22
265 n = E.b / hypot_(x * E.a_b, y * E.a_b, z)
266 r = h + n * E.a2_b2
268 x *= r
269 y *= r
270 z *= h + n
272 if Cartesian is None:
273 r = self.Ecef(D).reverse(x, y, z, M=True)
274 else:
275 kwds = _xkwds(Cartesian_kwds, datum=D) # h=0
276 r = Cartesian(x, y, z, **kwds)
277 return self._xnamed(r)
279 @deprecated_method
280 def to2ll(self): # PYCHOK no cover
281 '''DEPRECATED, use property L{latlon}.
283 @return: A L{LatLon2Tuple}C{(lat, lon)}.
284 '''
285 return self.latlon
287 @deprecated_method
288 def to3llh(self, height=None): # PYCHOK no cover
289 '''DEPRECATED, use property C{latlonheight} or C{latlon.to3Tuple(B{height})}.
291 @kwarg height: Optional height, overriding this
292 n-vector's height (C{meter}).
294 @return: A L{LatLon3Tuple}C{(lat, lon, height)}.
296 @raise ValueError: Invalid B{C{height}}.
297 '''
298 return self.latlonheight if height in (None, self.h) else \
299 self.latlon.to3Tuple(height)
301 def toLatLon(self, height=None, LatLon=None, datum=None, **LatLon_kwds):
302 '''Convert this n-vector to an C{Nvector}-based geodetic point.
304 @kwarg height: Optional height, overriding this n-vector's
305 height (C{meter}).
306 @kwarg LatLon: Optional class to return the geodetic point
307 (C{LatLon}) or C{None}.
308 @kwarg datum: Optional, spherical datum (C{Datum}).
309 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
310 arguments, ignored if C{B{LatLon} is None}.
312 @return: The geodetic point (C{LatLon}) or if C{B{LatLon} is None},
313 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M,
314 datum)} with C{C} and C{M} if available.
316 @raise TypeError: Invalid B{C{LatLon}} or B{C{LatLon_kwds}}
317 argument.
319 @raise ValueError: Invalid B{C{height}}.
320 '''
321 d = _spherical_datum(datum or self.datum, name=self.name)
322 h = self.h if height is None else Height(height)
323 # use self.Cartesian(Cartesian=None) for better accuracy of the height
324 # than self.Ecef(d).forward(self.lat, self.lon, height=h, M=True)
325 if LatLon is None:
326 r = self.toCartesian(h=h, Cartesian=None, datum=d)
327 else:
328 kwds = _xkwds(LatLon_kwds, height=h, datum=d)
329 r = self._xnamed(LatLon(self.lat, self.lon, **kwds))
330 return r
332 def toStr(self, prec=5, fmt=Fmt.PAREN, sep=_COMMASPACE_): # PYCHOK expected
333 '''Return a string representation of this n-vector.
335 Height component is only included if non-zero.
337 @kwarg prec: Number of (decimal) digits, unstripped (C{int}).
338 @kwarg fmt: Enclosing backets format (C{str}).
339 @kwarg sep: Optional separator between components (C{str}).
341 @return: Comma-separated C{"(x, y, z [, h])"} enclosed in
342 B{C{fmt}} brackets (C{str}).
343 '''
344 t = Vector3d.toStr(self, prec=prec, fmt=NN, sep=sep)
345 if self.h:
346 t = sep.join((t, self.hStr()))
347 return (fmt % (t,)) if fmt else t
349 def toVector3d(self, norm=True):
350 '''Convert this n-vector to a 3-D vector, I{ignoring height}.
352 @kwarg norm: Normalize the 3-D vector (C{bool}).
354 @return: The (normalized) vector (L{Vector3d}).
355 '''
356 v = Vector3d.unit(self) if norm else self
357 return Vector3d(v.x, v.y, v.z, name=self.name)
359 @deprecated_method
360 def to4xyzh(self, h=None): # PYCHOK no cover
361 '''DEPRECATED, use property L{xyzh} or C{xyz.to4Tuple(B{h})}.'''
362 return self.xyzh if h in (None, self.h) else Vector4Tuple(
363 self.x, self.y, self.z, h, name=self.name)
365 def unit(self, ll=None):
366 '''Normalize this n-vector to unit length.
368 @kwarg ll: Optional, original latlon (C{LatLon}).
370 @return: Normalized vector (C{Nvector}).
371 '''
372 return _xattrs(Vector3d.unit(self, ll=ll), _under(_h_))
374 @Property_RO
375 def xyzh(self):
376 '''Get this n-vector's components (L{Vector4Tuple}C{(x, y, z, h)})
377 '''
378 return self.xyz.to4Tuple(self.h)
381NorthPole = NvectorBase(0, 0, +1, name=_NorthPole_) # North pole (C{Nvector})
382SouthPole = NvectorBase(0, 0, -1, name=_SouthPole_) # South pole (C{Nvector})
385class _N_vector_(NvectorBase):
386 '''(INTERNAL) Minimal, low-overhead C{n-vector}.
387 '''
388 def __init__(self, x, y, z, h=0, **name):
389 self._x, self._y, self._z = x, y, z
390 if h:
391 self._h = h
392 if name:
393 self.name = name
396class LatLonNvectorBase(LatLonBase):
397 '''(INTERNAL) Base class for n-vector-based ellipsoidal and
398 spherical C{LatLon} classes.
399 '''
401 def _update(self, updated, *attrs, **setters): # PYCHOK _Nv=None
402 '''(INTERNAL) Zap cached attributes if updated.
404 @see: C{ellipsoidalNvector.LatLon} and C{sphericalNvector.LatLon}
405 for the special case of B{C{_Nv}}.
406 '''
407 if updated:
408 _Nv, setters = _xkwds_pop2(setters, _Nv=None)
409 if _Nv is not None:
410 if _Nv._fromll is not None:
411 _Nv._fromll = None
412 self._Nv = None
413 LatLonBase._update(self, updated, *attrs, **setters)
415# def distanceTo(self, other, **kwds): # PYCHOK no cover
416# '''I{Must be overloaded}.'''
417# self._notOverloaded(other, **kwds)
419 def intersections2(self, radius1, other, radius2, **kwds): # PYCHOK expected
420 '''B{Not implemented}, throws a C{NotImplementedError} always.'''
421 self._notImplemented(radius1, other, radius2, **kwds)
423 def others(self, *other, **name_other_up):
424 '''Refined class comparison.
426 @arg other: The other instance (C{LatLonNvectorBase}).
427 @kwarg name_other_up: Overriding C{name=other} and C{up=1}
428 keyword arguments.
430 @return: The B{C{other}} if compatible.
432 @raise TypeError: Incompatible B{C{other}} C{type}.
433 '''
434 if other:
435 other0 = other[0]
436 if isinstance(other0, (self.__class__, LatLonNvectorBase)): # XXX NvectorBase?
437 return other0
439 other, name, up = _xother3(self, other, **name_other_up)
440 if not isinstance(other, (self.__class__, LatLonNvectorBase)): # XXX NvectorBase?
441 LatLonBase.others(self, other, name=name, up=up + 1)
442 return other
444 def toNvector(self, **Nvector_and_kwds): # PYCHOK signature
445 '''Convert this point to C{Nvector} components, I{including height}.
447 @kwarg Nvector_and_kwds: Optional C{Nvector} class and C{Nvector} keyword arguments,
448 Specify C{B{Nvector}=...} to override this C{Nvector} class
449 or use C{B{Nvector}=None}.
451 @return: An C{Nvector} or if C{Nvector is None}, a L{Vector4Tuple}C{(x, y, z, h)}.
453 @raise TypeError: Invalid C{Nvector} or other B{C{Nvector_and_kwds}} item.
454 '''
455 return LatLonBase.toNvector(self, **_xkwds(Nvector_and_kwds, Nvector=NvectorBase))
457 def triangulate(self, bearing1, other, bearing2, height=None, wrap=False): # PYCHOK signature
458 '''Locate a point given this, an other point and the (initial) bearing
459 from this and the other point.
461 @arg bearing1: Bearing at this point (compass C{degrees360}).
462 @arg other: The other point (C{LatLon}).
463 @arg bearing2: Bearing at the other point (compass C{degrees360}).
464 @kwarg height: Optional height at the triangulated point, overriding
465 the mean height (C{meter}).
466 @kwarg wrap: If C{True}, use this and the B{C{other}} point
467 I{normalized} (C{bool}).
469 @return: Triangulated point (C{LatLon}).
471 @raise TypeError: Invalid B{C{other}} point.
473 @raise Valuerror: Points coincide.
474 '''
475 return _triangulate(self, bearing1, self.others(other), bearing2,
476 height=height, wrap=wrap, LatLon=self.classof)
478 def trilaterate(self, distance1, point2, distance2, point3, distance3,
479 radius=R_M, height=None, useZ=False, wrap=False):
480 '''Locate a point at given distances from this and two other points.
482 @arg distance1: Distance to this point (C{meter}, same units
483 as B{C{radius}}).
484 @arg point2: Second reference point (C{LatLon}).
485 @arg distance2: Distance to point2 (C{meter}, same units as
486 B{C{radius}}).
487 @arg point3: Third reference point (C{LatLon}).
488 @arg distance3: Distance to point3 (C{meter}, same units as
489 B{C{radius}}).
490 @kwarg radius: Mean earth radius (C{meter}).
491 @kwarg height: Optional height at trilaterated point, overriding
492 the mean height (C{meter}, same units as B{C{radius}}).
493 @kwarg useZ: Include Z component iff non-NaN, non-zero (C{bool}).
494 @kwarg wrap: If C{True}, use this, B{C{point2}} and B{C{point3}}
495 I{normalized} (C{bool}).
497 @return: Trilaterated point (C{LatLon}).
499 @raise IntersectionError: No intersection, trilateration failed.
501 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
503 @raise ValueError: Some B{C{points}} coincide or invalid B{C{distance1}},
504 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.
506 @see: U{Trilateration<https://WikiPedia.org/wiki/Trilateration>},
507 Veness' JavaScript U{Trilateration<https://www.Movable-Type.co.UK/
508 scripts/latlong-vectors.html>} and method C{LatLon.trilaterate5}
509 of other, non-C{Nvector LatLon} classes.
510 '''
511 return _trilaterate(self, distance1, self.others(point2=point2), distance2,
512 self.others(point3=point3), distance3,
513 radius=radius, height=height, useZ=useZ,
514 wrap=wrap, LatLon=self.classof)
516 def trilaterate5(self, distance1, point2, distance2, point3, distance3, # PYCHOK signature
517 area=False, eps=EPS1, radius=R_M, wrap=False):
518 '''B{Not implemented} for C{B{area}=True} and falls back to method
519 C{trilaterate} otherwise.
521 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
522 with a single trilaterated intersection C{minPoint I{is}
523 maxPoint}, C{min I{is} max} the nearest intersection
524 margin and count C{n = 1}.
526 @raise NotImplementedError: Keyword argument C{B{area}=True} not
527 (yet) supported.
529 @see: Method L{trilaterate} for other and more details.
530 '''
531 if area:
532 self._notImplemented(area=area)
534 t = _trilaterate(self, distance1, self.others(point2=point2), distance2,
535 self.others(point3=point3), distance3,
536 radius=radius, useZ=True, wrap=wrap,
537 LatLon=self.classof)
538 # ... and handle B{C{eps}} and C{IntersectionError}
539 # like function C{.latlonBase._trilaterate5}
540 d = self.distanceTo(t, radius=radius, wrap=wrap) # PYCHOK distanceTo
541 d = min(fabs(distance1 - d), fabs(distance2 - d), fabs(distance3 - d))
542 if d < eps: # min is max, minPoint is maxPoint
543 return Trilaterate5Tuple(d, t, d, t, 1) # n = 1
544 t = _SPACE_(_no_(_intersection_), Fmt.PAREN(min.__name__, Fmt.f(d, prec=3)))
545 raise IntersectionError(area=area, eps=eps, radius=radius, wrap=wrap, txt=t)
548def _nsumOf(nvs, h_None, Vector, Vector_kwds): # .sphericalNvector, .vector3d
549 '''(INTERNAL) Separated to allow callers to embellish exceptions.
550 '''
551 X, Y, Z, n = Fsum(), Fsum(), Fsum(), 0
552 H = Fsum() if h_None is None else n
553 for n, v in enumerate(nvs or ()): # one pass
554 X += v.x
555 Y += v.y
556 Z += v.z
557 H += v.h
558 if n < 1:
559 raise ValueError(_SPACE_(Fmt.PARENSPACED(len=n), _insufficient_))
561 x, y, z = map1(float, X, Y, Z)
562 h = H.fover(n) if h_None is None else h_None
563 return Vector3Tuple(x, y, z).to4Tuple(h) if Vector is None else \
564 Vector(x, y, z, **_xkwds(Vector_kwds, h=h))
567def sumOf(nvectors, Vector=None, h=None, **Vector_kwds):
568 '''Return the I{vectorial} sum of two or more n-vectors.
570 @arg nvectors: Vectors to be added (C{Nvector}[]).
571 @kwarg Vector: Optional class for the vectorial sum (C{Nvector})
572 or C{None}.
573 @kwarg h: Optional height, overriding the mean height (C{meter}).
574 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword
575 arguments, ignored if C{B{Vector} is None}.
577 @return: Vectorial sum (B{C{Vector}}) or a L{Vector4Tuple}C{(x, y,
578 z, h)} if B{C{Vector}} is C{None}.
580 @raise VectorError: No B{C{nvectors}}.
581 '''
582 try:
583 return _nsumOf(nvectors, h, Vector, Vector_kwds)
584 except (TypeError, ValueError) as x:
585 raise VectorError(nvectors=nvectors, Vector=Vector, cause=x)
588def _triangulate(point1, bearing1, point2, bearing2, height=None,
589 wrap=False, **LatLon_and_kwds):
590 # (INTERNAL) Locate a point given two known points and initial
591 # bearings from those points, see C{LatLon.triangulate} above
593 def _gc(p, b, _i_):
594 n = p.toNvector()
595 de = NorthPole.cross(n, raiser=_pole_).unit() # east vector @ n
596 dn = n.cross(de) # north vector @ n
597 s, c = sincos2d(Bearing(b, name=_bearing_ + _i_))
598 dest = de.times(s)
599 dnct = dn.times(c)
600 d = dnct.plus(dest) # direction vector @ n
601 return n.cross(d) # great circle point + bearing
603 if wrap:
604 point2 = _unrollon(point1, point2, wrap=wrap)
605 if _isequalTo(point1, point2, eps=EPS):
606 raise _ValueError(points=point2, wrap=wrap, txt=_coincident_)
608 gc1 = _gc(point1, bearing1, _1_) # great circle p1 + b1
609 gc2 = _gc(point2, bearing2, _2_) # great circle p2 + b2
611 n = gc1.cross(gc2, raiser=_point_) # n-vector of intersection point
612 h = point1._havg(point2, h=height)
613 kwds = _xkwds(LatLon_and_kwds, height=h)
614 return n.toLatLon(**kwds) # Nvector(n.x, n.y, n.z).toLatLon(...)
617def _trilaterate(point1, distance1, point2, distance2, point3, distance3,
618 radius=R_M, height=None, useZ=False,
619 wrap=False, **LatLon_and_kwds):
620 # (INTERNAL) Locate a point at given distances from
621 # three other points, see LatLon.triangulate above
623 def _nr2(p, d, r, _i_, *qs): # .toNvector and angular distance squared
624 for q in qs:
625 if _isequalTo(p, q, eps=EPS):
626 raise _ValueError(points=p, txt=_coincident_)
627 return p.toNvector(), (Scalar(d, name=_distance_ + _i_) / r)**2
629 p1, r = point1, Radius_(radius)
630 p2, p3, _ = _unrollon3(p1, point2, point3, wrap)
632 n1, r12 = _nr2(p1, distance1, r, _1_)
633 n2, r22 = _nr2(p2, distance2, r, _2_, p1)
634 n3, r32 = _nr2(p3, distance3, r, _3_, p1, p2)
636 # the following uses x,y coordinate system with origin at n1, x axis n1->n2
637 y = n3.minus(n1)
638 x = n2.minus(n1)
639 z = None
641 d = x.length # distance n1->n2
642 if d > EPS_2: # and y.length > EPS_2:
643 X = x.unit() # unit vector in x direction n1->n2
644 i = X.dot(y) # signed magnitude of x component of n1->n3
645 Y = y.minus(X.times(i)).unit() # unit vector in y direction
646 j = Y.dot(y) # signed magnitude of y component of n1->n3
647 if fabs(j) > EPS_2:
648 # courtesy of U{Carlos Freitas<https://GitHub.com/mrJean1/PyGeodesy/issues/33>}
649 x = fsumf_(r12, -r22, d**2) / (d * _2_0) # n1->intersection x- and ...
650 y = fsumf_(r12, -r32, i**2, j**2, x * i * _N_2_0) / (j * _2_0) # ... y-component
651 # courtesy of U{AleixDev<https://GitHub.com/mrJean1/PyGeodesy/issues/43>}
652 z = fsumf_(max(r12, r22, r32), -(x**2), -(y**2)) # XXX not just r12!
653 if z > EPS:
654 n = n1.plus(X.times(x)).plus(Y.times(y))
655 if useZ: # include Z component
656 Z = X.cross(Y) # unit vector perpendicular to plane
657 n = n.plus(Z.times(sqrt(z)))
658 if height is None:
659 h = fidw((point1.height, point2.height, point3.height),
660 map1(fabs, distance1, distance2, distance3))
661 else:
662 h = Height(height)
663 kwds = _xkwds(LatLon_and_kwds, height=h)
664 return n.toLatLon(**kwds) # Nvector(n.x, n.y, n.z).toLatLon(...)
666 # no intersection, d < EPS_2 or fabs(j) < EPS_2 or z < EPS
667 t = _SPACE_(_no_, _intersection_, NN)
668 raise IntersectionError(point1=point1, distance1=distance1,
669 point2=point2, distance2=distance2,
670 point3=point3, distance3=distance3,
671 txt=unstr(t, z=z, useZ=useZ, wrap=wrap))
674__all__ += _ALL_DOCS(LatLonNvectorBase, NvectorBase, sumOf) # classes
676# **) MIT License
677#
678# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
679#
680# Permission is hereby granted, free of charge, to any person obtaining a
681# copy of this software and associated documentation files (the "Software"),
682# to deal in the Software without restriction, including without limitation
683# the rights to use, copy, modify, merge, publish, distribute, sublicense,
684# and/or sell copies of the Software, and to permit persons to whom the
685# Software is furnished to do so, subject to the following conditions:
686#
687# The above copyright notice and this permission notice shall be included
688# in all copies or substantial portions of the Software.
689#
690# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
691# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
692# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
693# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
694# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
695# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
696# OTHER DEALINGS IN THE SOFTWARE.