Coverage for pygeodesy/sphericalNvector.py: 92%
320 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-07-12 13:40 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-07-12 13:40 -0400
2# -*- coding: utf-8 -*-
4u'''Spherical, C{N-vector}-based geodesy.
6N-vector-based classes geodetic (lat-/longitude) L{LatLon}, geocentric
7(ECEF) L{Cartesian} and L{Nvector} and functions L{areaOf}, L{intersection},
8L{meanOf}, L{nearestOn3}, L{perimeterOf}, L{sumOf}, L{triangulate} and
9L{trilaterate}, I{all spherical}.
11Pure Python implementation of n-vector-based spherical geodetic (lat-/longitude)
12methods, transcoded from JavaScript originals by I{(C) Chris Veness 2011-2016},
13published under the same MIT Licence**. See U{Vector-based geodesy
14<https://www.Movable-Type.co.UK/scripts/latlong-vectors.html>} and
15U{Module latlon-nvector-spherical
16<https://www.Movable-Type.co.UK/scripts/geodesy/docs/module-latlon-nvector-spherical.html>}.
18Tools for working with points and lines on (a spherical model of) the
19earth’s surface using using n-vectors rather than the more common
20spherical trigonometry. N-vectors make many calculations much simpler,
21and easier to follow, compared with the trigonometric equivalents.
23Based on Kenneth Gade’s U{‘Non-singular Horizontal Position Representation’
24<https://www.NavLab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf>},
25The Journal of Navigation (2010), vol 63, nr 3, pp 395-417.
27Note that the formulations below take x => 0°N,0°E, y => 0°N,90°E and
28z => 90°N while Gade uses x => 90°N, y => 0°N,90°E, z => 0°N,0°E.
30Also note that on a spherical earth model, an n-vector is equivalent
31to a normalised version of an (ECEF) cartesian coordinate.
32'''
33# make sure int/int division yields float quosient, see .basics
34from __future__ import division as _; del _ # PYCHOK semicolon
36from pygeodesy.basics import isscalar, _xinstanceof
37from pygeodesy.constants import EPS, EPS0, PI, PI2, PI_2, R_M, \
38 _0_0, _0_5, _1_0
39# from pygeodesy.datums import Datums # from .sphericalBase
40from pygeodesy.errors import PointsError, VectorError, _xError, _xkwds
41from pygeodesy.fmath import fmean, fsum
42# from pygeodesy.fsums import fsum # from .fmath
43from pygeodesy.interns import _composite_, _end_, _Nv00_, _other_, \
44 _point_, _pole_
45from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER
46# from pygeodesy.named import notImplemented # from .nvectorBase
47# from pygeodesy.namedTuples import NearestOn3Tuple # from .points
48from pygeodesy.nvectorBase import LatLonNvectorBase, NorthPole, \
49 notImplemented, NvectorBase, _nsumOf, \
50 _triangulate, _trilaterate
51from pygeodesy.points import NearestOn3Tuple, ispolar # PYCHOK exported
52from pygeodesy.props import deprecated_function, deprecated_method
53from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, \
54 _intersecant2, LatLonSphericalBase, \
55 Datums
56from pygeodesy.units import Bearing, Bearing_, Radius, Scalar
57from pygeodesy.utily import atan2, degrees360, fabs, sincos2, sincos2_, \
58 sincos2d, _unrollon, _Wrap
60# from math import atan2, fabs # from utily
62__all__ = _ALL_LAZY.sphericalNvector
63__version__ = '23.06.12'
65_lines_ = 'lines'
68class Cartesian(CartesianSphericalBase):
69 '''Extended to convert geocentric, L{Cartesian} points to
70 L{Nvector} and n-vector-based, spherical L{LatLon}.
71 '''
73 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon
74 '''Convert this cartesian to an C{Nvector}-based geodetic point.
76 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
77 arguments, like C{datum}. Use C{B{LatLon}=...}
78 to override this L{LatLon} class or specify
79 C{B{LatLon}=None}.
81 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set
82 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
83 C, M, datum)} with C{C} and C{M} if available.
85 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
86 '''
87 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
88 return CartesianSphericalBase.toLatLon(self, **kwds)
90 def toNvector(self, **Nvector_and_kwds): # PYCHOK Datums.WGS84
91 '''Convert this cartesian to L{Nvector} components, I{including height}.
93 @kwarg Nvector_and_kwds: Optional L{Nvector} and L{Nvector} keyword
94 arguments, like C{datum}. Use C{B{Nvector}=...}
95 to override this L{Nvector} class or specify
96 C{B{Nvector}=None}.
98 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}}
99 is set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}
101 @raise TypeError: Invalid B{C{Nvector_and_kwds}} argument.
102 '''
103 # ll = CartesianBase.toLatLon(self, LatLon=LatLon,
104 # datum=datum or self.datum)
105 # kwds = _xkwds(kwds, Nvector=Nvector)
106 # return ll.toNvector(**kwds)
107 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector, datum=self.datum)
108 return CartesianSphericalBase.toNvector(self, **kwds)
111class LatLon(LatLonNvectorBase, LatLonSphericalBase):
112 '''New n-vector based point on a spherical earth model.
114 Tools for working with points, lines and paths on (a spherical
115 model of) the earth's surface using vector-based methods.
117 @example:
119 >>> from sphericalNvector import LatLon
120 >>> p = LatLon(52.205, 0.119)
121 '''
122 _Nv = None # cached_toNvector L{Nvector})
124 def _update(self, updated, *attrs, **setters): # PYCHOK args
125 '''(INTERNAL) Zap cached attributes if updated.
126 '''
127 if updated: # reset caches
128 LatLonNvectorBase._update(self, updated, _Nv=self._Nv) # special case
129 LatLonSphericalBase._update(self, updated, *attrs, **setters)
131 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
132 '''Compute the (signed) distance from the start to the closest
133 point on the great circle line defined by a start and an
134 end point.
136 That is, if a perpendicular is drawn from this point to the
137 great circle line, the along-track distance is the distance
138 from the start point to the point where the perpendicular
139 crosses the line.
141 @arg start: Start point of great circle line (L{LatLon}).
142 @arg end: End point of great circle line (L{LatLon}) or
143 initial bearing from start point (compass
144 C{degrees360}).
145 @kwarg radius: Mean earth radius (C{meter}).
146 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
147 the B{C{start}} and B{C{end}} points (C{bool}).
149 @return: Distance along the great circle line (positive if
150 after the start toward the end point of the line
151 or negative if before the start point).
153 @raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.
155 @raise Valuerror: Some points coincide.
157 @example:
159 >>> p = LatLon(53.2611, -0.7972)
161 >>> s = LatLon(53.3206, -1.7297)
162 >>> e = LatLon(53.1887, 0.1334)
163 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58
164 '''
165 p = self.others(start=start)
166 n = self.toNvector()
168 gc, _, _ = self._gc3(p, end, _end_, wrap=wrap)
169 a = gc.cross(n).cross(gc) # along-track point gc × p × gc
170 return start.toNvector().angleTo(a, vSign=gc) * radius
172 @deprecated_method
173 def bearingTo(self, other, **unused): # PYCHOK no cover
174 '''DEPRECATED, use method L{initialBearingTo}.
175 '''
176 return self.initialBearingTo(other)
178 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
179 '''Compute the (signed) distance from this point to great circle
180 defined by a start and end point.
182 @arg start: Start point of great circle line (L{LatLon}).
183 @arg end: End point of great circle line (L{LatLon}) or
184 initial bearing from start point (compass
185 C{degrees360}).
186 @kwarg radius: Mean earth radius (C{meter}).
187 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
188 the B{C{start}} and B{C{end}} points (C{bool}).
190 @return: Distance to great circle (negative if to the
191 left or positive if to the right of the line).
193 @raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.
195 @raise Valuerror: Some points coincide.
197 @example:
199 >>> p = LatLon(53.2611, -0.7972)
201 >>> s = LatLon(53.3206, -1.7297)
202 >>> d = p.crossTrackDistanceTo(s, 96) # -305.7
204 >>> e = LatLon(53.1887, 0.1334)
205 >>> d = p.crossTrackDistanceTo(s, e) # -307.5
206 '''
207 p = self.others(start=start)
208 n = self.toNvector()
210 gc, _, _ = self._gc3(p, end, _end_, wrap=wrap)
211 return (gc.angleTo(n) - PI_2) * radius
213 def destination(self, distance, bearing, radius=R_M, height=None):
214 '''Locate the destination from this point after having
215 travelled the given distance on the given bearing.
217 @arg distance: Distance travelled (C{meter}, same units
218 as B{C{radius}}).
219 @arg bearing: Bearing from this point (compass C{degrees360}).
220 @kwarg radius: Mean earth radius (C{meter}).
221 @kwarg height: Optional height at destination, overriding the
222 default height (C{meter}, same units as B{C{radius}}).
224 @return: Destination point (L{LatLon}).
226 @raise Valuerror: Polar coincidence or invalid B{C{distance}},
227 B{C{bearing}}, B{C{radius}} or B{C{height}}.
229 @example:
231 >>> p = LatLon(51.4778, -0.0015)
232 >>> q = p.destination(7794, 300.7)
233 >>> q.toStr() # 51.513546°N, 000.098345°W
234 '''
235 b = Bearing_(bearing)
236 a = _angular(distance, radius, low=None)
237 sa, ca, sb, cb = sincos2_(a, b)
239 n = self.toNvector()
240 e = NorthPole.cross(n, raiser=_pole_).unit() # east vector at n
241 x = n.cross(e) # north vector at n
242 d = x.times(cb).plus(e.times(sb)) # direction vector @ n
243 n = n.times(ca).plus(d.times(sa))
244 return n.toLatLon(height=height, LatLon=self.classof) # Nvector(n.x, n.y, n.z).toLatLon(...)
246 def distanceTo(self, other, radius=R_M, wrap=False):
247 '''Compute the distance from this to an other point.
249 @arg other: The other point (L{LatLon}).
250 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
251 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
252 the B{C{other}} point (C{bool}).
254 @return: Distance between this and the B{C{other}} point
255 (C{meter}, same units as B{C{radius}} or C{radians}
256 if B{C{radius}} is C{None}).
258 @raise TypeError: Invalid B{C{other}} point.
260 @example:
262 >>> p = LatLon(52.205, 0.119)
263 >>> q = LatLon(48.857, 2.351);
264 >>> d = p.distanceTo(q) # 404.3 Km
265 '''
266 p = self.others(other)
267 if wrap:
268 p = _unrollon(self, p)
269 n = p.toNvector()
270 r = fabs(self.toNvector().angleTo(n, wrap=wrap))
271 return r if radius is None else (Radius(radius) * r)
273# @Property_RO
274# def Ecef(self):
275# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
276# '''
277# return _ALL_MODS.ecef.EcefKarney
279 def _gc3(self, start, end, namend, raiser=_point_, wrap=False):
280 '''(INTERNAL) Return great circle, start and end Nvectors.
281 '''
282 s = start.toNvector()
283 if isscalar(end): # bearing
284 gc = s.greatCircle(end)
285 e = None
286 else: # point
287 p = self.others(end, name=namend)
288 if wrap:
289 p = _unrollon(start, p, wrap=wrap)
290 e = p.toNvector()
291 gc = s.cross(e, raiser=raiser) # XXX .unit()?
292 return gc, s, e
294 def greatCircle(self, bearing):
295 '''Compute the vector normal to great circle obtained by
296 heading on the given bearing from this point.
298 Direction of vector is such that initial bearing vector
299 b = c × n, where n is an n-vector representing this point.
301 @arg bearing: Bearing from this point (compass C{degrees360}).
303 @return: N-vector representing the great circle (L{Nvector}).
304 '''
305 t = Bearing_(bearing)
306 a, b = self.philam
308 sa, ca, sb, cb, st, ct = sincos2_(a, b, t)
309 return Nvector(sb * ct - sa * cb * st,
310 -cb * ct - sa * sb * st,
311 ca * st, name=self.name) # XXX .unit()
313 def greatCircleTo(self, other, wrap=False):
314 '''Compute the vector normal to great circle obtained by
315 heading from this to an other point or on a given bearing.
317 Direction of vector is such that initial bearing vector
318 b = c × n, where n is an n-vector representing this point.
320 @arg other: The other point (L{LatLon}) or the bearing from
321 this point (compass C{degrees360}).
322 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
323 the B{C{other}} point (C{bool}).
325 @return: N-vector representing the great circle (L{Nvector}).
327 @raise TypeError: The B{C{other}} point is not L{LatLon}.
329 @raise Valuerror: Points coincide.
331 @example:
333 >>> p = LatLon(53.3206, -1.7297)
334 >>> gc = p.greatCircle(96.0)
335 >>> gc.toStr() # (-0.79408, 0.12856, 0.59406)
337 >>> q = LatLon(53.1887, 0.1334)
338 >>> g = p.greatCircleTo(q)
339 >>> g.toStr() # (-0.79408, 0.12859, 0.59406)
340 '''
341 gc, _, _ = self._gc3(self, other, _other_, wrap=wrap)
342 return gc.unit()
344 def initialBearingTo(self, other, wrap=False, **unused): # raiser=...
345 '''Compute the initial bearing (forward azimuth) from this
346 to an other point.
348 @arg other: The other point (L{LatLon}).
349 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
350 the B{C{other}} point (C{bool}).
352 @return: Initial bearing (compass C{degrees360}).
354 @raise Crosserror: This point coincides with the B{C{other}}
355 point or the C{NorthPole}, provided
356 L{pygeodesy.crosserrors} is C{True}.
358 @raise TypeError: The B{C{other}} point is not L{LatLon}.
360 @example:
362 >>> p1 = LatLon(52.205, 0.119)
363 >>> p2 = LatLon(48.857, 2.351)
364 >>> b = p1.initialBearingTo(p2) # 156.2
365 '''
366 n = self.toNvector()
367 p = self.others(other)
368 if wrap:
369 p = _unrollon(self, p, wrap=wrap)
370 p = p.toNvector()
371 # see <https://MathForum.org/library/drmath/view/55417.html>
372# gc1 = self.greatCircleTo(other)
373 gc1 = n.cross(p, raiser=_point_) # .unit()
374# gc2 = self.greatCircleTo(NorthPole)
375 gc2 = n.cross(NorthPole, raiser=_pole_) # .unit()
376 return degrees360(gc1.angleTo(gc2, vSign=n))
378 def intermediateChordTo(self, other, fraction, height=None, wrap=False):
379 '''Locate the point projected from the point at given fraction
380 on a straight line (chord) between this and an other point.
382 @arg other: The other point (L{LatLon}).
383 @arg fraction: Fraction between both points (float, between
384 0.0 for this and 1.0 for the other point).
385 @kwarg height: Optional height at the intermediate point,
386 overriding the fractional height (C{meter}).
387 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
388 the B{C{other}} point (C{bool}).
390 @return: Intermediate point (L{LatLon}).
392 @raise TypeError: The B{C{other}} point is not L{LatLon}.
394 @example:
396 >>> p = LatLon(52.205, 0.119)
397 >>> q = LatLon(48.857, 2.351)
398 >>> i = p.intermediateChordTo(q, 0.25) # 51.3723°N, 000.7072°E
399 '''
400 n = self.toNvector()
401 p = self.others(other)
402 if wrap:
403 p = _unrollon(self, p, wrap=wrap)
405 f = Scalar(fraction=fraction)
406 i = p.toNvector().times(f).plus(n.times(1 - f))
407# i = p.toNvector() * f + self.toNvector() * (1 - f))
409 h = self._havg(other, f=f, h=height)
410 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)
412 def intermediateTo(self, other, fraction, height=None, wrap=False):
413 '''Locate the point at a given fraction between this and an
414 other point.
416 @arg other: The other point (L{LatLon}).
417 @arg fraction: Fraction between both points (C{float}, between
418 0.0 for this and 1.0 for the other point).
419 @kwarg height: Optional height at the intermediate point,
420 overriding the fractional height (C{meter}).
421 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
422 the B{C{other}} point (C{bool}).
424 @return: Intermediate point (L{LatLon}).
426 @raise TypeError: The B{C{other}} point is not L{LatLon}.
428 @raise Valuerror: Points coincide or invalid B{C{height}}.
430 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
432 @example:
434 >>> p = LatLon(52.205, 0.119)
435 >>> q = LatLon(48.857, 2.351)
436 >>> i = p.intermediateTo(q, 0.25) # 51.3721°N, 000.7074°E
437 '''
438 n = self.toNvector()
439 p = self.others(other)
440 if wrap:
441 p = _unrollon(self, p, wrap=wrap)
442 p = p.toNvector()
443 f = Scalar(fraction=fraction)
445 x = n.cross(p, raiser=_point_)
446 d = x.unit().cross(n) # unit(n × p) × n
447 # angular distance α, tan(α) = |n × p| / n ⋅ p
448 s, c = sincos2(atan2(x.length, n.dot(p)) * f) # interpolated
449 i = n.times(c).plus(d.times(s)) # n * cosα + d * sinα
451 h = self._havg(other, f=f, h=height)
452 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)
454 def intersection(self, end1, start2, end2, height=None, wrap=False):
455 '''Locate the intersection point of two lines each defined
456 by two points or a start point and bearing from North.
458 @arg end1: End point of the first line (L{LatLon}) or the
459 initial bearing at this point (compass C{degrees360}).
460 @arg start2: Start point of the second line (L{LatLon}).
461 @arg end2: End point of the second line (L{LatLon}) or the
462 initial bearing at the second point (compass
463 C{degrees}).
464 @kwarg height: Optional height at the intersection point,
465 overriding the mean height (C{meter}).
466 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll all
467 start and end points (C{bool}).
469 @return: The intersection point (L{LatLon}).
471 @raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}}
472 point is not L{LatLon}.
474 @raise ValueError: Intersection is ambiguous or infinite or
475 the lines are parallel, coincident or null.
477 @see: Function L{sphericalNvector.intersection} and method
478 L{intersection2}.
480 @example:
482 >>> s = LatLon(51.8853, 0.2545)
483 >>> e = LatLon(49.0034, 2.5735)
484 >>> i = s.intersection(108.55, e, 32.44) # 50.9076°N, 004.5086°E
485 '''
486 return intersection(self, end1, start2, end2, height=height,
487 wrap=wrap, LatLon=self.classof)
489 def intersection2(self, end1, start2, end2, height=None, wrap=False):
490 '''Locate the intersections of two (great circle) lines each defined
491 by two points or by a start point and an (initial) bearing.
493 @arg end1: End point of the first line (L{LatLon}) or the
494 initial bearing at this point (compass C{degrees360}).
495 @arg start2: Start point of the second line (L{LatLon}).
496 @arg end2: End point of the second line (L{LatLon}) or the
497 initial bearing at the second start point (compass
498 C{degrees360}).
499 @kwarg height: Optional height at the intersection and antipodal
500 point, overriding the mean height (C{meter}).
501 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
502 B{C{start2}} and both B{C{end*}} points (C{bool}).
504 @return: 2-Tuple C{(intersection, antipode)}, each a B{C{LatLon}}.
506 @raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}}
507 point is not L{LatLon}.
509 @raise ValueError: Intersection is ambiguous or infinite or
510 the lines are parallel, coincident or null.
512 @see: Function L{sphericalNvector.intersection2} and method
513 L{intersection}.
514 '''
515 return intersection2(self, end1, start2, end2, height=height,
516 wrap=wrap, LatLon=self.classof)
518 def isenclosedBy(self, points, wrap=False):
519 '''Check whether a (convex) polygon or composite encloses this point.
521 @arg points: The polygon points or composite (L{LatLon}[],
522 L{BooleanFHP} or L{BooleanGH}).
523 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
524 B{C{points}} (C{bool}).
526 @return: C{True} if this point is inside the polygon or composite,
527 C{False} otherwise.
529 @raise PointsError: Insufficient number of B{C{points}}.
531 @raise TypeError: Some B{C{points}} are not L{LatLon}.
533 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
534 and L{pygeodesy.ispolar} especially if the B{C{points}} may
535 enclose a pole or wrap around the earth I{longitudinally}.
536 '''
537 if _MODS.booleans.isBoolean(points):
538 return points._encloses(self.lat, self.lon, wrap=wrap)
540 # sum subtended angles of each edge (using n0, the
541 # normal vector to this point for sign of α)
542 def _subtangles(ps, w):
543 Ps = self.PointsIter(ps, loop=1, wrap=w)
544 n0 = self.toNvector()
545 _m0 = n0.minus
546 p1 = Ps[0]
547 vs1 = _m0(p1.toNvector())
548 for p2 in Ps.iterate(closed=True):
549 if w and not Ps.looped:
550 p2 = _unrollon(p1, p2)
551 p1 = p2
552 vs2 = _m0(p2.toNvector())
553 yield vs1.angleTo(vs2, vSign=n0) # PYCHOK false
554 vs1 = vs2
556 # Note, this method uses angle summation test: on a plane,
557 # angles for an enclosed point will sum to 360°, angles for
558 # an exterior point will sum to 0°. On a sphere, enclosed
559 # point angles will sum to less than 360° (due to spherical
560 # excess), exterior point angles will be small but non-zero.
561 s = fsum(_subtangles(points, wrap), floats=True) # normal vector
562 # XXX are winding number optimisations equally applicable to
563 # spherical surface?
564 return fabs(s) > PI
566 @deprecated_method
567 def isEnclosedBy(self, points): # PYCHOK no cover
568 '''DEPRECATED, use method C{isenclosedBy}.'''
569 return self.isenclosedBy(points)
571 def iswithin(self, point1, point2, wrap=False):
572 '''Check whether this point is between two other points.
574 If this point is not on the great circle arc defined by
575 both points, return whether it is within the area bound
576 by perpendiculars to the great circle at each point (in
577 the same hemispere).
579 @arg point1: Start point of the arc (L{LatLon}).
580 @arg point2: End point of the arc (L{LatLon}).
581 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
582 B{C{point1}} and B{C{point2}} (C{bool}).
584 @return: C{True} if this point is within the (great circle)
585 arc, C{False} otherwise.
587 @raise TypeError: If B{C{point1}} or B{C{point2}} is not
588 L{LatLon}.
589 '''
590 p1 = self.others(point1=point1)
591 p2 = self.others(point2=point2)
592 if wrap:
593 p1 = _Wrap.point(p1)
594 p2 = _unrollon(p1, p2, wrap=wrap)
595 n, n1, n2 = (_.toNvector() for _ in (self, p1, p2))
597 # corner case, null arc
598 if n1.isequalTo(n2):
599 return n.isequalTo(n1) or n.isequalTo(n2) # PYCHOK returns
601 if n.dot(n1) < 0 or n.dot(n2) < 0: # different hemisphere
602 return False # PYCHOK returns
604 # get vectors representing d0=p0->p1 and d2=p2->p1 and the
605 # dot product d0⋅d2 tells us if p0 is on the p2 side of p1 or
606 # on the other side (similarly for d0=p0->p2 and d1=p1->p2
607 # and dot product d0⋅d1 and p0 on the p1 side of p2 or not)
608 return n.minus(n1).dot(n2.minus(n1)) >= 0 and \
609 n.minus(n2).dot(n1.minus(n2)) >= 0
611 @deprecated_method
612 def isWithin(self, point1, point2): # PYCHOK no cover
613 '''DEPRECATED, use method C{iswithin}.'''
614 return self.iswithin(point1, point2)
616 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
617 '''Find the midpoint between this and an other point.
619 @arg other: The other point (L{LatLon}).
620 @kwarg height: Optional height at the midpoint, overriding
621 the mean height (C{meter}).
622 @kwarg fraction: Midpoint location from this point (C{scalar}),
623 may be negative or greater than 1.0.
624 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
625 B{C{other}} point (C{bool}).
627 @return: Midpoint (L{LatLon}).
629 @raise TypeError: The B{C{other}} point is not L{LatLon}.
631 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
633 @example:
635 >>> p1 = LatLon(52.205, 0.119)
636 >>> p2 = LatLon(48.857, 2.351)
637 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
638 '''
639 if fraction is _0_5:
640 p = self.others(other)
641 if wrap:
642 p = _unrollon(self, p, wrap=wrap)
643 m = self.toNvector().plus(p.toNvector())
644 h = self._havg(other, f=fraction, h=height)
645 r = m.toLatLon(height=h, LatLon=self.classof)
646 else:
647 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
648 return r
650 def nearestOn(self, point1, point2, height=None, within=True, wrap=False):
651 '''Locate the point on the great circle arc between two points
652 closest to this point.
654 @arg point1: Start point of the arc (L{LatLon}).
655 @arg point2: End point of the arc (L{LatLon}).
656 @kwarg height: Optional height, overriding the mean height for
657 the point within the arc (C{meter}), or C{None}
658 to interpolate the height.
659 @kwarg within: If C{True}, return the closest point between both
660 given points, otherwise the closest point
661 elsewhere on the great circle arc (C{bool}).
662 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
663 B{C{point1}} and B{C{point2}} (C{bool}).
665 @return: Closest point on the arc (L{LatLon}).
667 @raise NotImplementedError: Keyword argument C{B{wrap}=True}
668 not supported.
670 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
672 @example:
674 >>> s1 = LatLon(51.0, 1.0)
675 >>> s2 = LatLon(51.0, 2.0)
677 >>> s = LatLon(51.0, 1.9)
678 >>> p = s.nearestOn(s1, s2) # 51.0004°N, 001.9000°E
680 >>> d = p.distanceTo(s) # 42.71 m
682 >>> s = LatLon(51.0, 2.1)
683 >>> p = s.nearestOn(s1, s2) # 51.0000°N, 002.0000°E
684 '''
685 p1 = self.others(point1=point1)
686 p2 = self.others(point2=point2)
687 if wrap:
688 p1 = _Wrap.point(p1)
689 p2 = _unrollon(p1, p2, wrap=wrap)
690 p0 = self
692 if p0.iswithin(p1, p2) and not p1.isequalTo(p2, EPS):
693 # closer to arc than to its endpoints,
694 # find the closest point on the arc
695 gc1 = p1.toNvector().cross(p2.toNvector())
696 gc2 = p0.toNvector().cross(gc1)
697 n = gc1.cross(gc2)
699 elif within: # for backward compatibility, XXX unwrapped
700 return point1 if (self.distanceTo(point1) <
701 self.distanceTo(point2)) else point2
703 else: # handle beyond arc extent by .vector3d.nearestOn
704 n1 = p1.toNvector()
705 n2 = p2.toNvector()
706 n = p0.toNvector().nearestOn(n1, n2, within=False)
707 if n is n1:
708 return p1 # is point1
709 elif n is n2:
710 return p2 # is point2 if not wrap
712 p = n.toLatLon(height=height or 0, LatLon=self.classof)
713 if height in (None, False): # interpolate height within extent
714 d = p1.distanceTo(p2)
715 f = (p1.distanceTo(p) / d) if d > EPS0 else _0_5
716 p.height = p1._havg(p2, f=max(_0_0, min(f, _1_0)))
717 return p
719 # @deprecated_method
720 def nearestOn2(self, points, **closed_radius_height): # PYCHOK no cover
721 '''DEPRECATED, use method L{sphericalNvector.LatLon.nearestOn3}.
723 @return: ... 2-Tuple C{(closest, distance)} of the C{closest}
724 point (L{LatLon}) on the polygon and the C{distance}
725 to that point from this point ...
726 '''
727 r = self.nearestOn3(points, **closed_radius_height)
728 return r.closest, r.distance
730 def nearestOn3(self, points, closed=False, radius=R_M, height=None, wrap=False):
731 '''Locate the point on a path or polygon (with great circle
732 arcs joining consecutive points) closest to this point.
734 The closest point is either on within the extent of any great
735 circle arc or the nearest of the arc's end points.
737 @arg points: The path or polygon points (L{LatLon}[]).
738 @kwarg closed: Optionally, close the polygon (C{bool}).
739 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
740 @kwarg height: Optional height, overriding the mean height
741 for a point within the arc (C{meter}).
742 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
743 the B{C{points}} (C{bool}).
745 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of
746 the C{closest} point (L{LatLon}), the C{distance}
747 between this and the C{closest} point in C{meter},
748 same units as B{C{radius}} or in C{radians} if
749 B{C{radius}} is C{None} and the C{angle} from this
750 to the C{closest} point in compass C{degrees360}.
752 @raise TypeError: Some B{C{points}} are not C{LatLon}.
754 @raise ValueError: No B{C{points}}.
755 '''
756 Ps = self.PointsIter(points, loop=1, wrap=wrap)
757 _r = self.distanceTo
758 _n = self.nearestOn
760 c = p1 = Ps[0]
761 r = _r(c, radius=None) # radians
762 for p2 in Ps.iterate(closed=closed):
763 if wrap and not Ps.looped:
764 p2 = _unrollon(p1, p2)
765 p = _n(p1, p2, height=height)
766 d = _r(p, radius=None) # radians
767 if d < r:
768 c, r = p, d
769 p1 = p2
770 d = r if radius is None else (Radius(radius) * r)
771 return NearestOn3Tuple(c, d, degrees360(r))
773 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian, datum=None
774 '''Convert this point to C{Nvector}-based cartesian (ECEF) coordinates.
776 @kwarg Cartesian_and_kwds: Optional L{Cartesian} and L{Cartesian} keyword
777 arguments, like C{datum}. Use C{B{Cartesian}=...}
778 to override this L{Cartesian} class or specify
779 C{B{Cartesian}=None}.
781 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} is
782 set to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
783 C, M, datum)} with C{C} and C{M} if available.
785 @raise TypeError: Invalid B{C{Cartesian_and_kwds}} argument.
786 '''
787 kwds = _xkwds(Cartesian_and_kwds, Cartesian=Cartesian, datum=self.datum)
788 return LatLonSphericalBase.toCartesian(self, **kwds)
790 def toNvector(self, **Nvector_and_kwds): # PYCHOK signature
791 '''Convert this point to L{Nvector} components, I{including height}.
793 @kwarg Nvector_and_kwds: Optional L{Nvector} and L{Nvector} keyword
794 arguments. Specify C{B{Nvector}=...} to
795 override this L{Nvector} class or use
796 C{B{Nvector}=None}.
798 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}} is
799 set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}.
801 @raise TypeError: Invalid B{C{Nvector_and_kwds}} argument.
803 @example:
805 >>> p = LatLon(45, 45)
806 >>> n = p.toNvector()
807 >>> n.toStr() # [0.50, 0.50, 0.70710]
808 '''
809 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector)
810 return LatLonNvectorBase.toNvector(self, **kwds)
813class Nvector(NvectorBase):
814 '''An n-vector is a position representation using a (unit) vector
815 normal to the earth's surface. Unlike lat-/longitude points,
816 n-vectors have no singularities or discontinuities.
818 For many applications, n-vectors are more convenient to work
819 with than other position representations like lat-/longitude,
820 earth-centred earth-fixed (ECEF) vectors, UTM coordinates, etc.
822 On a spherical model earth, an n-vector is equivalent to an
823 earth-centred earth-fixed (ECEF) vector.
825 Note commonality with L{ellipsoidalNvector.Nvector}.
826 '''
827 _datum = Datums.Sphere # default datum (L{Datum})
829 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian
830 '''Convert this n-vector to C{Nvector}-based cartesian
831 (ECEF) coordinates.
833 @kwarg Cartesian_and_kwds: Optional L{Cartesian} and L{Cartesian} keyword
834 arguments, like C{h}. Use C{B{Cartesian}=...}
835 to override this L{Cartesian} class or specify
836 C{B{Cartesian}=None}.
838 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} is
839 set to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
840 C, M, datum)} with C{C} and C{M} if available.
842 @raise TypeError: Invalid B{C{Cartesian_and_kwds}} argument.
843 '''
844 kwds = _xkwds(Cartesian_and_kwds, h=self.h, Cartesian=Cartesian)
845 return NvectorBase.toCartesian(self, **kwds) # class or .classof
847 def toLatLon(self, **LatLon_and_kwds): # PYCHOK height=None, LatLon=LatLon
848 '''Convert this n-vector to an C{Nvector}-based geodetic point.
850 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
851 arguments, like C{height}. Use C{B{LatLon}=...}
852 to override this L{LatLon} class or specify
853 C{B{LatLon}=None}.
855 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set
856 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
857 C, M, datum)} with C{C} and C{M} if available.
859 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
861 @raise ValueError: Invalid B{C{height}}.
862 '''
863 kwds = _xkwds(LatLon_and_kwds, height=self.h, LatLon=LatLon)
864 return NvectorBase.toLatLon(self, **kwds) # class or .classof
866 def greatCircle(self, bearing):
867 '''Compute the n-vector normal to great circle obtained by
868 heading on given compass bearing from this point as its
869 n-vector.
871 Direction of vector is such that initial bearing vector
872 b = c × p.
874 @arg bearing: Initial compass bearing (C{degrees}).
876 @return: N-vector representing great circle (L{Nvector}).
878 @raise Valuerror: Polar coincidence.
880 @example:
882 >>> n = LatLon(53.3206, -1.7297).toNvector()
883 >>> gc = n.greatCircle(96.0) # [-0.794, 0.129, 0.594]
884 '''
885 s, c = sincos2d(Bearing(bearing))
887 e = NorthPole.cross(self, raiser=_pole_) # easting
888 n = self.cross(e, raiser=_point_) # northing
890 e = e.times(c / e.length)
891 n = n.times(s / n.length)
892 return n.minus(e)
895_Nvll = LatLon(_0_0, _0_0, name=_Nv00_) # reference instance (L{LatLon})
898def areaOf(points, radius=R_M, wrap=False):
899 '''Calculate the area of a (spherical) polygon or composite (with
900 great circle arcs joining consecutive points).
902 @arg points: The polygon points or clips (C{LatLon}[],
903 L{BooleanFHP} or L{BooleanGH}).
904 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
905 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
906 B{C{points}} (C{bool}).
908 @return: Polygon area (C{meter} I{squared} , same units as
909 B{C{radius}}, or C{radians} if B{C{radius}} is C{None}).
911 @raise PointsError: Insufficient number of B{C{points}}.
913 @raise TypeError: Some B{C{points}} are not L{LatLon}.
915 @see: Functions L{pygeodesy.areaOf}, L{sphericalTrigonometry.areaOf}
916 and L{ellipsoidalKarney.areaOf}.
918 @example:
920 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
921 >>> areaOf(b) # 8666058750.718977
922 '''
923 def _interangles(ps, w): # like .karney._polygon
924 Ps = _Nvll.PointsIter(ps, loop=2, wrap=w)
925 # use vector to 1st point as plane normal for sign of α
926 n0 = Ps[0].toNvector()
928 v2 = Ps[0]._N_vector # XXX v2 == no?
929 p1 = Ps[1]
930 v1 = p1._N_vector
931 gc = v2.cross(v1)
932 for p2 in Ps.iterate(closed=True):
933 if w and not Ps.looped:
934 p2 = _unrollon(p1, p2)
935 p1 = p2
936 v2 = p2._N_vector
937 gc1 = v1.cross(v2)
938 v1 = v2
939 yield gc.angleTo(gc1, vSign=n0)
940 gc = gc1
942 if _MODS.booleans.isBoolean(points):
943 r = points._sum2(LatLon, areaOf, radius=None, wrap=wrap)
944 else:
945 # sum interior angles: depending on whether polygon is cw or ccw,
946 # angle between edges is π−α or π+α, where α is angle between
947 # great-circle vectors; so sum α, then take n·π − |Σα| (cannot
948 # use Σ(π−|α|) as concave polygons would fail)
949 s = fsum(_interangles(points, wrap), floats=True)
950 # using Girard’s theorem: A = [Σθᵢ − (n−2)·π]·R²
951 # (PI2 - abs(s) == (n*PI - abs(s)) - (n-2)*PI)
952 r = fabs(PI2 - fabs(s))
953 return r if radius is None else (r * Radius(radius)**2)
956def intersecant2(center, circle, point, bearing, radius=R_M, exact=False,
957 height=None, wrap=False): # was=True
958 '''Compute the intersections of a circle and a line.
960 @arg center: Center of the circle (L{LatLon}).
961 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
962 or a point on the circle (L{LatLon}).
963 @arg point: A point in- or outside the circle (L{LatLon}).
964 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or
965 a second point on the line (L{LatLon}).
966 @kwarg radius: Mean earth radius (C{meter}, conventionally).
967 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
968 destination and distance, if C{False} use the basic
969 rhumb methods (C{bool}) or if C{None} use the I{great
970 circle} methods.
971 @kwarg height: Optional height for the intersection points (C{meter},
972 conventionally) or C{None}.
973 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}}
974 and the B{C{circle}} and B{C{bearing}} points (C{bool}).
976 @return: 2-Tuple of the intersection points (representing a chord), each
977 an instance of this class. For a tangent line, both points are
978 the same instance, the B{C{point}} or wrapped or I{normalized}.
980 @raise IntersectionError: The circle and line do not intersect.
982 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
983 B{C{circle}} or B{C{bearing}} invalid.
985 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
986 B{C{exact}} or B{C{height}}.
987 '''
988 c = _Nvll.others(center=center)
989 p = _Nvll.others(point=point)
990 try:
991 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact,
992 height=height, wrap=wrap)
993 except (TypeError, ValueError) as x:
994 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing,
995 exact=exact, wrap=wrap)
998def intersection(start1, end1, start2, end2, height=None, wrap=False,
999 LatLon=LatLon, **LatLon_kwds):
1000 '''Locate the intersections of two (great circle) lines each defined
1001 by two points or by a start point and an (initial) bearing.
1003 @arg start1: Start point of the first line (L{LatLon}).
1004 @arg end1: End point of the first line (L{LatLon}) or the initial
1005 bearing at the first start point (compass C{degrees360}).
1006 @arg start2: Start point of the second line (L{LatLon}).
1007 @arg end2: End point of the second line (L{LatLon}) or the initial
1008 bearing at the second start point (compass C{degrees360}).
1009 @kwarg height: Optional height at the intersection point,
1010 overriding the mean height (C{meter}).
1011 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{start2}}
1012 and both B{C{end*}} points (C{bool}).
1013 @kwarg LatLon: Optional class to return the intersection point
1014 (L{LatLon}).
1015 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1016 arguments, ignored if C{B{LatLon} is None}.
1018 @return: The intersection point (B{C{LatLon}}) or if C{B{LatLon}
1019 is None}, a cartesian L{Ecef9Tuple}C{(x, y, z, lat, lon,
1020 height, C, M, datum)} with C{C} and C{M} if available.
1022 @raise TypeError: If B{C{start*}} or B{C{end*}} is not L{LatLon}.
1024 @raise ValueError: Intersection is ambiguous or infinite or
1025 the lines are parallel, coincident or null.
1027 @see: Function L{sphericalNvector.intersection2}.
1029 @example:
1031 >>> p = LatLon(51.8853, 0.2545)
1032 >>> q = LatLon(49.0034, 2.5735)
1033 >>> i = intersection(p, 108.55, q, 32.44) # 50.9076°N, 004.5086°E
1034 '''
1035 i, _, h = _intersect3(start1, end1, start2, end2, height, wrap)
1036 kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)
1037 return i.toLatLon(**kwds)
1040def intersection2(start1, end1, start2, end2, height=None, wrap=False,
1041 LatLon=LatLon, **LatLon_kwds):
1042 '''Locate the intersections of two (great circle) lines each defined
1043 by two points or by a start point and an (initial) bearing.
1045 @arg start1: Start point of the first line (L{LatLon}).
1046 @arg end1: End point of the first line (L{LatLon}) or the
1047 initial bearing at the first start point
1048 (compass C{degrees360}).
1049 @arg start2: Start point of the second line (L{LatLon}).
1050 @arg end2: End point of the second line (L{LatLon}) or the
1051 initial bearing at the second start point
1052 (compass C{degrees360}).
1053 @kwarg height: Optional height at the intersection and antipodal
1054 point, overriding the mean height (C{meter}).
1055 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{start2}}
1056 and both B{C{end*}} points (C{bool}).
1057 @kwarg LatLon: Optional class to return the intersection and
1058 antipodal points (L{LatLon}).
1059 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1060 arguments, ignored if C{B{LatLon} is None}.
1062 @return: 2-Tuple C{(intersection, antipode)}, each a (B{C{LatLon}})
1063 or if C{B{LatLon} is None}, a cartesian L{Ecef9Tuple}C{(x,
1064 y, z, lat, lon, height, C, M, datum)} with C{C} and C{M}
1065 if available.
1067 @raise TypeError: If B{C{start*}} or B{C{end*}} is not L{LatLon}.
1069 @raise ValueError: Intersection is ambiguous or infinite or
1070 the lines are parallel, coincident or null.
1072 @see: Function L{sphericalNvector.intersection}.
1073 '''
1074 i, a, h = _intersect3(start1, end1, start2, end2, height, wrap)
1075 kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)
1076 return i.toLatLon(**kwds), a.toLatLon(**kwds)
1079def _intersect3(start1, end1, start2, end2, height, wrap):
1080 '''(INTERNAL) Return the intersection and antipodal points for
1081 functions C{intersection} and C{intersection2}.
1082 '''
1083 p1 = _Nvll.others(start1=start1)
1084 p2 = _Nvll.others(start2=start2)
1085 if wrap:
1086 p2 = _unrollon(p1, p2, wrap=wrap)
1087 # If gc1 and gc2 are great circles through start and end points
1088 # (or defined by start point and bearing), then the candidate
1089 # intersections are simply gc1 × gc2 and gc2 × gc1. Most of the
1090 # work is deciding the correct intersection point to select! If
1091 # bearing is given, that determines the intersection, but if both
1092 # lines are defined by start/end points, take closer intersection.
1093 gc1, s1, e1 = _Nvll._gc3(p1, end1, 'end1', wrap=wrap)
1094 gc2, s2, e2 = _Nvll._gc3(p2, end2, 'end2', wrap=wrap)
1096 hs = start1.height, start2.height
1097 # there are two (antipodal) candidate intersection
1098 # points ... we have to choose the one to return
1099 i1 = gc1.cross(gc2, raiser=_lines_)
1100 i2 = gc2.cross(gc1, raiser=_lines_)
1102 # selection of intersection point depends on how
1103 # lines are defined (by bearings or endpoints)
1104 if e1 and e2: # endpoint+endpoint
1105 d = sumOf((s1, s2, e1, e2)).dot(i1)
1106 hs += end1.height, end2.height
1107 elif e1 and not e2: # endpoint+bearing
1108 # gc2 x v2 . i1 +ve means v2 bearing points to i1
1109 d = gc2.cross(s2).dot(i1)
1110 hs += end1.height,
1111 elif e2 and not e1: # bearing+endpoint
1112 # gc1 x v1 . i1 +ve means v1 bearing points to i1
1113 d = gc1.cross(s1).dot(i1)
1114 hs += end2.height,
1115 else: # bearing+bearing
1116 # if gc x v . i1 is +ve, initial bearing is
1117 # towards i1, otherwise towards antipodal i2
1118 d1 = gc1.cross(s1).dot(i1) # +ve means p1 bearing points to i1
1119 d2 = gc2.cross(s2).dot(i1) # +ve means p2 bearing points to i1
1120 if d1 > 0 and d2 > 0:
1121 d = 1 # both point to i1
1122 elif d1 < 0 and d2 < 0:
1123 d = -1 # both point to i2
1124 else: # d1, d2 opposite signs
1125 # intersection is at further-away intersection point,
1126 # take opposite intersection from mid- point of v1
1127 # and v2 [is this always true?] XXX changed to always
1128 # get intersection p1 bearing points to, aka being
1129 # located "after" p1 along the bearing at p1, like
1130 # function .sphericalTrigonometry._intersect and
1131 # .ellipsoidalBaseDI._intersect3
1132 d = d1 # neg(s1.plus(s2).dot(i1))
1134 h = fmean(hs) if height is None else height
1135 return (i1, i2, h) if d > 0 else (i2, i1, h)
1138def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):
1139 '''Compute the I{geographic} mean of the supplied points.
1141 @arg points: Array of points to be averaged (L{LatLon}[]).
1142 @kwarg height: Optional height, overriding the mean height
1143 (C{meter}).
1144 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{points}} (C{bool}).
1145 @kwarg LatLon: Optional class to return the mean point (L{LatLon}).
1146 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1147 arguments, ignored if C{B{LatLon} is None}.
1149 @return: Point at geographic mean and mean height (B{C{LatLon}}).
1151 @raise PointsError: Insufficient number of B{C{points}} or
1152 some B{C{points}} are not C{LatLon}.
1153 '''
1154 def _N_vs(ps, w):
1155 Ps = _Nvll.PointsIter(ps, wrap=w)
1156 for p in Ps.iterate(closed=False):
1157 yield p._N_vector
1159 try:
1160 # geographic mean
1161 n = _nsumOf(_N_vs(points, wrap), height, Nvector, {})
1162 except (TypeError, ValueError) as x:
1163 raise PointsError(points=points, wrap=wrap, LatLon=LatLon, cause=x)
1164 return n.toLatLon(**_xkwds(LatLon_kwds, LatLon=LatLon, height=n.h,
1165 name=meanOf.__name__))
1168@deprecated_function
1169def nearestOn2(point, points, **closed_radius_height): # PYCHOK no cover
1170 '''DEPRECATED, use method L{sphericalNvector.nearestOn3}.
1172 @return: ... 2-Tuple C{(closest, distance)} of the C{closest}
1173 point (L{LatLon}) on the polygon and the C{distance}
1174 between the C{closest} and the given B{C{point}} ...
1175 '''
1176 r = nearestOn3(point, points, **closed_radius_height)
1177 return r.closest, r.distance
1180def nearestOn3(point, points, closed=False, radius=R_M, height=None, wrap=False):
1181 '''Locate the point on a polygon (with great circle arcs joining
1182 consecutive points) closest to an other point.
1184 If the given point is between the end points of a great circle
1185 arc, the closest point is on that arc. Otherwise, the closest
1186 point is the nearest of the arc's end points.
1188 @arg point: The other, reference point (L{LatLon}).
1189 @arg points: The polygon points (L{LatLon}[]).
1190 @kwarg closed: Optionally, close the polygon (C{bool}).
1191 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1192 @kwarg height: Optional height, overriding the mean height for
1193 a point within the (great circle) arc (C{meter}).
1194 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1195 B{C{points}} (C{bool}).
1197 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of
1198 the C{closest} point (L{LatLon}) on the polygon, the
1199 C{distance} and the C{angle} between the C{closest}
1200 and the given B{C{point}}. The C{distance} is in
1201 C{meter}, same units as B{C{radius}} or in C{radians}
1202 if B{C{radius}} is C{None}, the C{angle} is in compass
1203 C{degrees360}.
1205 @raise PointsError: Insufficient number of B{C{points}}.
1207 @raise TypeError: Some B{C{points}} or B{C{point}} not C{LatLon}.
1208 '''
1209 _xinstanceof(LatLon, point=point)
1211 return point.nearestOn3(points, closed=closed, radius=radius,
1212 height=height, wrap=wrap)
1215def perimeterOf(points, closed=False, radius=R_M, wrap=False):
1216 '''Compute the perimeter of a (spherical) polygon or composite (with
1217 great circle arcs joining consecutive points).
1219 @arg points: The polygon points (L{LatLon}[]).
1220 @kwarg closed: Optionally, close the polygon (C{bool}).
1221 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1222 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1223 B{C{points}} (C{bool}).
1225 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1226 or C{radians} if B{C{radius}} is C{None}).
1228 @raise PointsError: Insufficient number of B{C{points}}.
1230 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1232 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1233 C{B{points}} a composite.
1235 @see: Functions L{pygeodesy.perimeterOf}, L{ellipsoidalKarney.perimeterOf}
1236 and L{sphericalTrigonometry.perimeterOf}.
1237 '''
1238 def _rads(ps, c, w): # angular edge lengths in radians
1239 Ps = _Nvll.PointsIter(ps, loop=1, wrap=w)
1240 p1 = Ps[0]
1241 v1 = p1._N_vector
1242 for p2 in Ps.iterate(closed=c):
1243 if w and not (c and Ps.looped):
1244 p2 = _unrollon(p1, p2)
1245 p1 = p2
1246 v2 = p2._N_vector
1247 yield v1.angleTo(v2)
1248 v1 = v2
1250 if _MODS.booleans.isBoolean(points):
1251 if not closed:
1252 notImplemented(None, closed=closed, points=_composite_)
1253 r = points._sum2(LatLon, perimeterOf, closed=True, radius=None, wrap=wrap)
1254 else:
1255 r = fsum(_rads(points, closed, wrap), floats=True)
1256 return r if radius is None else (Radius(radius) * r)
1259def sumOf(nvectors, Vector=Nvector, h=None, **Vector_kwds):
1260 '''Return the I{vectorial} sum of two or more n-vectors.
1262 @arg nvectors: Vectors to be added (L{Nvector}[]).
1263 @kwarg Vector: Optional class for the vectorial sum (L{Nvector}).
1264 @kwarg h: Optional height, overriding the mean height (C{meter}).
1265 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword arguments.
1267 @return: Vectorial sum (B{C{Vector}}).
1269 @raise VectorError: No B{C{nvectors}}.
1270 '''
1271 try:
1272 return _nsumOf(nvectors, h, Vector, Vector_kwds)
1273 except (TypeError, ValueError) as x:
1274 raise VectorError(nvectors=nvectors, Vector=Vector, cause=x)
1277def triangulate(point1, bearing1, point2, bearing2,
1278 height=None, wrap=False,
1279 LatLon=LatLon, **LatLon_kwds):
1280 '''Locate a point given two known points and the initial bearings
1281 from those points.
1283 @arg point1: First reference point (L{LatLon}).
1284 @arg bearing1: Bearing at the first point (compass C{degrees360}).
1285 @arg point2: Second reference point (L{LatLon}).
1286 @arg bearing2: Bearing at the second point (compass C{degrees360}).
1287 @kwarg height: Optional height at the triangulated point, overriding
1288 the mean height (C{meter}).
1289 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{point2}}
1290 (C{bool}).
1291 @kwarg LatLon: Optional class to return the triangulated point (L{LatLon}).
1292 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1293 arguments, ignored if C{B{LatLon} is None}.
1295 @return: Triangulated point (B{C{LatLon}}).
1297 @raise TypeError: If B{C{point1}} or B{C{point2}} is not L{LatLon}.
1299 @raise Valuerror: Points coincide.
1301 @example:
1303 >>> p = LatLon("47°18.228'N","002°34.326'W") # Basse Castouillet
1304 >>> q = LatLon("47°18.664'N","002°31.717'W") # Basse Hergo
1305 >>> t = triangulate(p, 7, q, 295) # 47.323667°N, 002.568501°W'
1306 '''
1307 return _triangulate(_Nvll.others(point1=point1), bearing1,
1308 _Nvll.others(point2=point2), bearing2,
1309 height=height, wrap=wrap,
1310 LatLon=LatLon, **LatLon_kwds)
1313def trilaterate(point1, distance1, point2, distance2, point3, distance3, # PYCHOK args
1314 radius=R_M, height=None, useZ=False, wrap=False,
1315 LatLon=LatLon, **LatLon_kwds):
1316 '''Locate a point at given distances from three other points.
1318 @arg point1: First point (L{LatLon}).
1319 @arg distance1: Distance to the first point (C{meter}, same units
1320 as B{C{radius}}).
1321 @arg point2: Second point (L{LatLon}).
1322 @arg distance2: Distance to the second point (C{meter}, same units
1323 as B{C{radius}}).
1324 @arg point3: Third point (L{LatLon}).
1325 @arg distance3: Distance to the third point (C{meter}, same units
1326 as B{C{radius}}).
1327 @kwarg radius: Mean earth radius (C{meter}).
1328 @kwarg height: Optional height at the trilaterated point, overriding
1329 the IDW height (C{meter}, same units as B{C{radius}}).
1330 @kwarg useZ: Include Z component iff non-NaN, non-zero (C{bool}).
1331 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{point2}}
1332 and B{C{point3}} (C{bool}).
1333 @kwarg LatLon: Optional class to return the trilaterated point (L{LatLon}).
1334 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments,
1335 ignored if C{B{LatLon} is None}.
1337 @return: Trilaterated point (B{C{LatLon}}).
1339 @raise IntersectionError: No intersection, trilateration failed.
1341 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or B{C{point3}}.
1343 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}},
1344 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.
1346 @see: U{Trilateration<https://WikiPedia.org/wiki/Trilateration>}.
1347 '''
1348 return _trilaterate(_Nvll.others(point1=point1), distance1,
1349 _Nvll.others(point2=point2), distance2,
1350 _Nvll.others(point3=point3), distance3,
1351 radius=radius, height=height, useZ=useZ,
1352 wrap=wrap, LatLon=LatLon, **LatLon_kwds)
1355__all__ += _ALL_OTHER(Cartesian, LatLon, Nvector, # classes
1356 areaOf, # functions
1357 intersecant2, intersection, intersection2, ispolar,
1358 meanOf,
1359 nearestOn2, nearestOn3,
1360 perimeterOf,
1361 sumOf,
1362 triangulate, trilaterate)
1364# **) MIT License
1365#
1366# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1367#
1368# Permission is hereby granted, free of charge, to any person obtaining a
1369# copy of this software and associated documentation files (the "Software"),
1370# to deal in the Software without restriction, including without limitation
1371# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1372# and/or sell copies of the Software, and to permit persons to whom the
1373# Software is furnished to do so, subject to the following conditions:
1374#
1375# The above copyright notice and this permission notice shall be included
1376# in all copies or substantial portions of the Software.
1377#
1378# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1379# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1380# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1381# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1382# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1383# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1384# OTHER DEALINGS IN THE SOFTWARE.