Coverage for pygeodesy/sphericalNvector.py: 92%
320 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-12 12:31 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-12 12:31 -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, _r2m
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.30'
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}) or C{None}.
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 (C{radians}
150 if C{B{radius} is None} else C{meter}, same units
151 as B{C{radius}}), positive if "after" the start
152 toward the end point of the line or negative if
153 "before" the start point.
155 @raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.
157 @raise Valuerror: Some points coincide.
159 @example:
161 >>> p = LatLon(53.2611, -0.7972)
163 >>> s = LatLon(53.3206, -1.7297)
164 >>> e = LatLon(53.1887, 0.1334)
165 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58
166 '''
167 p = self.others(start=start)
168 n = self.toNvector()
170 gc, _, _ = self._gc3(p, end, _end_, wrap=wrap)
171 a = gc.cross(n).cross(gc) # along-track point gc × p × gc
172 return _r2m(start.toNvector().angleTo(a, vSign=gc), radius)
174 @deprecated_method
175 def bearingTo(self, other, **unused): # PYCHOK no cover
176 '''DEPRECATED, use method L{initialBearingTo}.
177 '''
178 return self.initialBearingTo(other)
180 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
181 '''Compute the (signed) distance from this point to great circle
182 defined by a start and end point.
184 @arg start: Start point of great circle line (L{LatLon}).
185 @arg end: End point of great circle line (L{LatLon}) or
186 initial bearing from start point (compass
187 C{degrees360}).
188 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
189 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
190 the B{C{start}} and B{C{end}} points (C{bool}).
192 @return: Distance to great circle (C{radians} if C{B{radius}
193 is None} else C{meter}, same units as B{C{radius}}),
194 negative if to the left or positive if to the right
195 of the line .
197 @raise TypeError: If B{C{start}} or B{C{end}} point is not L{LatLon}.
199 @raise Valuerror: Some points coincide.
201 @example:
203 >>> p = LatLon(53.2611, -0.7972)
205 >>> s = LatLon(53.3206, -1.7297)
206 >>> d = p.crossTrackDistanceTo(s, 96) # -305.7
208 >>> e = LatLon(53.1887, 0.1334)
209 >>> d = p.crossTrackDistanceTo(s, e) # -307.5
210 '''
211 p = self.others(start=start)
212 n = self.toNvector()
214 gc, _, _ = self._gc3(p, end, _end_, wrap=wrap)
215 return _r2m(gc.angleTo(n) - PI_2, radius)
217 def destination(self, distance, bearing, radius=R_M, height=None):
218 '''Locate the destination from this point after having
219 travelled the given distance on the given bearing.
221 @arg distance: Distance travelled (C{meter}, same units
222 as B{C{radius}}).
223 @arg bearing: Bearing from this point (compass C{degrees360}).
224 @kwarg radius: Mean earth radius (C{meter}).
225 @kwarg height: Optional height at destination, overriding the
226 default height (C{meter}, same units as B{C{radius}}).
228 @return: Destination point (L{LatLon}).
230 @raise Valuerror: Polar coincidence or invalid B{C{distance}},
231 B{C{bearing}}, B{C{radius}} or B{C{height}}.
233 @example:
235 >>> p = LatLon(51.4778, -0.0015)
236 >>> q = p.destination(7794, 300.7)
237 >>> q.toStr() # 51.513546°N, 000.098345°W
238 '''
239 b = Bearing_(bearing)
240 a = _angular(distance, radius, low=None)
241 sa, ca, sb, cb = sincos2_(a, b)
243 n = self.toNvector()
244 e = NorthPole.cross(n, raiser=_pole_).unit() # east vector at n
245 x = n.cross(e) # north vector at n
246 d = x.times(cb).plus(e.times(sb)) # direction vector @ n
247 n = n.times(ca).plus(d.times(sa))
248 return n.toLatLon(height=height, LatLon=self.classof) # Nvector(n.x, n.y, n.z).toLatLon(...)
250 def distanceTo(self, other, radius=R_M, wrap=False):
251 '''Compute the distance from this to an other point.
253 @arg other: The other point (L{LatLon}).
254 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
255 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
256 the B{C{other}} point (C{bool}).
258 @return: Distance between this and the B{C{other}} point
259 (C{meter}, same units as B{C{radius}} or C{radians}
260 if B{C{radius}} is C{None}).
262 @raise TypeError: Invalid B{C{other}} point.
264 @example:
266 >>> p = LatLon(52.205, 0.119)
267 >>> q = LatLon(48.857, 2.351);
268 >>> d = p.distanceTo(q) # 404.3 Km
269 '''
270 p = self.others(other)
271 if wrap:
272 p = _unrollon(self, p)
273 n = p.toNvector()
274 r = fabs(self.toNvector().angleTo(n, wrap=wrap))
275 return r if radius is None else (Radius(radius) * r)
277# @Property_RO
278# def Ecef(self):
279# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
280# '''
281# return _ALL_MODS.ecef.EcefKarney
283 def _gc3(self, start, end, namend, raiser=_point_, wrap=False):
284 '''(INTERNAL) Return great circle, start and end Nvectors.
285 '''
286 s = start.toNvector()
287 if isscalar(end): # bearing
288 gc = s.greatCircle(end)
289 e = None
290 else: # point
291 p = self.others(end, name=namend)
292 if wrap:
293 p = _unrollon(start, p, wrap=wrap)
294 e = p.toNvector()
295 gc = s.cross(e, raiser=raiser) # XXX .unit()?
296 return gc, s, e
298 def greatCircle(self, bearing):
299 '''Compute the vector normal to great circle obtained by
300 heading on the given bearing from this point.
302 Direction of vector is such that initial bearing vector
303 b = c × n, where n is an n-vector representing this point.
305 @arg bearing: Bearing from this point (compass C{degrees360}).
307 @return: N-vector representing the great circle (L{Nvector}).
308 '''
309 t = Bearing_(bearing)
310 a, b = self.philam
312 sa, ca, sb, cb, st, ct = sincos2_(a, b, t)
313 return Nvector(sb * ct - sa * cb * st,
314 -cb * ct - sa * sb * st,
315 ca * st, name=self.name) # XXX .unit()
317 def greatCircleTo(self, other, wrap=False):
318 '''Compute the vector normal to great circle obtained by
319 heading from this to an other point or on a given bearing.
321 Direction of vector is such that initial bearing vector
322 b = c × n, where n is an n-vector representing this point.
324 @arg other: The other point (L{LatLon}) or the bearing from
325 this point (compass C{degrees360}).
326 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
327 the B{C{other}} point (C{bool}).
329 @return: N-vector representing the great circle (L{Nvector}).
331 @raise TypeError: The B{C{other}} point is not L{LatLon}.
333 @raise Valuerror: Points coincide.
335 @example:
337 >>> p = LatLon(53.3206, -1.7297)
338 >>> gc = p.greatCircle(96.0)
339 >>> gc.toStr() # (-0.79408, 0.12856, 0.59406)
341 >>> q = LatLon(53.1887, 0.1334)
342 >>> g = p.greatCircleTo(q)
343 >>> g.toStr() # (-0.79408, 0.12859, 0.59406)
344 '''
345 gc, _, _ = self._gc3(self, other, _other_, wrap=wrap)
346 return gc.unit()
348 def initialBearingTo(self, other, wrap=False, **unused): # raiser=...
349 '''Compute the initial bearing (forward azimuth) from this
350 to an other point.
352 @arg other: The other point (L{LatLon}).
353 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
354 the B{C{other}} point (C{bool}).
356 @return: Initial bearing (compass C{degrees360}).
358 @raise Crosserror: This point coincides with the B{C{other}}
359 point or the C{NorthPole}, provided
360 L{pygeodesy.crosserrors} is C{True}.
362 @raise TypeError: The B{C{other}} point is not L{LatLon}.
364 @example:
366 >>> p1 = LatLon(52.205, 0.119)
367 >>> p2 = LatLon(48.857, 2.351)
368 >>> b = p1.initialBearingTo(p2) # 156.2
369 '''
370 n = self.toNvector()
371 p = self.others(other)
372 if wrap:
373 p = _unrollon(self, p, wrap=wrap)
374 p = p.toNvector()
375 # see <https://MathForum.org/library/drmath/view/55417.html>
376# gc1 = self.greatCircleTo(other)
377 gc1 = n.cross(p, raiser=_point_) # .unit()
378# gc2 = self.greatCircleTo(NorthPole)
379 gc2 = n.cross(NorthPole, raiser=_pole_) # .unit()
380 return degrees360(gc1.angleTo(gc2, vSign=n))
382 def intermediateChordTo(self, other, fraction, height=None, wrap=False):
383 '''Locate the point projected from the point at given fraction
384 on a straight line (chord) between this and an other point.
386 @arg other: The other point (L{LatLon}).
387 @arg fraction: Fraction between both points (float, between
388 0.0 for this and 1.0 for the other point).
389 @kwarg height: Optional height at the intermediate point,
390 overriding the fractional height (C{meter}).
391 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
392 the B{C{other}} point (C{bool}).
394 @return: Intermediate point (L{LatLon}).
396 @raise TypeError: The B{C{other}} point is not L{LatLon}.
398 @example:
400 >>> p = LatLon(52.205, 0.119)
401 >>> q = LatLon(48.857, 2.351)
402 >>> i = p.intermediateChordTo(q, 0.25) # 51.3723°N, 000.7072°E
403 '''
404 n = self.toNvector()
405 p = self.others(other)
406 if wrap:
407 p = _unrollon(self, p, wrap=wrap)
409 f = Scalar(fraction=fraction)
410 i = p.toNvector().times(f).plus(n.times(1 - f))
411# i = p.toNvector() * f + self.toNvector() * (1 - f))
413 h = self._havg(other, f=f, h=height)
414 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)
416 def intermediateTo(self, other, fraction, height=None, wrap=False):
417 '''Locate the point at a given fraction between this and an
418 other point.
420 @arg other: The other point (L{LatLon}).
421 @arg fraction: Fraction between both points (C{float}, between
422 0.0 for this and 1.0 for the other point).
423 @kwarg height: Optional height at the intermediate point,
424 overriding the fractional height (C{meter}).
425 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
426 the B{C{other}} point (C{bool}).
428 @return: Intermediate point (L{LatLon}).
430 @raise TypeError: The B{C{other}} point is not L{LatLon}.
432 @raise Valuerror: Points coincide or invalid B{C{height}}.
434 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
436 @example:
438 >>> p = LatLon(52.205, 0.119)
439 >>> q = LatLon(48.857, 2.351)
440 >>> i = p.intermediateTo(q, 0.25) # 51.3721°N, 000.7074°E
441 '''
442 n = self.toNvector()
443 p = self.others(other)
444 if wrap:
445 p = _unrollon(self, p, wrap=wrap)
446 p = p.toNvector()
447 f = Scalar(fraction=fraction)
449 x = n.cross(p, raiser=_point_)
450 d = x.unit().cross(n) # unit(n × p) × n
451 # angular distance α, tan(α) = |n × p| / n ⋅ p
452 s, c = sincos2(atan2(x.length, n.dot(p)) * f) # interpolated
453 i = n.times(c).plus(d.times(s)) # n * cosα + d * sinα
455 h = self._havg(other, f=f, h=height)
456 return i.toLatLon(height=h, LatLon=self.classof) # Nvector(i.x, i.y, i.z).toLatLon(...)
458 def intersection(self, end1, start2, end2, height=None, wrap=False):
459 '''Locate the intersection point of two lines each defined
460 by two points or a start point and bearing from North.
462 @arg end1: End point of the first line (L{LatLon}) or the
463 initial bearing at this point (compass C{degrees360}).
464 @arg start2: Start point of the second line (L{LatLon}).
465 @arg end2: End point of the second line (L{LatLon}) or the
466 initial bearing at the second point (compass
467 C{degrees}).
468 @kwarg height: Optional height at the intersection point,
469 overriding the mean height (C{meter}).
470 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll all
471 start and end points (C{bool}).
473 @return: The intersection point (L{LatLon}).
475 @raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}}
476 point is not L{LatLon}.
478 @raise ValueError: Intersection is ambiguous or infinite or
479 the lines are parallel, coincident or null.
481 @see: Function L{sphericalNvector.intersection} and method
482 L{intersection2}.
484 @example:
486 >>> s = LatLon(51.8853, 0.2545)
487 >>> e = LatLon(49.0034, 2.5735)
488 >>> i = s.intersection(108.55, e, 32.44) # 50.9076°N, 004.5086°E
489 '''
490 return intersection(self, end1, start2, end2, height=height,
491 wrap=wrap, LatLon=self.classof)
493 def intersection2(self, end1, start2, end2, height=None, wrap=False):
494 '''Locate the intersections of two (great circle) lines each defined
495 by two points or by a start point and an (initial) bearing.
497 @arg end1: End point of the first line (L{LatLon}) or the
498 initial bearing at this point (compass C{degrees360}).
499 @arg start2: Start point of the second line (L{LatLon}).
500 @arg end2: End point of the second line (L{LatLon}) or the
501 initial bearing at the second start point (compass
502 C{degrees360}).
503 @kwarg height: Optional height at the intersection and antipodal
504 point, overriding the mean height (C{meter}).
505 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
506 B{C{start2}} and both B{C{end*}} points (C{bool}).
508 @return: 2-Tuple C{(intersection, antipode)}, each a B{C{LatLon}}.
510 @raise TypeError: If B{C{start2}}, B{C{end1}} or B{C{end2}}
511 point is not L{LatLon}.
513 @raise ValueError: Intersection is ambiguous or infinite or
514 the lines are parallel, coincident or null.
516 @see: Function L{sphericalNvector.intersection2} and method
517 L{intersection}.
518 '''
519 return intersection2(self, end1, start2, end2, height=height,
520 wrap=wrap, LatLon=self.classof)
522 def isenclosedBy(self, points, wrap=False):
523 '''Check whether a (convex) polygon or composite encloses this point.
525 @arg points: The polygon points or composite (L{LatLon}[],
526 L{BooleanFHP} or L{BooleanGH}).
527 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
528 B{C{points}} (C{bool}).
530 @return: C{True} if this point is inside the polygon or composite,
531 C{False} otherwise.
533 @raise PointsError: Insufficient number of B{C{points}}.
535 @raise TypeError: Some B{C{points}} are not L{LatLon}.
537 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
538 and L{pygeodesy.ispolar} especially if the B{C{points}} may
539 enclose a pole or wrap around the earth I{longitudinally}.
540 '''
541 if _MODS.booleans.isBoolean(points):
542 return points._encloses(self.lat, self.lon, wrap=wrap)
544 # sum subtended angles of each edge (using n0, the
545 # normal vector to this point for sign of α)
546 def _subtangles(ps, w):
547 Ps = self.PointsIter(ps, loop=1, wrap=w)
548 n0 = self.toNvector()
549 _m0 = n0.minus
550 p1 = Ps[0]
551 vs1 = _m0(p1.toNvector())
552 for p2 in Ps.iterate(closed=True):
553 if w and not Ps.looped:
554 p2 = _unrollon(p1, p2)
555 p1 = p2
556 vs2 = _m0(p2.toNvector())
557 yield vs1.angleTo(vs2, vSign=n0) # PYCHOK false
558 vs1 = vs2
560 # Note, this method uses angle summation test: on a plane,
561 # angles for an enclosed point will sum to 360°, angles for
562 # an exterior point will sum to 0°. On a sphere, enclosed
563 # point angles will sum to less than 360° (due to spherical
564 # excess), exterior point angles will be small but non-zero.
565 s = fsum(_subtangles(points, wrap), floats=True) # normal vector
566 # XXX are winding number optimisations equally applicable to
567 # spherical surface?
568 return fabs(s) > PI
570 @deprecated_method
571 def isEnclosedBy(self, points): # PYCHOK no cover
572 '''DEPRECATED, use method C{isenclosedBy}.'''
573 return self.isenclosedBy(points)
575 def iswithin(self, point1, point2, wrap=False):
576 '''Check whether this point is between two other points.
578 If this point is not on the great circle arc defined by
579 both points, return whether it is within the area bound
580 by perpendiculars to the great circle at each point (in
581 the same hemispere).
583 @arg point1: Start point of the arc (L{LatLon}).
584 @arg point2: End point of the arc (L{LatLon}).
585 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
586 B{C{point1}} and B{C{point2}} (C{bool}).
588 @return: C{True} if this point is within the (great circle)
589 arc, C{False} otherwise.
591 @raise TypeError: If B{C{point1}} or B{C{point2}} is not
592 L{LatLon}.
593 '''
594 p1 = self.others(point1=point1)
595 p2 = self.others(point2=point2)
596 if wrap:
597 p1 = _Wrap.point(p1)
598 p2 = _unrollon(p1, p2, wrap=wrap)
599 n, n1, n2 = (_.toNvector() for _ in (self, p1, p2))
601 # corner case, null arc
602 if n1.isequalTo(n2):
603 return n.isequalTo(n1) or n.isequalTo(n2) # PYCHOK returns
605 if n.dot(n1) < 0 or n.dot(n2) < 0: # different hemisphere
606 return False # PYCHOK returns
608 # get vectors representing d0=p0->p1 and d2=p2->p1 and the
609 # dot product d0⋅d2 tells us if p0 is on the p2 side of p1 or
610 # on the other side (similarly for d0=p0->p2 and d1=p1->p2
611 # and dot product d0⋅d1 and p0 on the p1 side of p2 or not)
612 return n.minus(n1).dot(n2.minus(n1)) >= 0 and \
613 n.minus(n2).dot(n1.minus(n2)) >= 0
615 @deprecated_method
616 def isWithin(self, point1, point2): # PYCHOK no cover
617 '''DEPRECATED, use method C{iswithin}.'''
618 return self.iswithin(point1, point2)
620 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
621 '''Find the midpoint between this and an other point.
623 @arg other: The other point (L{LatLon}).
624 @kwarg height: Optional height at the midpoint, overriding
625 the mean height (C{meter}).
626 @kwarg fraction: Midpoint location from this point (C{scalar}),
627 may be negative or greater than 1.0.
628 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
629 B{C{other}} point (C{bool}).
631 @return: Midpoint (L{LatLon}).
633 @raise TypeError: The B{C{other}} point is not L{LatLon}.
635 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
637 @example:
639 >>> p1 = LatLon(52.205, 0.119)
640 >>> p2 = LatLon(48.857, 2.351)
641 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
642 '''
643 if fraction is _0_5:
644 p = self.others(other)
645 if wrap:
646 p = _unrollon(self, p, wrap=wrap)
647 m = self.toNvector().plus(p.toNvector())
648 h = self._havg(other, f=fraction, h=height)
649 r = m.toLatLon(height=h, LatLon=self.classof)
650 else:
651 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
652 return r
654 def nearestOn(self, point1, point2, height=None, within=True, wrap=False):
655 '''Locate the point on the great circle arc between two points
656 closest to this point.
658 @arg point1: Start point of the arc (L{LatLon}).
659 @arg point2: End point of the arc (L{LatLon}).
660 @kwarg height: Optional height, overriding the mean height for
661 the point within the arc (C{meter}), or C{None}
662 to interpolate the height.
663 @kwarg within: If C{True}, return the closest point between both
664 given points, otherwise the closest point
665 elsewhere on the great circle arc (C{bool}).
666 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
667 B{C{point1}} and B{C{point2}} (C{bool}).
669 @return: Closest point on the arc (L{LatLon}).
671 @raise NotImplementedError: Keyword argument C{B{wrap}=True}
672 not supported.
674 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
676 @example:
678 >>> s1 = LatLon(51.0, 1.0)
679 >>> s2 = LatLon(51.0, 2.0)
681 >>> s = LatLon(51.0, 1.9)
682 >>> p = s.nearestOn(s1, s2) # 51.0004°N, 001.9000°E
684 >>> d = p.distanceTo(s) # 42.71 m
686 >>> s = LatLon(51.0, 2.1)
687 >>> p = s.nearestOn(s1, s2) # 51.0000°N, 002.0000°E
688 '''
689 p1 = self.others(point1=point1)
690 p2 = self.others(point2=point2)
691 if wrap:
692 p1 = _Wrap.point(p1)
693 p2 = _unrollon(p1, p2, wrap=wrap)
694 p0 = self
696 if p0.iswithin(p1, p2) and not p1.isequalTo(p2, EPS):
697 # closer to arc than to its endpoints,
698 # find the closest point on the arc
699 gc1 = p1.toNvector().cross(p2.toNvector())
700 gc2 = p0.toNvector().cross(gc1)
701 n = gc1.cross(gc2)
703 elif within: # for backward compatibility, XXX unwrapped
704 return point1 if (self.distanceTo(point1) <
705 self.distanceTo(point2)) else point2
707 else: # handle beyond arc extent by .vector3d.nearestOn
708 n1 = p1.toNvector()
709 n2 = p2.toNvector()
710 n = p0.toNvector().nearestOn(n1, n2, within=False)
711 if n is n1:
712 return p1 # is point1
713 elif n is n2:
714 return p2 # is point2 if not wrap
716 p = n.toLatLon(height=height or 0, LatLon=self.classof)
717 if height in (None, False): # interpolate height within extent
718 d = p1.distanceTo(p2)
719 f = (p1.distanceTo(p) / d) if d > EPS0 else _0_5
720 p.height = p1._havg(p2, f=max(_0_0, min(f, _1_0)))
721 return p
723 # @deprecated_method
724 def nearestOn2(self, points, **closed_radius_height): # PYCHOK no cover
725 '''DEPRECATED, use method L{sphericalNvector.LatLon.nearestOn3}.
727 @return: ... 2-Tuple C{(closest, distance)} of the C{closest}
728 point (L{LatLon}) on the polygon and the C{distance}
729 to that point from this point ...
730 '''
731 r = self.nearestOn3(points, **closed_radius_height)
732 return r.closest, r.distance
734 def nearestOn3(self, points, closed=False, radius=R_M, height=None, wrap=False):
735 '''Locate the point on a path or polygon (with great circle
736 arcs joining consecutive points) closest to this point.
738 The closest point is either on within the extent of any great
739 circle arc or the nearest of the arc's end points.
741 @arg points: The path or polygon points (L{LatLon}[]).
742 @kwarg closed: Optionally, close the polygon (C{bool}).
743 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
744 @kwarg height: Optional height, overriding the mean height
745 for a point within the arc (C{meter}).
746 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
747 the B{C{points}} (C{bool}).
749 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of
750 the C{closest} point (L{LatLon}), the C{distance}
751 between this and the C{closest} point in C{meter},
752 same units as B{C{radius}} or in C{radians} if
753 B{C{radius}} is C{None} and the C{angle} from this
754 to the C{closest} point in compass C{degrees360}.
756 @raise TypeError: Some B{C{points}} are not C{LatLon}.
758 @raise ValueError: No B{C{points}}.
759 '''
760 Ps = self.PointsIter(points, loop=1, wrap=wrap)
761 _r = self.distanceTo
762 _n = self.nearestOn
764 c = p1 = Ps[0]
765 r = _r(c, radius=None) # radians
766 for p2 in Ps.iterate(closed=closed):
767 if wrap and not Ps.looped:
768 p2 = _unrollon(p1, p2)
769 p = _n(p1, p2, height=height)
770 d = _r(p, radius=None) # radians
771 if d < r:
772 c, r = p, d
773 p1 = p2
774 d = r if radius is None else (Radius(radius) * r)
775 return NearestOn3Tuple(c, d, degrees360(r))
777 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian, datum=None
778 '''Convert this point to C{Nvector}-based cartesian (ECEF) coordinates.
780 @kwarg Cartesian_and_kwds: Optional L{Cartesian} and L{Cartesian} keyword
781 arguments, like C{datum}. Use C{B{Cartesian}=...}
782 to override this L{Cartesian} class or specify
783 C{B{Cartesian}=None}.
785 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} is
786 set to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
787 C, M, datum)} with C{C} and C{M} if available.
789 @raise TypeError: Invalid B{C{Cartesian_and_kwds}} argument.
790 '''
791 kwds = _xkwds(Cartesian_and_kwds, Cartesian=Cartesian, datum=self.datum)
792 return LatLonSphericalBase.toCartesian(self, **kwds)
794 def toNvector(self, **Nvector_and_kwds): # PYCHOK signature
795 '''Convert this point to L{Nvector} components, I{including height}.
797 @kwarg Nvector_and_kwds: Optional L{Nvector} and L{Nvector} keyword
798 arguments. Specify C{B{Nvector}=...} to
799 override this L{Nvector} class or use
800 C{B{Nvector}=None}.
802 @return: The C{n-vector} components (L{Nvector}) or if B{C{Nvector}} is
803 set to C{None}, a L{Vector4Tuple}C{(x, y, z, h)}.
805 @raise TypeError: Invalid B{C{Nvector_and_kwds}} argument.
807 @example:
809 >>> p = LatLon(45, 45)
810 >>> n = p.toNvector()
811 >>> n.toStr() # [0.50, 0.50, 0.70710]
812 '''
813 kwds = _xkwds(Nvector_and_kwds, Nvector=Nvector)
814 return LatLonNvectorBase.toNvector(self, **kwds)
817class Nvector(NvectorBase):
818 '''An n-vector is a position representation using a (unit) vector
819 normal to the earth's surface. Unlike lat-/longitude points,
820 n-vectors have no singularities or discontinuities.
822 For many applications, n-vectors are more convenient to work
823 with than other position representations like lat-/longitude,
824 earth-centred earth-fixed (ECEF) vectors, UTM coordinates, etc.
826 On a spherical model earth, an n-vector is equivalent to an
827 earth-centred earth-fixed (ECEF) vector.
829 Note commonality with L{ellipsoidalNvector.Nvector}.
830 '''
831 _datum = Datums.Sphere # default datum (L{Datum})
833 def toCartesian(self, **Cartesian_and_kwds): # PYCHOK Cartesian=Cartesian
834 '''Convert this n-vector to C{Nvector}-based cartesian
835 (ECEF) coordinates.
837 @kwarg Cartesian_and_kwds: Optional L{Cartesian} and L{Cartesian} keyword
838 arguments, like C{h}. Use C{B{Cartesian}=...}
839 to override this L{Cartesian} class or specify
840 C{B{Cartesian}=None}.
842 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} is
843 set to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
844 C, M, datum)} with C{C} and C{M} if available.
846 @raise TypeError: Invalid B{C{Cartesian_and_kwds}} argument.
847 '''
848 kwds = _xkwds(Cartesian_and_kwds, h=self.h, Cartesian=Cartesian)
849 return NvectorBase.toCartesian(self, **kwds) # class or .classof
851 def toLatLon(self, **LatLon_and_kwds): # PYCHOK height=None, LatLon=LatLon
852 '''Convert this n-vector to an C{Nvector}-based geodetic point.
854 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
855 arguments, like C{height}. Use C{B{LatLon}=...}
856 to override this L{LatLon} class or specify
857 C{B{LatLon}=None}.
859 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is set
860 to C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
861 C, M, datum)} with C{C} and C{M} if available.
863 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
865 @raise ValueError: Invalid B{C{height}}.
866 '''
867 kwds = _xkwds(LatLon_and_kwds, height=self.h, LatLon=LatLon)
868 return NvectorBase.toLatLon(self, **kwds) # class or .classof
870 def greatCircle(self, bearing):
871 '''Compute the n-vector normal to great circle obtained by
872 heading on given compass bearing from this point as its
873 n-vector.
875 Direction of vector is such that initial bearing vector
876 b = c × p.
878 @arg bearing: Initial compass bearing (C{degrees}).
880 @return: N-vector representing great circle (L{Nvector}).
882 @raise Valuerror: Polar coincidence.
884 @example:
886 >>> n = LatLon(53.3206, -1.7297).toNvector()
887 >>> gc = n.greatCircle(96.0) # [-0.794, 0.129, 0.594]
888 '''
889 s, c = sincos2d(Bearing(bearing))
891 e = NorthPole.cross(self, raiser=_pole_) # easting
892 n = self.cross(e, raiser=_point_) # northing
894 e = e.times(c / e.length)
895 n = n.times(s / n.length)
896 return n.minus(e)
899_Nvll = LatLon(_0_0, _0_0, name=_Nv00_) # reference instance (L{LatLon})
902def areaOf(points, radius=R_M, wrap=False):
903 '''Calculate the area of a (spherical) polygon or composite (with
904 great circle arcs joining consecutive points).
906 @arg points: The polygon points or clips (C{LatLon}[],
907 L{BooleanFHP} or L{BooleanGH}).
908 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
909 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
910 B{C{points}} (C{bool}).
912 @return: Polygon area (C{meter} I{squared} , same units as
913 B{C{radius}}, or C{radians} if B{C{radius}} is C{None}).
915 @raise PointsError: Insufficient number of B{C{points}}.
917 @raise TypeError: Some B{C{points}} are not L{LatLon}.
919 @see: Functions L{pygeodesy.areaOf}, L{sphericalTrigonometry.areaOf}
920 and L{ellipsoidalKarney.areaOf}.
922 @example:
924 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
925 >>> areaOf(b) # 8666058750.718977
926 '''
927 def _interangles(ps, w): # like .karney._polygon
928 Ps = _Nvll.PointsIter(ps, loop=2, wrap=w)
929 # use vector to 1st point as plane normal for sign of α
930 n0 = Ps[0].toNvector()
932 v2 = Ps[0]._N_vector # XXX v2 == no?
933 p1 = Ps[1]
934 v1 = p1._N_vector
935 gc = v2.cross(v1)
936 for p2 in Ps.iterate(closed=True):
937 if w and not Ps.looped:
938 p2 = _unrollon(p1, p2)
939 p1 = p2
940 v2 = p2._N_vector
941 gc1 = v1.cross(v2)
942 v1 = v2
943 yield gc.angleTo(gc1, vSign=n0)
944 gc = gc1
946 if _MODS.booleans.isBoolean(points):
947 r = points._sum2(LatLon, areaOf, radius=None, wrap=wrap)
948 else:
949 # sum interior angles: depending on whether polygon is cw or ccw,
950 # angle between edges is π−α or π+α, where α is angle between
951 # great-circle vectors; so sum α, then take n·π − |Σα| (cannot
952 # use Σ(π−|α|) as concave polygons would fail)
953 s = fsum(_interangles(points, wrap), floats=True)
954 # using Girard’s theorem: A = [Σθᵢ − (n−2)·π]·R²
955 # (PI2 - abs(s) == (n*PI - abs(s)) - (n-2)*PI)
956 r = fabs(PI2 - fabs(s))
957 return r if radius is None else (r * Radius(radius)**2)
960def intersecant2(center, circle, point, bearing, radius=R_M, exact=False,
961 height=None, wrap=False): # was=True
962 '''Compute the intersections of a circle and a line.
964 @arg center: Center of the circle (L{LatLon}).
965 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
966 or a point on the circle (L{LatLon}).
967 @arg point: A point in- or outside the circle (L{LatLon}).
968 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or
969 a second point on the line (L{LatLon}).
970 @kwarg radius: Mean earth radius (C{meter}, conventionally).
971 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
972 destination and distance, if C{False} use the basic
973 rhumb methods (C{bool}) or if C{None} use the I{great
974 circle} methods.
975 @kwarg height: Optional height for the intersection points (C{meter},
976 conventionally) or C{None}.
977 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}}
978 and the B{C{circle}} and B{C{bearing}} points (C{bool}).
980 @return: 2-Tuple of the intersection points (representing a chord), each
981 an instance of this class. For a tangent line, both points are
982 the same instance, the B{C{point}} or wrapped or I{normalized}.
984 @raise IntersectionError: The circle and line do not intersect.
986 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
987 B{C{circle}} or B{C{bearing}} invalid.
989 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
990 B{C{exact}} or B{C{height}}.
991 '''
992 c = _Nvll.others(center=center)
993 p = _Nvll.others(point=point)
994 try:
995 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact,
996 height=height, wrap=wrap)
997 except (TypeError, ValueError) as x:
998 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing,
999 exact=exact, wrap=wrap)
1002def intersection(start1, end1, start2, end2, height=None, wrap=False,
1003 LatLon=LatLon, **LatLon_kwds):
1004 '''Locate the intersections of two (great circle) lines each defined
1005 by two points or by a start point and an (initial) bearing.
1007 @arg start1: Start point of the first line (L{LatLon}).
1008 @arg end1: End point of the first line (L{LatLon}) or the initial
1009 bearing at the first start point (compass C{degrees360}).
1010 @arg start2: Start point of the second line (L{LatLon}).
1011 @arg end2: End point of the second line (L{LatLon}) or the initial
1012 bearing at the second start point (compass C{degrees360}).
1013 @kwarg height: Optional height at the intersection point,
1014 overriding the mean height (C{meter}).
1015 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{start2}}
1016 and both B{C{end*}} points (C{bool}).
1017 @kwarg LatLon: Optional class to return the intersection point
1018 (L{LatLon}).
1019 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1020 arguments, ignored if C{B{LatLon} is None}.
1022 @return: The intersection point (B{C{LatLon}}) or if C{B{LatLon}
1023 is None}, a cartesian L{Ecef9Tuple}C{(x, y, z, lat, lon,
1024 height, C, M, datum)} with C{C} and C{M} if available.
1026 @raise TypeError: If B{C{start*}} or B{C{end*}} is not L{LatLon}.
1028 @raise ValueError: Intersection is ambiguous or infinite or
1029 the lines are parallel, coincident or null.
1031 @see: Function L{sphericalNvector.intersection2}.
1033 @example:
1035 >>> p = LatLon(51.8853, 0.2545)
1036 >>> q = LatLon(49.0034, 2.5735)
1037 >>> i = intersection(p, 108.55, q, 32.44) # 50.9076°N, 004.5086°E
1038 '''
1039 i, _, h = _intersect3(start1, end1, start2, end2, height, wrap)
1040 kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)
1041 return i.toLatLon(**kwds)
1044def intersection2(start1, end1, start2, end2, height=None, wrap=False,
1045 LatLon=LatLon, **LatLon_kwds):
1046 '''Locate the intersections of two (great circle) lines each defined
1047 by two points or by a start point and an (initial) bearing.
1049 @arg start1: Start point of the first line (L{LatLon}).
1050 @arg end1: End point of the first line (L{LatLon}) or the
1051 initial bearing at the first start point
1052 (compass C{degrees360}).
1053 @arg start2: Start point of the second line (L{LatLon}).
1054 @arg end2: End point of the second line (L{LatLon}) or the
1055 initial bearing at the second start point
1056 (compass C{degrees360}).
1057 @kwarg height: Optional height at the intersection and antipodal
1058 point, overriding the mean height (C{meter}).
1059 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{start2}}
1060 and both B{C{end*}} points (C{bool}).
1061 @kwarg LatLon: Optional class to return the intersection and
1062 antipodal points (L{LatLon}).
1063 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1064 arguments, ignored if C{B{LatLon} is None}.
1066 @return: 2-Tuple C{(intersection, antipode)}, each a (B{C{LatLon}})
1067 or if C{B{LatLon} is None}, a cartesian L{Ecef9Tuple}C{(x,
1068 y, z, lat, lon, height, C, M, datum)} with C{C} and C{M}
1069 if available.
1071 @raise TypeError: If B{C{start*}} or B{C{end*}} is not L{LatLon}.
1073 @raise ValueError: Intersection is ambiguous or infinite or
1074 the lines are parallel, coincident or null.
1076 @see: Function L{sphericalNvector.intersection}.
1077 '''
1078 i, a, h = _intersect3(start1, end1, start2, end2, height, wrap)
1079 kwds = _xkwds(LatLon_kwds, height=h, LatLon=LatLon)
1080 return i.toLatLon(**kwds), a.toLatLon(**kwds)
1083def _intersect3(start1, end1, start2, end2, height, wrap):
1084 '''(INTERNAL) Return the intersection and antipodal points for
1085 functions C{intersection} and C{intersection2}.
1086 '''
1087 p1 = _Nvll.others(start1=start1)
1088 p2 = _Nvll.others(start2=start2)
1089 if wrap:
1090 p2 = _unrollon(p1, p2, wrap=wrap)
1091 # If gc1 and gc2 are great circles through start and end points
1092 # (or defined by start point and bearing), then the candidate
1093 # intersections are simply gc1 × gc2 and gc2 × gc1. Most of the
1094 # work is deciding the correct intersection point to select! If
1095 # bearing is given, that determines the intersection, but if both
1096 # lines are defined by start/end points, take closer intersection.
1097 gc1, s1, e1 = _Nvll._gc3(p1, end1, 'end1', wrap=wrap)
1098 gc2, s2, e2 = _Nvll._gc3(p2, end2, 'end2', wrap=wrap)
1100 hs = start1.height, start2.height
1101 # there are two (antipodal) candidate intersection
1102 # points ... we have to choose the one to return
1103 i1 = gc1.cross(gc2, raiser=_lines_)
1104 i2 = gc2.cross(gc1, raiser=_lines_)
1106 # selection of intersection point depends on how
1107 # lines are defined (by bearings or endpoints)
1108 if e1 and e2: # endpoint+endpoint
1109 d = sumOf((s1, s2, e1, e2)).dot(i1)
1110 hs += end1.height, end2.height
1111 elif e1 and not e2: # endpoint+bearing
1112 # gc2 x v2 . i1 +ve means v2 bearing points to i1
1113 d = gc2.cross(s2).dot(i1)
1114 hs += end1.height,
1115 elif e2 and not e1: # bearing+endpoint
1116 # gc1 x v1 . i1 +ve means v1 bearing points to i1
1117 d = gc1.cross(s1).dot(i1)
1118 hs += end2.height,
1119 else: # bearing+bearing
1120 # if gc x v . i1 is +ve, initial bearing is
1121 # towards i1, otherwise towards antipodal i2
1122 d1 = gc1.cross(s1).dot(i1) # +ve means p1 bearing points to i1
1123 d2 = gc2.cross(s2).dot(i1) # +ve means p2 bearing points to i1
1124 if d1 > 0 and d2 > 0:
1125 d = 1 # both point to i1
1126 elif d1 < 0 and d2 < 0:
1127 d = -1 # both point to i2
1128 else: # d1, d2 opposite signs
1129 # intersection is at further-away intersection point,
1130 # take opposite intersection from mid- point of v1
1131 # and v2 [is this always true?] XXX changed to always
1132 # get intersection p1 bearing points to, aka being
1133 # located "after" p1 along the bearing at p1, like
1134 # function .sphericalTrigonometry._intersect and
1135 # .ellipsoidalBaseDI._intersect3
1136 d = d1 # neg(s1.plus(s2).dot(i1))
1138 h = fmean(hs) if height is None else height
1139 return (i1, i2, h) if d > 0 else (i2, i1, h)
1142def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):
1143 '''Compute the I{geographic} mean of the supplied points.
1145 @arg points: Array of points to be averaged (L{LatLon}[]).
1146 @kwarg height: Optional height, overriding the mean height
1147 (C{meter}).
1148 @kwarg wrap: If C{True}, wrap or I{normalize} B{C{points}} (C{bool}).
1149 @kwarg LatLon: Optional class to return the mean point (L{LatLon}).
1150 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1151 arguments, ignored if C{B{LatLon} is None}.
1153 @return: Point at geographic mean and mean height (B{C{LatLon}}).
1155 @raise PointsError: Insufficient number of B{C{points}} or
1156 some B{C{points}} are not C{LatLon}.
1157 '''
1158 def _N_vs(ps, w):
1159 Ps = _Nvll.PointsIter(ps, wrap=w)
1160 for p in Ps.iterate(closed=False):
1161 yield p._N_vector
1163 try:
1164 # geographic mean
1165 n = _nsumOf(_N_vs(points, wrap), height, Nvector, {})
1166 except (TypeError, ValueError) as x:
1167 raise PointsError(points=points, wrap=wrap, LatLon=LatLon, cause=x)
1168 return n.toLatLon(**_xkwds(LatLon_kwds, LatLon=LatLon, height=n.h,
1169 name=meanOf.__name__))
1172@deprecated_function
1173def nearestOn2(point, points, **closed_radius_height): # PYCHOK no cover
1174 '''DEPRECATED, use method L{sphericalNvector.nearestOn3}.
1176 @return: ... 2-Tuple C{(closest, distance)} of the C{closest}
1177 point (L{LatLon}) on the polygon and the C{distance}
1178 between the C{closest} and the given B{C{point}} ...
1179 '''
1180 r = nearestOn3(point, points, **closed_radius_height)
1181 return r.closest, r.distance
1184def nearestOn3(point, points, closed=False, radius=R_M, height=None, wrap=False):
1185 '''Locate the point on a polygon (with great circle arcs joining
1186 consecutive points) closest to an other point.
1188 If the given point is between the end points of a great circle
1189 arc, the closest point is on that arc. Otherwise, the closest
1190 point is the nearest of the arc's end points.
1192 @arg point: The other, reference point (L{LatLon}).
1193 @arg points: The polygon points (L{LatLon}[]).
1194 @kwarg closed: Optionally, close the polygon (C{bool}).
1195 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1196 @kwarg height: Optional height, overriding the mean height for
1197 a point within the (great circle) arc (C{meter}).
1198 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1199 B{C{points}} (C{bool}).
1201 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of
1202 the C{closest} point (L{LatLon}) on the polygon, the
1203 C{distance} and the C{angle} between the C{closest}
1204 and the given B{C{point}}. The C{distance} is in
1205 C{meter}, same units as B{C{radius}} or in C{radians}
1206 if B{C{radius}} is C{None}, the C{angle} is in compass
1207 C{degrees360}.
1209 @raise PointsError: Insufficient number of B{C{points}}.
1211 @raise TypeError: Some B{C{points}} or B{C{point}} not C{LatLon}.
1212 '''
1213 _xinstanceof(LatLon, point=point)
1215 return point.nearestOn3(points, closed=closed, radius=radius,
1216 height=height, wrap=wrap)
1219def perimeterOf(points, closed=False, radius=R_M, wrap=False):
1220 '''Compute the perimeter of a (spherical) polygon or composite (with
1221 great circle arcs joining consecutive points).
1223 @arg points: The polygon points (L{LatLon}[]).
1224 @kwarg closed: Optionally, close the polygon (C{bool}).
1225 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1226 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1227 B{C{points}} (C{bool}).
1229 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1230 or C{radians} if B{C{radius}} is C{None}).
1232 @raise PointsError: Insufficient number of B{C{points}}.
1234 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1236 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1237 C{B{points}} a composite.
1239 @see: Functions L{pygeodesy.perimeterOf}, L{ellipsoidalKarney.perimeterOf}
1240 and L{sphericalTrigonometry.perimeterOf}.
1241 '''
1242 def _rads(ps, c, w): # angular edge lengths in radians
1243 Ps = _Nvll.PointsIter(ps, loop=1, wrap=w)
1244 p1 = Ps[0]
1245 v1 = p1._N_vector
1246 for p2 in Ps.iterate(closed=c):
1247 if w and not (c and Ps.looped):
1248 p2 = _unrollon(p1, p2)
1249 p1 = p2
1250 v2 = p2._N_vector
1251 yield v1.angleTo(v2)
1252 v1 = v2
1254 if _MODS.booleans.isBoolean(points):
1255 if not closed:
1256 notImplemented(None, closed=closed, points=_composite_)
1257 r = points._sum2(LatLon, perimeterOf, closed=True, radius=None, wrap=wrap)
1258 else:
1259 r = fsum(_rads(points, closed, wrap), floats=True)
1260 return r if radius is None else (Radius(radius) * r)
1263def sumOf(nvectors, Vector=Nvector, h=None, **Vector_kwds):
1264 '''Return the I{vectorial} sum of two or more n-vectors.
1266 @arg nvectors: Vectors to be added (L{Nvector}[]).
1267 @kwarg Vector: Optional class for the vectorial sum (L{Nvector}).
1268 @kwarg h: Optional height, overriding the mean height (C{meter}).
1269 @kwarg Vector_kwds: Optional, additional B{C{Vector}} keyword arguments.
1271 @return: Vectorial sum (B{C{Vector}}).
1273 @raise VectorError: No B{C{nvectors}}.
1274 '''
1275 try:
1276 return _nsumOf(nvectors, h, Vector, Vector_kwds)
1277 except (TypeError, ValueError) as x:
1278 raise VectorError(nvectors=nvectors, Vector=Vector, cause=x)
1281def triangulate(point1, bearing1, point2, bearing2,
1282 height=None, wrap=False,
1283 LatLon=LatLon, **LatLon_kwds):
1284 '''Locate a point given two known points and the initial bearings
1285 from those points.
1287 @arg point1: First reference point (L{LatLon}).
1288 @arg bearing1: Bearing at the first point (compass C{degrees360}).
1289 @arg point2: Second reference point (L{LatLon}).
1290 @arg bearing2: Bearing at the second point (compass C{degrees360}).
1291 @kwarg height: Optional height at the triangulated point, overriding
1292 the mean height (C{meter}).
1293 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{point2}}
1294 (C{bool}).
1295 @kwarg LatLon: Optional class to return the triangulated point (L{LatLon}).
1296 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1297 arguments, ignored if C{B{LatLon} is None}.
1299 @return: Triangulated point (B{C{LatLon}}).
1301 @raise TypeError: If B{C{point1}} or B{C{point2}} is not L{LatLon}.
1303 @raise Valuerror: Points coincide.
1305 @example:
1307 >>> p = LatLon("47°18.228'N","002°34.326'W") # Basse Castouillet
1308 >>> q = LatLon("47°18.664'N","002°31.717'W") # Basse Hergo
1309 >>> t = triangulate(p, 7, q, 295) # 47.323667°N, 002.568501°W'
1310 '''
1311 return _triangulate(_Nvll.others(point1=point1), bearing1,
1312 _Nvll.others(point2=point2), bearing2,
1313 height=height, wrap=wrap,
1314 LatLon=LatLon, **LatLon_kwds)
1317def trilaterate(point1, distance1, point2, distance2, point3, distance3, # PYCHOK args
1318 radius=R_M, height=None, useZ=False, wrap=False,
1319 LatLon=LatLon, **LatLon_kwds):
1320 '''Locate a point at given distances from three other points.
1322 @arg point1: First point (L{LatLon}).
1323 @arg distance1: Distance to the first point (C{meter}, same units
1324 as B{C{radius}}).
1325 @arg point2: Second point (L{LatLon}).
1326 @arg distance2: Distance to the second point (C{meter}, same units
1327 as B{C{radius}}).
1328 @arg point3: Third point (L{LatLon}).
1329 @arg distance3: Distance to the third point (C{meter}, same units
1330 as B{C{radius}}).
1331 @kwarg radius: Mean earth radius (C{meter}).
1332 @kwarg height: Optional height at the trilaterated point, overriding
1333 the IDW height (C{meter}, same units as B{C{radius}}).
1334 @kwarg useZ: Include Z component iff non-NaN, non-zero (C{bool}).
1335 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{point2}}
1336 and B{C{point3}} (C{bool}).
1337 @kwarg LatLon: Optional class to return the trilaterated point (L{LatLon}).
1338 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments,
1339 ignored if C{B{LatLon} is None}.
1341 @return: Trilaterated point (B{C{LatLon}}).
1343 @raise IntersectionError: No intersection, trilateration failed.
1345 @raise TypeError: Invalid B{C{point1}}, B{C{point2}} or B{C{point3}}.
1347 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}},
1348 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.
1350 @see: U{Trilateration<https://WikiPedia.org/wiki/Trilateration>}.
1351 '''
1352 return _trilaterate(_Nvll.others(point1=point1), distance1,
1353 _Nvll.others(point2=point2), distance2,
1354 _Nvll.others(point3=point3), distance3,
1355 radius=radius, height=height, useZ=useZ,
1356 wrap=wrap, LatLon=LatLon, **LatLon_kwds)
1359__all__ += _ALL_OTHER(Cartesian, LatLon, Nvector, # classes
1360 areaOf, # functions
1361 intersecant2, intersection, intersection2, ispolar,
1362 meanOf,
1363 nearestOn2, nearestOn3,
1364 perimeterOf,
1365 sumOf,
1366 triangulate, trilaterate)
1368# **) MIT License
1369#
1370# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1371#
1372# Permission is hereby granted, free of charge, to any person obtaining a
1373# copy of this software and associated documentation files (the "Software"),
1374# to deal in the Software without restriction, including without limitation
1375# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1376# and/or sell copies of the Software, and to permit persons to whom the
1377# Software is furnished to do so, subject to the following conditions:
1378#
1379# The above copyright notice and this permission notice shall be included
1380# in all copies or substantial portions of the Software.
1381#
1382# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1383# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1384# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1385# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1386# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1387# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1388# OTHER DEALINGS IN THE SOFTWARE.