Coverage for pygeodesy/sphericalTrigonometry.py: 93%
393 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-12-29 12:35 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2023-12-29 12:35 -0500
2# -*- coding: utf-8 -*-
4u'''Spherical, C{trigonometry}-based geodesy.
6Trigonometric classes geodetic (lat-/longitude) L{LatLon} and
7geocentric (ECEF) L{Cartesian} and functions L{areaOf}, L{intersection},
8L{intersections2}, L{isPoleEnclosedBy}, L{meanOf}, L{nearestOn3} and
9L{perimeterOf}, I{all spherical}.
11Pure Python implementation of geodetic (lat-/longitude) methods using
12spherical trigonometry, transcoded from JavaScript originals by
13I{(C) Chris Veness 2011-2016} published under the same MIT Licence**, see
14U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}.
15'''
16# make sure int/int division yields float quotient, see .basics
17from __future__ import division as _; del _ # PYCHOK semicolon
19from pygeodesy.basics import copysign0, map1, signOf
20from pygeodesy.constants import EPS, EPS1, EPS4, PI, PI2, PI_2, PI_4, R_M, \
21 isnear0, isnear1, isnon0, _0_0, _0_5, \
22 _1_0, _2_0, _90_0
23from pygeodesy.datums import _ellipsoidal_datum, _mean_radius
24from pygeodesy.errors import _AssertionError, CrossError, crosserrors, \
25 _TypeError, _ValueError, IntersectionError, \
26 _xError, _xkwds, _xkwds_get, _xkwds_pop
27from pygeodesy.fmath import favg, fdot, fmean, hypot
28from pygeodesy.fsums import Fsum, fsum, fsumf_
29from pygeodesy.formy import antipode_, bearing_, _bearingTo2, excessAbc_, \
30 excessGirard_, excessLHuilier_, opposing_, _radical2, \
31 vincentys_
32from pygeodesy.interns import _1_, _2_, _coincident_, _composite_, _colinear_, \
33 _concentric_, _convex_, _end_, _infinite_, _invalid_,\
34 _line_, _near_, _not_, _null_, _parallel_, _point_, \
35 _SPACE_, _too_
36from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER
37# from pygeodesy.named import notImplemented # from .points
38# from pygeodesy.nvectorBase import NvectorBase, sumOf # _MODE
39from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, NearestOn3Tuple, \
40 Triangle7Tuple, Triangle8Tuple
41from pygeodesy.points import ispolar, nearestOn5 as _nearestOn5, \
42 Fmt as _Fmt, notImplemented # XXX shadowed
43from pygeodesy.props import deprecated_function, deprecated_method
44from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \
45 _intersecant2, LatLonSphericalBase, \
46 _rads3, _radians2m, _trilaterate5
47# from pygeodesy.streprs import Fmt as _Fmt # from .points XXX shadowed
48from pygeodesy.units import Bearing_, Height, _isDegrees, _isRadius, Lam_, \
49 Phi_, Radius_, Scalar
50from pygeodesy.utily import acos1, asin1, atan1d, atan2d, degrees90, degrees180, \
51 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \
52 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI
53from pygeodesy.vector3d import sumOf, Vector3d
55from math import asin, atan2, cos, degrees, fabs, radians, sin
57__all__ = _ALL_LAZY.sphericalTrigonometry
58__version__ = '23.12.18'
60_PI_EPS4 = PI - EPS4
61if _PI_EPS4 >= PI:
62 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4)
65class Cartesian(CartesianSphericalBase):
66 '''Extended to convert geocentric, L{Cartesian} points to
67 spherical, geodetic L{LatLon}.
68 '''
70 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon
71 '''Convert this cartesian point to a C{spherical} geodetic point.
73 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
74 arguments. Use C{B{LatLon}=...} to override
75 this L{LatLon} class or specify C{B{LatLon}=None}.
77 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is C{None},
78 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
79 with C{C} and C{M} if available.
81 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
82 '''
83 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
84 return CartesianSphericalBase.toLatLon(self, **kwds)
87class LatLon(LatLonSphericalBase):
88 '''New point on a spherical earth model, based on trigonometry formulae.
89 '''
91 def _ab1_ab2_db5(self, other, wrap):
92 '''(INTERNAL) Helper for several methods.
93 '''
94 a1, b1 = self.philam
95 a2, b2 = self.others(other, up=2).philam
96 if wrap:
97 a2, b2 = _Wrap.philam(a2, b2)
98 db, b2 = unrollPI(b1, b2, wrap=wrap)
99 else: # unrollPI shortcut
100 db = b2 - b1
101 return a1, b1, a2, b2, db
103 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
104 '''Compute the (signed) distance from the start to the closest
105 point on the great circle line defined by a start and an
106 end point.
108 That is, if a perpendicular is drawn from this point to the
109 great circle line, the along-track distance is the distance
110 from the start point to the point where the perpendicular
111 crosses the line.
113 @arg start: Start point of the great circle line (L{LatLon}).
114 @arg end: End point of the great circle line (L{LatLon}).
115 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
116 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
117 the B{C{start}} and B{C{end}} point (C{bool}).
119 @return: Distance along the great circle line (C{radians}
120 if C{B{radius} is None} or C{meter}, same units
121 as B{C{radius}}), positive if I{after} the
122 B{C{start}} toward the B{C{end}} point of the
123 line, I{negative} if before or C{0} if at the
124 B{C{start}} point.
126 @raise TypeError: Invalid B{C{start}} or B{C{end}} point.
128 @raise ValueError: Invalid B{C{radius}}.
129 '''
130 r, x, b = self._a_x_b3(start, end, radius, wrap)
131 cx = cos(x)
132 return _0_0 if isnear0(cx) else \
133 _radians2m(copysign0(acos1(cos(r) / cx), cos(b)), radius)
135 def _a_x_b3(self, start, end, radius, wrap):
136 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo.
137 '''
138 s = self.others(start=start)
139 e = self.others(end=end)
140 s, e, w = _unrollon3(self, s, e, wrap)
142 r = Radius_(radius)
143 r = s.distanceTo(self, r, wrap=w) / r
145 b = radians(s.initialBearingTo(self, wrap=w)
146 - s.initialBearingTo(e, wrap=w))
147 x = asin(sin(r) * sin(b))
148 return r, x, -b
150 @deprecated_method
151 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover
152 '''DEPRECATED, use method L{initialBearingTo}.
153 '''
154 return self.initialBearingTo(other, wrap=wrap, raiser=raiser)
156 def crossingParallels(self, other, lat, wrap=False):
157 '''Return the pair of meridians at which a great circle defined
158 by this and an other point crosses the given latitude.
160 @arg other: The other point defining great circle (L{LatLon}).
161 @arg lat: Latitude at the crossing (C{degrees}).
162 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
163 B{C{other}} point (C{bool}).
165 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
166 C{None} if the great circle doesn't reach B{C{lat}}.
167 '''
168 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
169 sa, ca, sa1, ca1, \
170 sa2, ca2, sdb, cdb = sincos2_(radians(lat), a1, a2, db)
171 sa1 *= ca2 * ca
173 x = sa1 * sdb
174 y = sa1 * cdb - ca1 * sa2 * ca
175 z = ca1 * sdb * ca2 * sa
177 h = hypot(x, y)
178 if h < EPS or fabs(z) > h: # PYCHOK no cover
179 return None # great circle doesn't reach latitude
181 m = atan2(-y, x) + b1 # longitude at max latitude
182 d = acos1(z / h) # delta longitude to intersections
183 return degrees180(m - d), degrees180(m + d)
185 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
186 '''Compute the (signed) distance from this point to
187 the great circle defined by a start and an end point.
189 @arg start: Start point of the great circle line (L{LatLon}).
190 @arg end: End point of the great circle line (L{LatLon}).
191 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
192 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
193 the B{C{start}} and B{C{end}} point (C{bool}).
195 @return: Distance to the great circle (C{radians} if
196 B{C{radius}} or C{meter}, same units as
197 B{C{radius}}), I{negative} if to the left or
198 I{positive} if to the right of the line.
200 @raise TypeError: If B{C{start}} or B{C{end}} is not L{LatLon}.
202 @raise ValueError: Invalid B{C{radius}}.
203 '''
204 _, x, _ = self._a_x_b3(start, end, radius, wrap)
205 return _radians2m(x, radius)
207 def destination(self, distance, bearing, radius=R_M, height=None):
208 '''Locate the destination from this point after having
209 travelled the given distance on the given initial bearing.
211 @arg distance: Distance travelled (C{meter}, same units as
212 B{C{radius}}).
213 @arg bearing: Bearing from this point (compass C{degrees360}).
214 @kwarg radius: Mean earth radius (C{meter}).
215 @kwarg height: Optional height at destination (C{meter}, same
216 units a B{C{radius}}).
218 @return: Destination point (L{LatLon}).
220 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
221 B{C{radius}} or B{C{height}}.
222 '''
223 a, b = self.philam
224 r, t = _m2radians(distance, radius, low=None), Bearing_(bearing)
226 a, b = _destination2(a, b, r, t)
227 h = self._heigHt(height)
228 return self.classof(degrees90(a), degrees180(b), height=h)
230 def distanceTo(self, other, radius=R_M, wrap=False):
231 '''Compute the (angular) distance from this to an other point.
233 @arg other: The other point (L{LatLon}).
234 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
235 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
236 the B{C{other}} point (C{bool}).
238 @return: Distance between this and the B{C{other}} point
239 (C{meter}, same units as B{C{radius}} or
240 C{radians} if B{C{radius}} is C{None}).
242 @raise TypeError: The B{C{other}} point is not L{LatLon}.
244 @raise ValueError: Invalid B{C{radius}}.
245 '''
246 a1, _, a2, _, db = self._ab1_ab2_db5(other, wrap)
247 return _radians2m(vincentys_(a2, a1, db), radius)
249# @Property_RO
250# def Ecef(self):
251# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
252# '''
253# return _MODS.ecef.EcefKarney
255 def greatCircle(self, bearing, Vector=Vector3d, **Vector_kwds):
256 '''Compute the vector normal to great circle obtained by heading
257 on the given initial bearing from this point.
259 Direction of vector is such that initial bearing vector
260 b = c × n, where n is an n-vector representing this point.
262 @arg bearing: Bearing from this point (compass C{degrees360}).
263 @kwarg Vector: Vector class to return the great circle,
264 overriding the default L{Vector3d}.
265 @kwarg Vector_kwds: Optional, additional keyword argunents
266 for B{C{Vector}}.
268 @return: Vector representing great circle (C{Vector}).
270 @raise ValueError: Invalid B{C{bearing}}.
271 '''
272 a, b = self.philam
273 sa, ca, sb, cb, st, ct = sincos2_(a, b, Bearing_(bearing))
275 return Vector(sb * ct - cb * sa * st,
276 -cb * ct - sb * sa * st,
277 ca * st, **Vector_kwds) # XXX .unit()?
279 def initialBearingTo(self, other, wrap=False, raiser=False):
280 '''Compute the initial bearing (forward azimuth) from this
281 to an other point.
283 @arg other: The other point (spherical L{LatLon}).
284 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
285 the B{C{other}} point (C{bool}).
286 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}),
287 use C{B{raiser}=True} for behavior like
288 C{sphericalNvector.LatLon.initialBearingTo}.
290 @return: Initial bearing (compass C{degrees360}).
292 @raise CrossError: If this and the B{C{other}} point coincide,
293 provided both B{C{raiser}} is C{True} and
294 L{pygeodesy.crosserrors} is C{True}.
296 @raise TypeError: The B{C{other}} point is not L{LatLon}.
297 '''
298 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
299 # XXX behavior like sphericalNvector.LatLon.initialBearingTo
300 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(db)) < EPS:
301 raise CrossError(_point_, self, other=other, wrap=wrap, txt=_coincident_)
303 return degrees(bearing_(a1, b1, a2, b2, final=False))
305 def intermediateTo(self, other, fraction, height=None, wrap=False):
306 '''Locate the point at given fraction between (or along) this
307 and an other point.
309 @arg other: The other point (L{LatLon}).
310 @arg fraction: Fraction between both points (C{scalar},
311 0.0 at this and 1.0 at the other point).
312 @kwarg height: Optional height, overriding the intermediate
313 height (C{meter}).
314 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
315 B{C{other}} point (C{bool}).
317 @return: Intermediate point (L{LatLon}).
319 @raise TypeError: The B{C{other}} point is not L{LatLon}.
321 @raise ValueError: Invalid B{C{fraction}} or B{C{height}}.
323 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
324 '''
325 p = self
326 f = Scalar(fraction=fraction)
327 if not isnear0(f):
328 p = p.others(other)
329 if wrap:
330 p = _Wrap.point(p)
331 if not isnear1(f): # and not near0
332 a1, b1 = self.philam
333 a2, b2 = p.philam
334 db, b2 = unrollPI(b1, b2, wrap=wrap)
335 r = vincentys_(a2, a1, db)
336 sr = sin(r)
337 if isnon0(sr):
338 sa1, ca1, sa2, ca2, \
339 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2)
341 t = f * r
342 a = sin(r - t) # / sr superflous
343 b = sin( t) # / sr superflous
345 x = a * ca1 * cb1 + b * ca2 * cb2
346 y = a * ca1 * sb1 + b * ca2 * sb2
347 z = a * sa1 + b * sa2
349 a = atan1d(z, hypot(x, y))
350 b = atan2d(y, x)
352 else: # PYCHOK no cover
353 a = degrees90( favg(a1, a2, f=f)) # coincident
354 b = degrees180(favg(b1, b2, f=f))
356 h = self._havg(other, f=f, h=height)
357 p = self.classof(a, b, height=h)
358 return p
360 def intersection(self, end1, other, end2, height=None, wrap=False):
361 '''Compute the intersection point of two lines, each defined by
362 two points or a start point and bearing from North.
364 @arg end1: End point of this line (L{LatLon}) or the initial
365 bearing at this point (compass C{degrees360}).
366 @arg other: Start point of the other line (L{LatLon}).
367 @arg end2: End point of the other line (L{LatLon}) or the
368 initial bearing at the B{C{other}} point (compass
369 C{degrees360}).
370 @kwarg height: Optional height for intersection point,
371 overriding the mean height (C{meter}).
372 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
373 B{C{start2}} and both B{C{end*}} points (C{bool}).
375 @return: The intersection point (L{LatLon}). An alternate
376 intersection point might be the L{antipode} to
377 the returned result.
379 @raise IntersectionError: Ambiguous or infinite intersection
380 or colinear, parallel or otherwise
381 non-intersecting lines.
383 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}}
384 or B{C{end2}} not C{scalar} nor L{LatLon}.
386 @raise ValueError: Invalid B{C{height}} or C{null} line.
387 '''
388 try:
389 s2 = self.others(other)
390 return _intersect(self, end1, s2, end2, height=height, wrap=wrap,
391 LatLon=self.classof)
392 except (TypeError, ValueError) as x:
393 raise _xError(x, start1=self, end1=end1,
394 other=other, end2=end2, wrap=wrap)
396 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0,
397 height=None, wrap=True):
398 '''Compute the intersection points of two circles, each defined
399 by a center point and radius.
401 @arg rad1: Radius of the this circle (C{meter} or C{radians},
402 see B{C{radius}}).
403 @arg other: Center point of the other circle (L{LatLon}).
404 @arg rad2: Radius of the other circle (C{meter} or C{radians},
405 see B{C{radius}}).
406 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
407 B{C{rad2}} and B{C{eps}} are given in C{radians}).
408 @kwarg eps: Required overlap (C{meter} or C{radians}, see
409 B{C{radius}}).
410 @kwarg height: Optional height for the intersection points (C{meter},
411 conventionally) or C{None} for the I{"radical height"}
412 at the I{radical line} between both centers.
413 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
414 B{C{other}} point (C{bool}).
416 @return: 2-Tuple of the intersection points, each a L{LatLon}
417 instance. For abutting circles, both intersection
418 points are the same instance, aka the I{radical center}.
420 @raise IntersectionError: Concentric, antipodal, invalid or
421 non-intersecting circles.
423 @raise TypeError: If B{C{other}} is not L{LatLon}.
425 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
426 B{C{eps}} or B{C{height}}.
427 '''
428 try:
429 c2 = self.others(other)
430 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps,
431 height=height, wrap=wrap,
432 LatLon=self.classof)
433 except (TypeError, ValueError) as x:
434 raise _xError(x, center=self, rad1=rad1,
435 other=other, rad2=rad2, wrap=wrap)
437 @deprecated_method
438 def isEnclosedBy(self, points): # PYCHOK no cover
439 '''DEPRECATED, use method C{isenclosedBy}.'''
440 return self.isenclosedBy(points)
442 def isenclosedBy(self, points, wrap=False):
443 '''Check whether a (convex) polygon or composite encloses this point.
445 @arg points: The polygon points or composite (L{LatLon}[],
446 L{BooleanFHP} or L{BooleanGH}).
447 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
448 B{C{points}} (C{bool}).
450 @return: C{True} if this point is inside the polygon or
451 composite, C{False} otherwise.
453 @raise PointsError: Insufficient number of B{C{points}}.
455 @raise TypeError: Some B{C{points}} are not L{LatLon}.
457 @raise ValueError: Invalid B{C{points}}, non-convex polygon.
459 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
460 and L{pygeodesy.ispolar} especially if the B{C{points}} may
461 enclose a pole or wrap around the earth I{longitudinally}.
462 '''
463 if _MODS.booleans.isBoolean(points):
464 return points._encloses(self.lat, self.lon, wrap=wrap)
466 Ps = self.PointsIter(points, loop=2, dedup=True, wrap=wrap)
467 n0 = self._N_vector
469 v2 = Ps[0]._N_vector
470 p1 = Ps[1]
471 v1 = p1._N_vector
472 # check whether this point on same side of all
473 # polygon edges (to the left or right depending
474 # on the anti-/clockwise polygon direction)
475 gc1 = v2.cross(v1)
476 t0 = gc1.angleTo(n0) > PI_2
477 s0 = None
478 # get great-circle vector for each edge
479 for i, p2 in Ps.enumerate(closed=True):
480 if wrap and not Ps.looped:
481 p2 = _unrollon(p1, p2)
482 p1 = p2
483 v2 = p2._N_vector
484 gc = v1.cross(v2)
485 t = gc.angleTo(n0) > PI_2
486 if t != t0: # different sides of edge i
487 return False # outside
489 # check for convex polygon: angle between
490 # gc vectors, signed by direction of n0
491 # (otherwise the test above is not reliable)
492 s = signOf(gc1.angleTo(gc, vSign=n0))
493 if s != s0:
494 if s0 is None:
495 s0 = s
496 else:
497 t = _Fmt.SQUARE(points=i)
498 raise _ValueError(t, p2, wrap=wrap, txt=_not_(_convex_))
499 gc1, v1 = gc, v2
501 return True # inside
503 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
504 '''Find the midpoint between this and an other point.
506 @arg other: The other point (L{LatLon}).
507 @kwarg height: Optional height for midpoint, overriding
508 the mean height (C{meter}).
509 @kwarg fraction: Midpoint location from this point (C{scalar}),
510 may be negative or greater than 1.0.
511 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
512 B{C{other}} point (C{bool}).
514 @return: Midpoint (L{LatLon}).
516 @raise TypeError: The B{C{other}} point is not L{LatLon}.
518 @raise ValueError: Invalid B{C{height}}.
520 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
521 '''
522 if fraction is _0_5:
523 # see <https://MathForum.org/library/drmath/view/51822.html>
524 a1, b, a2, _, db = self._ab1_ab2_db5(other, wrap)
525 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db)
527 x = ca2 * cdb + ca1
528 y = ca2 * sdb
530 a = atan1d(sa1 + sa2, hypot(x, y))
531 b = degrees180(b + atan2(y, x))
533 h = self._havg(other, h=height)
534 r = self.classof(a, b, height=h)
535 else:
536 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
537 return r
539 def nearestOn(self, point1, point2, radius=R_M, **wrap_adjust_limit):
540 '''Locate the point between two points closest to this point.
542 Distances are approximated by function L{pygeodesy.equirectangular_},
543 subject to the supplied B{C{options}}.
545 @arg point1: Start point (L{LatLon}).
546 @arg point2: End point (L{LatLon}).
547 @kwarg radius: Mean earth radius (C{meter}).
548 @kwarg wrap_adjust_limit: Optional keyword arguments for functions
549 L{sphericalTrigonometry.nearestOn3} and
550 L{pygeodesy.equirectangular_},
552 @return: Closest point on the great circle line (L{LatLon}).
554 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
555 see function L{pygeodesy.equirectangular_}.
557 @raise NotImplementedError: Keyword argument C{B{within}=False}
558 is not (yet) supported.
560 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
562 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
564 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}
565 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
566 '''
567 # remove kwarg B{C{within}} if present
568 w = _xkwds_pop(wrap_adjust_limit, within=True)
569 if not w:
570 notImplemented(self, within=w)
572# # UNTESTED - handle C{B{within}=False} and C{B{within}=True}
573# wrap = _xkwds_get(options, wrap=False)
574# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap)
575# if fabs(a) < EPS or (within and a < EPS):
576# return point1
577# d = point1.distanceTo(point2, radius=radius, wrap=wrap)
578# if isnear0(d):
579# return point1 # or point2
580# elif fabs(d - a) < EPS or (a + EPS) > d:
581# return point2
582# f = a / d
583# if within:
584# if f > EPS1:
585# return point2
586# elif f < EPS:
587# return point1
588# return point1.intermediateTo(point2, f, wrap=wrap)
590 # without kwarg B{C{within}}, use backward compatible .nearestOn3
591 return self.nearestOn3([point1, point2], closed=False, radius=radius,
592 **wrap_adjust_limit)[0]
594 @deprecated_method
595 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover
596 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
598 @return: ... 2-Tuple C{(closest, distance)} of the closest
599 point (L{LatLon}) on the polygon and the distance
600 to that point from this point in C{meter}, same
601 units of B{C{radius}}.
602 '''
603 r = self.nearestOn3(points, closed=closed, radius=radius, **options)
604 return r.closest, r.distance
606 def nearestOn3(self, points, closed=False, radius=R_M, **wrap_adjust_limit):
607 '''Locate the point on a polygon closest to this point.
609 Distances are approximated by function L{pygeodesy.equirectangular_},
610 subject to the supplied B{C{options}}.
612 @arg points: The polygon points (L{LatLon}[]).
613 @kwarg closed: Optionally, close the polygon (C{bool}).
614 @kwarg radius: Mean earth radius (C{meter}).
615 @kwarg wrap_adjust_limit: Optional keyword arguments for function
616 L{sphericalTrigonometry.nearestOn3} and
617 L{pygeodesy.equirectangular_},
619 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
620 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_}
621 C{distance} between this and the C{closest} point converted to
622 C{meter}, same units as B{C{radius}}. The C{angle} from this
623 to the C{closest} point is in compass C{degrees360}, like
624 function L{pygeodesy.compassAngle}.
626 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
627 see function L{pygeodesy.equirectangular_}.
629 @raise PointsError: Insufficient number of B{C{points}}.
631 @raise TypeError: Some B{C{points}} are not C{LatLon}.
633 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
635 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_}
636 and L{pygeodesy.nearestOn5}.
637 '''
638 return nearestOn3(self, points, closed=closed, radius=radius,
639 LatLon=self.classof, **wrap_adjust_limit)
641 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
642 '''Convert this point to C{Karney}-based cartesian (ECEF)
643 coordinates.
645 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}}
646 and other keyword arguments, ignored
647 if C{B{Cartesian} is None}. Use
648 C{B{Cartesian}=...} to override
649 this L{Cartesian} class or specify
650 C{B{Cartesian}=None}.
652 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
653 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
654 C, M, datum)} with C{C} and C{M} if available.
656 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
657 '''
658 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
659 return LatLonSphericalBase.toCartesian(self, **kwds)
661 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
662 '''Compute the angles, sides and area of a spherical triangle.
664 @arg otherB: Second triangle point (C{LatLon}).
665 @arg otherC: Third triangle point (C{LatLon}).
666 @kwarg radius: Mean earth radius, ellipsoid or datum
667 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},
668 L{Datum} or L{a_f2Tuple}) or C{None}.
669 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
670 B{C{otherB}} and B{C{otherC}} points (C{bool}).
672 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if
673 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A,
674 a, B, b, C, c, D, E)}.
676 @see: Function L{triangle7} and U{Spherical trigonometry
677 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
678 '''
679 B = self.others(otherB=otherB)
680 C = self.others(otherC=otherC)
681 B, C, _ = _unrollon3(self, B, C, wrap)
683 r = self.philam + B.philam + C.philam
684 t = triangle8_(*r, wrap=wrap)
685 return self._xnamed(_t7Tuple(t, radius))
687 def triangulate(self, bearing1, other, bearing2, **height_wrap):
688 '''Locate a point given this, an other point and the (initial) bearing
689 at this and at the other point.
691 @arg bearing1: Bearing at this point (compass C{degrees360}).
692 @arg other: The other point (C{LatLon}).
693 @arg bearing2: Bearing at the other point (compass C{degrees360}).
694 @kwarg height_wrap_tol: Optional keyword arguments C{B{height}=None},
695 C{B{wrap}=False}, see method L{intersection}.
697 @return: Triangulated point (C{LatLon}).
699 @see: Method L{intersection} for further details.
700 '''
701 if _isDegrees(bearing1) and _isDegrees(bearing2):
702 return self.intersection(bearing1, other, bearing2, **height_wrap)
703 raise _TypeError(bearing1=bearing1, bearing2=bearing2, **height_wrap)
705 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
706 area=True, eps=EPS1, radius=R_M, wrap=False):
707 '''Trilaterate three points by I{area overlap} or I{perimeter
708 intersection} of three corresponding circles.
710 @arg distance1: Distance to this point (C{meter}, same units
711 as B{C{radius}}).
712 @arg point2: Second center point (C{LatLon}).
713 @arg distance2: Distance to point2 (C{meter}, same units as
714 B{C{radius}}).
715 @arg point3: Third center point (C{LatLon}).
716 @arg distance3: Distance to point3 (C{meter}, same units as
717 B{C{radius}}).
718 @kwarg area: If C{True} compute the area overlap, otherwise the
719 perimeter intersection of the circles (C{bool}).
720 @kwarg eps: The required I{minimal overlap} for C{B{area}=True}
721 or the I{intersection margin} for C{B{area}=False}
722 (C{meter}, same units as B{C{radius}}).
723 @kwarg radius: Mean earth radius (C{meter}, conventionally).
724 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
725 B{C{point2}} and B{C{point3}} (C{bool}).
727 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
728 with C{min} and C{max} in C{meter}, same units as B{C{eps}},
729 the corresponding trilaterated points C{minPoint} and
730 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number
731 of trilatered points found for the given B{C{eps}}.
733 If only a single trilaterated point is found, C{min I{is}
734 max}, C{minPoint I{is} maxPoint} and C{n = 1}.
736 For C{B{area}=True}, C{min} and C{max} are the smallest
737 respectively largest I{radial} overlap found.
739 For C{B{area}=False}, C{min} and C{max} represent the
740 nearest respectively farthest intersection margin.
742 If C{B{area}=True} and all 3 circles are concentric, C{n =
743 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}}
744 with the smallest B{C{distance#}} C{min} and C{max} the
745 largest B{C{distance#}}.
747 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
748 insufficient overlap for C{B{area}=True} or
749 no intersection or all (near-)concentric for
750 C{B{area}=False}.
752 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
754 @raise ValueError: Coincident B{C{point2}} or B{C{point3}} or invalid
755 B{C{distance1}}, B{C{distance2}}, B{C{distance3}}
756 or B{C{radius}}.
757 '''
758 return _trilaterate5(self, distance1,
759 self.others(point2=point2), distance2,
760 self.others(point3=point3), distance3,
761 area=area, radius=radius, eps=eps, wrap=wrap)
764_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
767def areaOf(points, radius=R_M, wrap=False): # was=True
768 '''Calculate the area of a (spherical) polygon or composite
769 (with the pointsjoined by great circle arcs).
771 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
772 or L{BooleanGH}).
773 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
774 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
775 or C{None}.
776 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{points}}
777 (C{bool}).
779 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
780 or C{radians} if B{C{radius}} is C{None}).
782 @raise PointsError: Insufficient number of B{C{points}}.
784 @raise TypeError: Some B{C{points}} are not L{LatLon}.
786 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
788 @note: The area is based on I{Karney}'s U{'Area of a spherical
789 polygon'<https://MathOverflow.net/questions/97711/
790 the-area-of-spherical-polygons>}, 3rd Answer.
792 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
793 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
794 L{ellipsoidalKarney.areaOf}.
795 '''
796 if _MODS.booleans.isBoolean(points):
797 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
799 _at2, _t_2 = atan2, tan_2
800 _un, _w180 = unrollPI, wrap180
802 Ps = _T00.PointsIter(points, loop=1, wrap=wrap)
803 p1 = p2 = Ps[0]
804 a1, b1 = p1.philam
805 ta1, z1 = _t_2(a1), None
807 A = Fsum() # mean phi
808 R = Fsum() # see L{pygeodesy.excessKarney_}
809 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
810 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
811 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
812 D = Fsum()
813 for i, p2 in Ps.enumerate(closed=True):
814 a2, b2 = p2.philam
815 db, b2 = _un(b1, b2, wrap=wrap and not Ps.looped)
816 A += a2
817 ta2 = _t_2(a2)
818 tdb = _t_2(db, points=i)
819 R += _at2(tdb * (ta1 + ta2),
820 _1_0 + ta1 * ta2)
821 ta1, b1 = ta2, b2
823 if not p2.isequalTo(p1, eps=EPS):
824 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
825 if z1 is not None:
826 D += _w180(z - z1) # (z - z1 + 540) ...
827 D += _w180(z2 - z) # (z2 - z + 540) % 360 - 180
828 p1, z1 = p2, z2
830 R = abs(R * _2_0)
831 if abs(D) < _90_0: # ispolar(points)
832 R = abs(R - PI2)
833 if radius:
834 a = degrees(A.fover(len(A))) # mean lat
835 R *= _mean_radius(radius, a)**2
836 return float(R)
839def _destination2(a, b, r, t):
840 '''(INTERNAL) Destination lat- and longitude in C{radians}.
842 @arg a: Latitude (C{radians}).
843 @arg b: Longitude (C{radians}).
844 @arg r: Angular distance (C{radians}).
845 @arg t: Bearing (compass C{radians}).
847 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
848 '''
849 # see <https://www.EdWilliams.org/avform.htm#LL>
850 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
851 ca *= sr
853 a = asin1(ct * ca + cr * sa)
854 d = atan2(st * ca, cr - sa * sin(a))
855 # note, in EdWilliams.org/avform.htm W is + and E is -
856 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
859def _int3d2(s, end, wrap, _i_, Vector, hs):
860 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
861 # and similar logic in .ellipsoidalBaseDI._intersect3
862 a1, b1 = s.philam
864 if _isDegrees(end): # bearing, get pseudo-end point
865 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
866 else: # must be a point
867 s.others(end, name=_end_ + _i_)
868 hs.append(end.height)
869 a2, b2 = end.philam
870 if wrap:
871 a2, b2 = _Wrap.philam(a2, b2)
873 db, b2 = unrollPI(b1, b2, wrap=wrap)
874 if max(fabs(db), fabs(a2 - a1)) < EPS:
875 raise _ValueError(_SPACE_(_line_ + _i_, _null_))
876 # note, in EdWilliams.org/avform.htm W is + and E is -
877 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
878 -(b1 + b2) * _0_5)
879 cb21 *= sin(a1 - a2) # sa21
880 sb21 *= sin(a1 + a2) # sa12
881 x = Vector(sb12 * cb21 - cb12 * sb21,
882 cb12 * cb21 + sb12 * sb21,
883 cos(a1) * cos(a2) * sin(db)) # ll=start
884 return x.unit(), (db, (a2 - a1)) # negated d
887def _intdot(ds, a1, b1, a, b, wrap):
888 # compute dot product ds . (-b + b1, a - a1)
889 db, _ = unrollPI(b1, b, wrap=wrap)
890 return fdot(ds, db, a - a1)
893def intersecant2(center, circle, point, other, **radius_exact_height_wrap):
894 '''Compute the intersections of a circle and a (great circle) line given as
895 two points or as a point and bearing.
897 @arg center: Center of the circle (L{LatLon}).
898 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
899 or a point on the circle (L{LatLon}).
900 @arg point: A point on the (great circle) line (L{LatLon}).
901 @arg other: An other point on the (great circle) line (L{LatLon}) or
902 the bearing at the B{C{point}} (compass C{degrees360}).
903 @kwarg radius_exact_height_wrap: Optional keyword arguments, see
904 method L{LatLon.intersecant2} for further details.
906 @return: 2-Tuple of the intersection points (representing a chord), each
907 an instance of the B{C{point}} class. Both points are the same
908 instance if the (great circle) line is tangent to the circle.
910 @raise IntersectionError: The circle and line do not intersect.
912 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
913 B{C{circle}} or B{C{other}} invalid.
915 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},
916 B{C{exact}}, B{C{height}} or B{C{napieradius}}.
917 '''
918 c = _T00.others(center=center)
919 p = _T00.others(point=point)
920 try:
921 return _intersecant2(c, circle, p, other, **radius_exact_height_wrap)
922 except (TypeError, ValueError) as x:
923 raise _xError(x, center=center, circle=circle, point=point, other=other,
924 **radius_exact_height_wrap)
927def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
928 LatLon=None, **LatLon_kwds):
929 # (INTERNAL) Intersect two (spherical) lines, see L{intersection}
930 # above, separated to allow callers to embellish any exceptions
932 s1, s2 = start1, start2
933 if wrap:
934 s2 = _Wrap.point(s2)
935 hs = [s1.height, s2.height]
937 a1, b1 = s1.philam
938 a2, b2 = s2.philam
939 db, b2 = unrollPI(b1, b2, wrap=wrap)
940 r12 = vincentys_(a2, a1, db)
941 if fabs(r12) < EPS: # [nearly] coincident points
942 a, b = favg(a1, a2), favg(b1, b2)
944 # see <https://www.EdWilliams.org/avform.htm#Intersection>
945 elif _isDegrees(end1) and _isDegrees(end2): # both bearings
946 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
948 x1, x2 = (sr12 * ca1), (sr12 * ca2)
949 if isnear0(x1) or isnear0(x2):
950 raise IntersectionError(_parallel_)
951 # handle domain error for equivalent longitudes,
952 # see also functions asin_safe and acos_safe at
953 # <https://www.EdWilliams.org/avform.htm#Math>
954 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1)
955 t21, t23 = acos1((sa1 - sa2 * cr12) / x2), radiansPI2(end2)
956 if sin(db) > 0:
957 t21 = PI2 - t21
958 else:
959 t12 = PI2 - t12
960 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
961 wrapPI(t21 - t23)) # angle 1-2-3)
962 if isnear0(sx1) and isnear0(sx2):
963 raise IntersectionError(_infinite_)
964 sx3 = sx1 * sx2
965# XXX if sx3 < 0:
966# XXX raise ValueError(_ambiguous_)
967 x3 = acos1(cr12 * sx3 - cx2 * cx1)
968 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
970 a, b = _destination2(a1, b1, r13, t13)
971 # like .ellipsoidalBaseDI,_intersect3, if this intersection
972 # is "before" the first point, use the antipodal intersection
973 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
974 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
976 else: # end point(s) or bearing(s)
977 _N_vector_ = _MODS.nvectorBase._N_vector_
979 x1, d1 = _int3d2(s1, end1, wrap, _1_, _N_vector_, hs)
980 x2, d2 = _int3d2(s2, end2, wrap, _2_, _N_vector_, hs)
981 x = x1.cross(x2)
982 if x.length < EPS: # [nearly] colinear or parallel lines
983 raise IntersectionError(_colinear_)
984 a, b = x.philam
985 # choose intersection similar to sphericalNvector
986 if not (_intdot(d1, a1, b1, a, b, wrap) *
987 _intdot(d2, a2, b2, a, b, wrap)) > 0:
988 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
990 h = fmean(hs) if height is None else Height(height)
991 return _LL3Tuple(degrees90(a), degrees180(b), h,
992 intersection, LatLon, LatLon_kwds)
995def intersection(start1, end1, start2, end2, height=None, wrap=False,
996 LatLon=LatLon, **LatLon_kwds):
997 '''Compute the intersection point of two lines, each defined
998 by two points or a start point and bearing from North.
1000 @arg start1: Start point of the first line (L{LatLon}).
1001 @arg end1: End point of the first line (L{LatLon}) or
1002 the initial bearing at the first start point
1003 (compass C{degrees360}).
1004 @arg start2: Start point of the second line (L{LatLon}).
1005 @arg end2: End point of the second line (L{LatLon}) or
1006 the initial bearing at the second start point
1007 (compass C{degrees360}).
1008 @kwarg height: Optional height for the intersection point,
1009 overriding the mean height (C{meter}).
1010 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1011 B{C{start2}} and both B{C{end*}} points (C{bool}).
1012 @kwarg LatLon: Optional class to return the intersection
1013 point (L{LatLon}) or C{None}.
1014 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1015 arguments, ignored if C{B{LatLon} is None}.
1017 @return: The intersection point as a (B{C{LatLon}}) or if
1018 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon,
1019 height)}. An alternate intersection point might
1020 be the L{antipode} to the returned result.
1022 @raise IntersectionError: Ambiguous or infinite intersection
1023 or colinear, parallel or otherwise
1024 non-intersecting lines.
1026 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}}
1027 or B{C{end2}} point not L{LatLon}.
1029 @raise ValueError: Invalid B{C{height}} or C{null} line.
1030 '''
1031 s1 = _T00.others(start1=start1)
1032 s2 = _T00.others(start2=start2)
1033 try:
1034 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap,
1035 LatLon=LatLon, **LatLon_kwds)
1036 except (TypeError, ValueError) as x:
1037 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1040def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1041 height=None, wrap=False, # was=True
1042 LatLon=LatLon, **LatLon_kwds):
1043 '''Compute the intersection points of two circles each defined
1044 by a center point and a radius.
1046 @arg center1: Center of the first circle (L{LatLon}).
1047 @arg rad1: Radius of the first circle (C{meter} or C{radians},
1048 see B{C{radius}}).
1049 @arg center2: Center of the second circle (L{LatLon}).
1050 @arg rad2: Radius of the second circle (C{meter} or C{radians},
1051 see B{C{radius}}).
1052 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1053 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1054 @kwarg eps: Required overlap (C{meter} or C{radians}, see
1055 B{C{radius}}).
1056 @kwarg height: Optional height for the intersection points (C{meter},
1057 conventionally) or C{None} for the I{"radical height"}
1058 at the I{radical line} between both centers.
1059 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{center2}}
1060 (C{bool}).
1061 @kwarg LatLon: Optional class to return the intersection
1062 points (L{LatLon}) or C{None}.
1063 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1064 arguments, ignored if C{B{LatLon} is None}.
1066 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1067 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1068 lon, height)}. For abutting circles, both intersection
1069 points are the same instance, aka the I{radical center}.
1071 @raise IntersectionError: Concentric, antipodal, invalid or
1072 non-intersecting circles.
1074 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1076 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1077 B{C{eps}} or B{C{height}}.
1079 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1081 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1082 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1083 '''
1084 c1 = _T00.others(center1=center1)
1085 c2 = _T00.others(center2=center2)
1086 try:
1087 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1088 height=height, wrap=wrap,
1089 LatLon=LatLon, **LatLon_kwds)
1090 except (TypeError, ValueError) as x:
1091 raise _xError(x, center1=center1, rad1=rad1,
1092 center2=center2, rad2=rad2, wrap=wrap)
1095def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1096 height=None, too_d=None, wrap=False, # was=True
1097 LatLon=LatLon, **LatLon_kwds):
1098 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1099 # above, separated to allow callers to embellish any exceptions
1101 def _dest3(bearing, h):
1102 a, b = _destination2(a1, b1, r1, bearing)
1103 return _LL3Tuple(degrees90(a), degrees180(b), h,
1104 intersections2, LatLon, LatLon_kwds)
1106 a1, b1 = c1.philam
1107 a2, b2 = c2.philam
1108 if wrap:
1109 a2, b2 = _Wrap.philam(a2, b2)
1111 r1, r2, f = _rads3(rad1, rad2, radius)
1112 if f: # swapped radii, swap centers
1113 a1, a2 = a2, a1 # PYCHOK swap!
1114 b1, b2 = b2, b1 # PYCHOK swap!
1116 db, b2 = unrollPI(b1, b2, wrap=wrap)
1117 d = vincentys_(a2, a1, db) # radians
1118 if d < max(r1 - r2, EPS):
1119 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1121 r = eps if radius is None else (m2radians(
1122 eps, radius=radius) if eps else _0_0)
1123 if r < _0_0:
1124 raise _ValueError(eps=r)
1126 x = fsumf_(r1, r2, -d) # overlap
1127 if x > max(r, EPS):
1128 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1129 x = sd * sr1
1130 if isnear0(x):
1131 raise _ValueError(_invalid_)
1132 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1134 elif x < r: # PYCHOK no cover
1135 t = (d * radius) if too_d is None else too_d
1136 raise IntersectionError(_too_(_Fmt.distant(t)))
1138 if height is None: # "radical height"
1139 f = _radical2(d, r1, r2).ratio
1140 h = Height(favg(c1.height, c2.height, f=f))
1141 else:
1142 h = Height(height)
1144 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1145 if x < EPS4: # externally ...
1146 r = _dest3(b, h)
1147 elif x > _PI_EPS4: # internally ...
1148 r = _dest3(b + PI, h)
1149 else:
1150 return _dest3(b + x, h), _dest3(b - x, h)
1151 return r, r # ... abutting circles
1154@deprecated_function
1155def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1156 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1157 '''
1158 return ispolar(points, wrap=wrap)
1161def _LL3Tuple(lat, lon, height, where, LatLon, LatLon_kwds):
1162 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1163 '''
1164 n = where.__name__
1165 if LatLon is None:
1166 r = LatLon3Tuple(lat, lon, height, name=n)
1167 else:
1168 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1169 r = LatLon(lat, lon, **kwds)
1170 return r
1173def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):
1174 '''Compute the I{geographic} mean of several points.
1176 @arg points: Points to be averaged (L{LatLon}[]).
1177 @kwarg height: Optional height at mean point, overriding the mean
1178 height (C{meter}).
1179 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}}
1180 (C{bool}).
1181 @kwarg LatLon: Optional class to return the mean point (L{LatLon})
1182 or C{None}.
1183 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1184 arguments, ignored if C{B{LatLon} is None}.
1186 @return: The geographic mean and height (B{C{LatLon}}) or a
1187 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}}
1188 is C{None}.
1190 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1192 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1193 '''
1194 def _N_vs(ps, w):
1195 Ps = _T00.PointsIter(ps, wrap=w)
1196 for p in Ps.iterate(closed=False):
1197 yield p._N_vector
1199 m = _MODS.nvectorBase
1200 # geographic, vectorial mean
1201 n = m.sumOf(_N_vs(points, wrap), h=height, Vector=m.NvectorBase)
1202 lat, lon, h = n.latlonheight
1203 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1206@deprecated_function
1207def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1208 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1210 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1211 point (L{LatLon}) on the polygon and the C{distance}
1212 between the C{closest} and the given B{C{point}}. The
1213 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1214 lon)} if B{C{LatLon}} is C{None} ...
1215 '''
1216 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1217 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1218 ll = LatLon2Tuple(ll.lat, ll.lon)
1219 return ll, d
1222def nearestOn3(point, points, closed=False, radius=R_M, wrap=False, adjust=True,
1223 limit=9, **LatLon_and_kwds):
1224 '''Locate the point on a path or polygon closest to a reference point.
1226 Distances are I{approximated} using function L{pygeodesy.equirectangular_},
1227 subject to the supplied B{C{options}}.
1229 @arg point: The reference point (L{LatLon}).
1230 @arg points: The path or polygon points (L{LatLon}[]).
1231 @kwarg closed: Optionally, close the polygon (C{bool}).
1232 @kwarg radius: Mean earth radius (C{meter}).
1233 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1234 B{C{points}} (C{bool}).
1235 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}).
1236 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}),
1237 default C{9 degrees} is about C{1,000 Kmeter} (for mean
1238 spherical earth radius L{R_KM}).
1239 @kwarg LatLon: Optional class to return the closest point (L{LatLon})
1240 or C{None}.
1241 @kwarg options: Optional keyword arguments for function
1242 L{pygeodesy.equirectangular_}.
1244 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1245 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1246 lon, height)} if B{C{LatLon}} is C{None}. The C{distance}
1247 is the L{pygeodesy.equirectangular_} distance between the
1248 C{closest} and the given B{C{point}} converted to C{meter},
1249 same units as B{C{radius}}. The C{angle} from the given
1250 B{C{point}} to the C{closest} is in compass C{degrees360},
1251 like function L{pygeodesy.compassAngle}. The C{height} is
1252 the (interpolated) height at the C{closest} point.
1254 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1255 see function L{pygeodesy.equirectangular_}.
1257 @raise PointsError: Insufficient number of B{C{points}}.
1259 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1261 @raise ValueError: Invalid B{C{radius}}.
1263 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}.
1264 '''
1265 t = _nearestOn5(point, points, closed=closed, wrap=wrap,
1266 adjust=adjust, limit=limit)
1267 d = degrees2m(t.distance, radius=radius)
1268 h = t.height
1269 n = nearestOn3.__name__
1271 kwds = _xkwds(LatLon_and_kwds, height=h, name=n)
1272 LL = _xkwds_pop(kwds, LatLon=LatLon)
1273 r = LatLon3Tuple(t.lat, t.lon, h, name=n) if LL is None else \
1274 LL(t.lat, t.lon, **kwds)
1275 return NearestOn3Tuple(r, d, t.angle, name=n)
1278def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1279 '''Compute the perimeter of a (spherical) polygon or composite
1280 (with great circle arcs joining the points).
1282 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1283 or L{BooleanGH}).
1284 @kwarg closed: Optionally, close the polygon (C{bool}).
1285 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1286 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1287 B{C{points}} (C{bool}).
1289 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1290 or C{radians} if B{C{radius}} is C{None}).
1292 @raise PointsError: Insufficient number of B{C{points}}.
1294 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1296 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1297 C{B{points}} a composite.
1299 @note: Distances are based on function L{pygeodesy.vincentys_}.
1301 @see: Functions L{perimeterOf<pygeodesy.perimeterOf>},
1302 L{sphericalNvector.perimeterOf} and L{ellipsoidalKarney.perimeterOf}.
1303 '''
1304 def _rads(ps, c, w): # angular edge lengths in radians
1305 Ps = _T00.PointsIter(ps, loop=1, wrap=w)
1306 a1, b1 = Ps[0].philam
1307 for p in Ps.iterate(closed=c):
1308 a2, b2 = p.philam
1309 db, b2 = unrollPI(b1, b2, wrap=w and not (c and Ps.looped))
1310 yield vincentys_(a2, a1, db)
1311 a1, b1 = a2, b2
1313 if _MODS.booleans.isBoolean(points):
1314 if not closed:
1315 raise _ValueError(closed=closed, points=_composite_)
1316 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1317 else:
1318 r = fsum(_rads(points, closed, wrap), floats=True)
1319 return _radians2m(r, radius)
1322def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1323 excess=excessAbc_,
1324 wrap=False):
1325 '''Compute the angles, sides, and area of a (spherical) triangle.
1327 @arg latA: First corner latitude (C{degrees}).
1328 @arg lonA: First corner longitude (C{degrees}).
1329 @arg latB: Second corner latitude (C{degrees}).
1330 @arg lonB: Second corner longitude (C{degrees}).
1331 @arg latC: Third corner latitude (C{degrees}).
1332 @arg lonC: Third corner longitude (C{degrees}).
1333 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1334 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1335 or C{None}.
1336 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1337 L{excessGirard_} or L{excessLHuilier_}).
1338 @kwarg wrap: If C{True}, wrap and L{pygeodesy.unroll180}
1339 longitudes (C{bool}).
1341 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1342 spherical angles C{A}, C{B} and C{C}, angular sides
1343 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1344 in I{square} C{meter} or same units as B{C{radius}}
1345 I{squared} or if C{B{radius}=0} or C{None}, a
1346 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in
1347 C{radians} with the I{spherical excess} C{E} as the
1348 C{unit area} in C{radians}.
1349 '''
1350 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA),
1351 Phi_(latB=latB), Lam_(lonB=lonB),
1352 Phi_(latC=latC), Lam_(lonC=lonC),
1353 excess=excess, wrap=wrap)
1354 return _t7Tuple(t, radius)
1357def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_,
1358 wrap=False):
1359 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1360 excess} of a (spherical) triangle.
1362 @arg phiA: First corner latitude (C{radians}).
1363 @arg lamA: First corner longitude (C{radians}).
1364 @arg phiB: Second corner latitude (C{radians}).
1365 @arg lamB: Second corner longitude (C{radians}).
1366 @arg phiC: Third corner latitude (C{radians}).
1367 @arg lamC: Third corner longitude (C{radians}).
1368 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1369 L{excessGirard_} or L{excessLHuilier_}).
1370 @kwarg wrap: If C{True}, L{pygeodesy.unrollPI} the
1371 longitudinal deltas (C{bool}).
1373 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1374 spherical angles C{A}, C{B} and C{C}, angular sides
1375 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1376 I{spherical excess} C{E}, all in C{radians}.
1377 '''
1378 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1379 d, _ = unrollPI(lamB, lamC, wrap=w)
1380 a = vincentys_(phiC, phiB, d)
1381 return a, (phiB, lamB, phiC, lamC, phiA, lamA) # rotate A, B, C
1383 def _A_r(a, sa, ca, sb, cb, sc, cc):
1384 s = sb * sc
1385 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1386 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2_'s
1388 # notation: side C{a} is oposite to corner C{A}, etc.
1389 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1390 b, r = _a_r(wrap, *r)
1391 c, _ = _a_r(wrap, *r)
1393 A, r = _A_r(a, *sincos2_(a, b, c))
1394 B, r = _A_r(b, *r)
1395 C, _ = _A_r(c, *r)
1397 D = fsumf_(PI2, -a, -b, -c) # deficit aka defect
1398 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else (
1399 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else
1400 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b))))
1402 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1405def _t7Tuple(t, radius):
1406 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1407 '''
1408 if radius: # not in (None, _0_0)
1409 r = radius if _isRadius(radius) else \
1410 _ellipsoidal_datum(radius).ellipsoid.Rmean
1411 A, B, C = map1(degrees, t.A, t.B, t.C)
1412 t = Triangle7Tuple(A, (r * t.a),
1413 B, (r * t.b),
1414 C, (r * t.c), t.E * r**2)
1415 return t
1418__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1419 areaOf, # functions
1420 intersecant2, intersection, intersections2, ispolar,
1421 isPoleEnclosedBy, # DEPRECATED, use ispolar
1422 meanOf,
1423 nearestOn2, nearestOn3,
1424 perimeterOf,
1425 sumOf, # XXX == vector3d.sumOf
1426 triangle7, triangle8_)
1428# **) MIT License
1429#
1430# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1431#
1432# Permission is hereby granted, free of charge, to any person obtaining a
1433# copy of this software and associated documentation files (the "Software"),
1434# to deal in the Software without restriction, including without limitation
1435# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1436# and/or sell copies of the Software, and to permit persons to whom the
1437# Software is furnished to do so, subject to the following conditions:
1438#
1439# The above copyright notice and this permission notice shall be included
1440# in all copies or substantial portions of the Software.
1441#
1442# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1443# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1444# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1445# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1446# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1447# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1448# OTHER DEALINGS IN THE SOFTWARE.