Coverage for pygeodesy/sphericalTrigonometry.py: 93%
387 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-11-12 16:17 -0500
« prev ^ index » next coverage.py v7.6.1, created at 2024-11-12 16:17 -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-2024} 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_pop2
27from pygeodesy.fmath import favg, fdot, 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_, \
34 _invalid_, _line_, _near_, _null_, _parallel_, \
35 _point_, _SPACE_, _too_
36from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER
37# from pygeodesy.nvectorBase import NvectorBase, sumOf # _MODE
38from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, NearestOn3Tuple, \
39 Triangle7Tuple, Triangle8Tuple
40from pygeodesy.points import ispolar, nearestOn5 as _nearestOn5, \
41 Fmt as _Fmt # XXX shadowed
42from pygeodesy.props import deprecated_function, deprecated_method
43from pygeodesy.sphericalBase import _m2radians, CartesianSphericalBase, \
44 _intersecant2, LatLonSphericalBase, \
45 _rads3, _radians2m, _trilaterate5
46# from pygeodesy.streprs import Fmt as _Fmt # from .points XXX shadowed
47from pygeodesy.units import Bearing_, Height, _isDegrees, _isRadius, Lamd, \
48 Phid, Radius_, Scalar
49from pygeodesy.utily import acos1, asin1, atan1d, atan2d, degrees90, degrees180, \
50 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \
51 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI
52from pygeodesy.vector3d import sumOf, Vector3d
54from math import asin, atan2, cos, degrees, fabs, radians, sin
56__all__ = _ALL_LAZY.sphericalTrigonometry
57__version__ = '24.11.07'
59_PI_EPS4 = PI - EPS4
60if _PI_EPS4 >= PI:
61 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4)
64class Cartesian(CartesianSphericalBase):
65 '''Extended to convert geocentric, L{Cartesian} points to
66 spherical, geodetic L{LatLon}.
67 '''
69 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon
70 '''Convert this cartesian point to a C{spherical} geodetic point.
72 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
73 arguments. Use C{B{LatLon}=...} to override
74 this L{LatLon} class or specify C{B{LatLon}=None}.
76 @return: The geodetic point (L{LatLon}) or if C{B{LatLon} is None},
77 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
78 with C{C} and C{M} if available.
80 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
81 '''
82 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
83 return CartesianSphericalBase.toLatLon(self, **kwds)
86class LatLon(LatLonSphericalBase):
87 '''New point on a spherical earth model, based on trigonometry formulae.
88 '''
90 def _ab1_ab2_db5(self, other, wrap):
91 '''(INTERNAL) Helper for several methods.
92 '''
93 a1, b1 = self.philam
94 a2, b2 = self.others(other, up=2).philam
95 if wrap:
96 a2, b2 = _Wrap.philam(a2, b2)
97 db, b2 = unrollPI(b1, b2, wrap=wrap)
98 else: # unrollPI shortcut
99 db = b2 - b1
100 return a1, b1, a2, b2, db
102 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
103 '''Compute the (signed) distance from the start to the closest
104 point on the great circle line defined by a start and an
105 end point.
107 That is, if a perpendicular is drawn from this point to the
108 great circle line, the along-track distance is the distance
109 from the start point to the point where the perpendicular
110 crosses the line.
112 @arg start: Start point of the great circle line (L{LatLon}).
113 @arg end: End point of the great circle line (L{LatLon}).
114 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
115 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
116 the B{C{start}} and B{C{end}} point (C{bool}).
118 @return: Distance along the great circle line (C{radians}
119 if C{B{radius} is None} or C{meter}, same units
120 as B{C{radius}}), positive if I{after} the
121 B{C{start}} toward the B{C{end}} point of the
122 line, I{negative} if before or C{0} if at the
123 B{C{start}} point.
125 @raise TypeError: Invalid B{C{start}} or B{C{end}} point.
127 @raise ValueError: Invalid B{C{radius}}.
128 '''
129 r, x, b = self._a_x_b3(start, end, radius, wrap)
130 cx = cos(x)
131 return _0_0 if isnear0(cx) else \
132 _radians2m(copysign0(acos1(cos(r) / cx), cos(b)), radius)
134 def _a_x_b3(self, start, end, radius, wrap):
135 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo.
136 '''
137 s = self.others(start=start)
138 e = self.others(end=end)
139 s, e, w = _unrollon3(self, s, e, wrap)
141 r = Radius_(radius)
142 r = s.distanceTo(self, r, wrap=w) / r
144 b = radians(s.initialBearingTo(self, wrap=w)
145 - s.initialBearingTo(e, wrap=w))
146 x = asin(sin(r) * sin(b))
147 return r, x, -b
149 @deprecated_method
150 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover
151 '''DEPRECATED, use method L{initialBearingTo}.
152 '''
153 return self.initialBearingTo(other, wrap=wrap, raiser=raiser)
155 def crossingParallels(self, other, lat, wrap=False):
156 '''Return the pair of meridians at which a great circle defined
157 by this and an other point crosses the given latitude.
159 @arg other: The other point defining great circle (L{LatLon}).
160 @arg lat: Latitude at the crossing (C{degrees}).
161 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
162 B{C{other}} point (C{bool}).
164 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
165 C{None} if the great circle doesn't reach B{C{lat}}.
166 '''
167 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
168 sa, ca, sa1, ca1, \
169 sa2, ca2, sdb, cdb = sincos2_(radians(lat), a1, a2, db)
170 sa1 *= ca2 * ca
172 x = sa1 * sdb
173 y = sa1 * cdb - ca1 * sa2 * ca
174 z = ca1 * sdb * ca2 * sa
176 h = hypot(x, y)
177 if h < EPS or fabs(z) > h: # PYCHOK no cover
178 return None # great circle doesn't reach latitude
180 m = atan2(-y, x) + b1 # longitude at max latitude
181 d = acos1(z / h) # delta longitude to intersections
182 return degrees180(m - d), degrees180(m + d)
184 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
185 '''Compute the (signed) distance from this point to a great
186 circle from a start to an end point.
188 @arg start: Start point of the great circle line (L{LatLon}).
189 @arg end: End point of the great circle line (L{LatLon}).
190 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
191 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
192 the B{C{start}} and B{C{end}} point (C{bool}).
194 @return: Distance to the great circle (C{radians} if
195 B{C{radius}} or C{meter}, same units as
196 B{C{radius}}), I{negative} if to the left or
197 I{positive} if to the right of the line.
199 @raise TypeError: If B{C{start}} or B{C{end}} is not L{LatLon}.
201 @raise ValueError: Invalid B{C{radius}}.
202 '''
203 _, x, _ = self._a_x_b3(start, end, radius, wrap)
204 return _radians2m(x, radius)
206 def destination(self, distance, bearing, radius=R_M, height=None):
207 '''Locate the destination from this point after having
208 travelled the given distance on a bearing from North.
210 @arg distance: Distance travelled (C{meter}, same units as
211 B{C{radius}}).
212 @arg bearing: Bearing from this point (compass C{degrees360}).
213 @kwarg radius: Mean earth radius (C{meter}).
214 @kwarg height: Optional height at destination (C{meter}, same
215 units a B{C{radius}}).
217 @return: Destination point (L{LatLon}).
219 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
220 B{C{radius}} or B{C{height}}.
221 '''
222 a, b = self.philam
223 r, t = _m2radians(distance, radius, low=None), Bearing_(bearing)
225 a, b = _destination2(a, b, r, t)
226 h = self._heigHt(height)
227 return self.classof(degrees90(a), degrees180(b), height=h)
229 def distanceTo(self, other, radius=R_M, wrap=False):
230 '''Compute the (angular) distance from this to an other point.
232 @arg other: The other point (L{LatLon}).
233 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
234 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
235 the B{C{other}} point (C{bool}).
237 @return: Distance between this and the B{C{other}} point
238 (C{meter}, same units as B{C{radius}} or
239 C{radians} if C{B{radius} is None}).
241 @raise TypeError: The B{C{other}} point is not L{LatLon}.
243 @raise ValueError: Invalid B{C{radius}}.
244 '''
245 a1, _, a2, _, db = self._ab1_ab2_db5(other, wrap)
246 return _radians2m(vincentys_(a2, a1, db), radius)
248# @Property_RO
249# def Ecef(self):
250# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
251# '''
252# return _MODS.ecef.EcefKarney
254 def greatCircle(self, bearing, Vector=Vector3d, **Vector_kwds):
255 '''Compute the vector normal to great circle obtained by heading
256 from this point on the bearing from North.
258 Direction of vector is such that initial bearing vector
259 b = c × n, where n is an n-vector representing this point.
261 @arg bearing: Bearing from this point (compass C{degrees360}).
262 @kwarg Vector: Vector class to return the great circle,
263 overriding the default L{Vector3d}.
264 @kwarg Vector_kwds: Optional, additional keyword argunents
265 for B{C{Vector}}.
267 @return: Vector representing great circle (C{Vector}).
269 @raise ValueError: Invalid B{C{bearing}}.
270 '''
271 a, b = self.philam
272 sa, ca, sb, cb, st, ct = sincos2_(a, b, Bearing_(bearing))
274 sa *= st
275 return Vector(fdot_(sb, ct, -cb, sa),
276 -fdot_(cb, ct, sb, sa),
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 and if B{C{raiser}} and L{crosserrors
294 <pygeodesy.crosserrors>} are both 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 = fdot_(a, ca1 * cb1, b, ca2 * cb2)
346 y = fdot_(a, ca1 * sb1, b, ca2 * sb2)
347 z = fdot_(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 a 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 a 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 other points closest to this point.
542 Distances are approximated by function L{pygeodesy.equirectangular4},
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.equirectangular4},
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.equirectangular4}.
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.equirectangular4} and L{pygeodesy.nearestOn5}
565 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
566 '''
567 # remove kwarg B{C{within}} if present
568 w, kwds = _xkwds_pop2(wrap_adjust_limit, within=True)
569 if not w:
570 self._notImplemented(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 **kwds)[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.equirectangular4},
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.equirectangular4},
619 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
620 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular4}
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.equirectangular4}.
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.equirectangular4}
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) coordinates.
644 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}} and other
645 keyword arguments, ignored if C{B{Cartesian} is
646 None}. Use C{B{Cartesian}=...} to override this
647 L{Cartesian} class or specify C{B{Cartesian}=None}.
649 @return: The cartesian point (L{Cartesian}) or if C{B{Cartesian} is None},
650 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with C{C}
651 and C{M} if available.
653 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
654 '''
655 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
656 return LatLonSphericalBase.toCartesian(self, **kwds)
658 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
659 '''Compute the angles, sides and area of a spherical triangle.
661 @arg otherB: Second triangle point (C{LatLon}).
662 @arg otherC: Third triangle point (C{LatLon}).
663 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, L{Ellipsoid},
664 L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) or C{None}.
665 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll points B{C{otherB}}
666 and B{C{otherC}} (C{bool}).
668 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if B{C{radius} is
669 None}, a L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)}.
671 @see: Function L{triangle7} and U{Spherical trigonometry
672 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
673 '''
674 B = self.others(otherB=otherB)
675 C = self.others(otherC=otherC)
676 B, C, _ = _unrollon3(self, B, C, wrap)
678 r = self.philam + B.philam + C.philam
679 t = triangle8_(*r, wrap=wrap)
680 return self._xnamed(_t7Tuple(t, radius))
682 def triangulate(self, bearing1, other, bearing2, **height_wrap):
683 '''Locate a point given this, an other point and a bearing from
684 North at both points.
686 @arg bearing1: Bearing at this point (compass C{degrees360}).
687 @arg other: The other point (C{LatLon}).
688 @arg bearing2: Bearing at the other point (compass C{degrees360}).
689 @kwarg height_wrap_tol: Optional keyword arguments C{B{height}=None},
690 C{B{wrap}=False}, see method L{intersection}.
692 @return: Triangulated point (C{LatLon}).
694 @see: Method L{intersection} for further details.
695 '''
696 if _isDegrees(bearing1) and _isDegrees(bearing2):
697 return self.intersection(bearing1, other, bearing2, **height_wrap)
698 raise _TypeError(bearing1=bearing1, bearing2=bearing2, **height_wrap)
700 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
701 area=True, eps=EPS1, radius=R_M, wrap=False):
702 '''Trilaterate three points by I{area overlap} or I{perimeter intersection}
703 of three corresponding circles.
705 @arg distance1: Distance to this point (C{meter}, same units as B{C{radius}}).
706 @arg point2: Second center point (C{LatLon}).
707 @arg distance2: Distance to point2 (C{meter}, same units as B{C{radius}}).
708 @arg point3: Third center point (C{LatLon}).
709 @arg distance3: Distance to point3 (C{meter}, same units as B{C{radius}}).
710 @kwarg area: If C{True}, compute the area overlap, otherwise the perimeter
711 intersection of the circles (C{bool}).
712 @kwarg eps: The required I{minimal overlap} for C{B{area}=True} or the
713 I{intersection margin} if C{B{area}=False} (C{meter}, same
714 units as B{C{radius}}).
715 @kwarg radius: Mean earth radius (C{meter}, conventionally).
716 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{point2}} and
717 B{C{point3}} (C{bool}).
719 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)} with
720 C{min} and C{max} in C{meter}, same units as B{C{eps}}, the
721 corresponding trilaterated points C{minPoint} and C{maxPoint}
722 as I{spherical} C{LatLon} and C{n}, the number of trilatered
723 points found for the given B{C{eps}}.
725 If only a single trilaterated point is found, C{min I{is} max},
726 C{minPoint I{is} maxPoint} and C{n = 1}.
728 For C{B{area}=True}, C{min} and C{max} are the smallest respectively
729 largest I{radial} overlap found.
731 For C{B{area}=False}, C{min} and C{max} represent the nearest
732 respectively farthest intersection margin.
734 If C{B{area}=True} and all 3 circles are concentric, C{n=0} and
735 C{minPoint} and C{maxPoint} are both the B{C{point#}} with the
736 smallest B{C{distance#}} C{min} and C{max} the largest B{C{distance#}}.
738 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
739 insufficient overlap for C{B{area}=True} or
740 no intersection or all (near-)concentric if
741 C{B{area}=False}.
743 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
745 @raise ValueError: Coincident B{C{point2}} or B{C{point3}} or invalid
746 B{C{distance1}}, B{C{distance2}}, B{C{distance3}}
747 or B{C{radius}}.
748 '''
749 return _trilaterate5(self, distance1,
750 self.others(point2=point2), distance2,
751 self.others(point3=point3), distance3,
752 area=area, radius=radius, eps=eps, wrap=wrap)
755_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
758def areaOf(points, radius=R_M, wrap=False): # was=True
759 '''Calculate the area of a (spherical) polygon or composite (with the
760 points joined by great circle arcs).
762 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
763 or L{BooleanGH}).
764 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
765 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
766 or C{None}.
767 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{points}}
768 (C{bool}).
770 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
771 or C{radians} if C{B{radius} is None}).
773 @raise PointsError: Insufficient number of B{C{points}}.
775 @raise TypeError: Some B{C{points}} are not L{LatLon}.
777 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
779 @note: The area is based on I{Karney}'s U{'Area of a spherical
780 polygon'<https://MathOverflow.net/questions/97711/
781 the-area-of-spherical-polygons>}, 3rd Answer.
783 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
784 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
785 L{ellipsoidalKarney.areaOf}.
786 '''
787 if _MODS.booleans.isBoolean(points):
788 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
790 _at2, _t_2 = atan2, tan_2
791 _un, _w180 = unrollPI, wrap180
793 Ps = _T00.PointsIter(points, loop=1, wrap=wrap)
794 p1 = p2 = Ps[0]
795 a1, b1 = p1.philam
796 ta1, z1 = _t_2(a1), None
798 A = Fsum() # mean phi
799 R = Fsum() # see L{pygeodesy.excessKarney_}
800 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
801 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
802 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
803 D = Fsum()
804 for i, p2 in Ps.enumerate(closed=True):
805 a2, b2 = p2.philam
806 db, b2 = _un(b1, b2, wrap=wrap and not Ps.looped)
807 A += a2
808 ta2 = _t_2(a2)
809 tdb = _t_2(db, points=i)
810 R += _at2(tdb * (ta1 + ta2),
811 _1_0 + ta1 * ta2)
812 ta1, b1 = ta2, b2
814 if not p2.isequalTo(p1, eps=EPS):
815 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
816 if z1 is not None:
817 D += _w180(z - z1) # (z - z1 + 540) ...
818 D += _w180(z2 - z) # (z2 - z + 540) % 360 - 180
819 p1, z1 = p2, z2
821 R = abs(R * _2_0)
822 if abs(D) < _90_0: # ispolar(points)
823 R = abs(R - PI2)
824 if radius:
825 a = degrees(A.fover(len(A))) # mean lat
826 R *= _mean_radius(radius, a)**2
827 return float(R)
830def _destination2(a, b, r, t):
831 '''(INTERNAL) Destination lat- and longitude in C{radians}.
833 @arg a: Latitude (C{radians}).
834 @arg b: Longitude (C{radians}).
835 @arg r: Angular distance (C{radians}).
836 @arg t: Bearing (compass C{radians}).
838 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
839 '''
840 # see <https://www.EdWilliams.org/avform.htm#LL>
841 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
842 ca *= sr
844 a = asin1(ct * ca + cr * sa)
845 d = atan2(st * ca, cr - sa * sin(a))
846 # note, in EdWilliams.org/avform.htm W is + and E is -
847 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
850def _int3d2(s, end, wrap, _i_, Vector, hs):
851 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
852 # and similar logic in .ellipsoidalBaseDI._intersect3
853 a1, b1 = s.philam
855 if _isDegrees(end): # bearing, get pseudo-end point
856 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
857 else: # must be a point
858 s.others(end, name=_end_ + _i_)
859 hs.append(end.height)
860 a2, b2 = end.philam
861 if wrap:
862 a2, b2 = _Wrap.philam(a2, b2)
864 db, b2 = unrollPI(b1, b2, wrap=wrap)
865 if max(fabs(db), fabs(a2 - a1)) < EPS:
866 raise _ValueError(_SPACE_(_line_ + _i_, _null_))
867 # note, in EdWilliams.org/avform.htm W is + and E is -
868 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
869 -(b1 + b2) * _0_5)
870 cb21 *= sin(a1 - a2) # sa21
871 sb21 *= sin(a1 + a2) # sa12
872 x = Vector(fdot_(sb12, cb21, -cb12, sb21),
873 fdot_(cb12, cb21, sb12, sb21),
874 cos(a1) * cos(a2) * sin(db)) # ll=start
875 return x.unit(), (db, (a2 - a1)) # negated d
878def _intdot(ds, a1, b1, a, b, wrap):
879 # compute dot product ds . (-b + b1, a - a1)
880 db, _ = unrollPI(b1, b, wrap=wrap)
881 return fdot(ds, db, a - a1)
884def intersecant2(center, circle, point, other, **radius_exact_height_wrap):
885 '''Compute the intersections of a circle and a (great circle) line given as
886 two points or as a point and a bearing from North.
888 @arg center: Center of the circle (L{LatLon}).
889 @arg circle: Radius of the circle (C{meter}, same units as the earth
890 B{C{radius}}) or a point on the circle (L{LatLon}).
891 @arg point: A point on the (great circle) line (L{LatLon}).
892 @arg other: An other point on the (great circle) line (L{LatLon}) or
893 the bearing at the B{C{point}} (compass C{degrees360}).
894 @kwarg radius_exact_height_wrap: Optional keyword arguments, see method
895 L{intersecant2<pygeodesy.sphericalBase.LatLonSphericalBase.
896 intersecant2>} for further details.
898 @return: 2-Tuple of the intersection points (representing a chord), each
899 an instance of the B{C{point}} class. Both points are the same
900 instance if the (great circle) line is tangent to the circle.
902 @raise IntersectionError: The circle and line do not intersect.
904 @raise TypeError: If B{C{center}}, B{C{point}}, B{C{circle}} or B{C{other}}
905 not L{LatLon}.
907 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},
908 B{C{exact}}, B{C{height}} or B{C{napieradius}}.
909 '''
910 c = _T00.others(center=center)
911 p = _T00.others(point=point)
912 try:
913 return _intersecant2(c, circle, p, other, **radius_exact_height_wrap)
914 except (TypeError, ValueError) as x:
915 raise _xError(x, center=center, circle=circle, point=point, other=other,
916 **radius_exact_height_wrap)
919def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
920 LatLon=LatLon, **LatLon_kwds):
921 # (INTERNAL) Intersect two (spherical) lines, see L{intersection}
922 # above, separated to allow callers to embellish any exceptions
924 s1, s2 = start1, start2
925 if wrap:
926 s2 = _Wrap.point(s2)
927 hs = [s1.height, s2.height]
929 a1, b1 = s1.philam
930 a2, b2 = s2.philam
931 db, b2 = unrollPI(b1, b2, wrap=wrap)
932 r12 = vincentys_(a2, a1, db)
933 if fabs(r12) < EPS: # [nearly] coincident points
934 a, b = favg(a1, a2), favg(b1, b2)
936 # see <https://www.EdWilliams.org/avform.htm#Intersection>
937 elif _isDegrees(end1) and _isDegrees(end2): # both bearings
938 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
940 x1, x2 = (sr12 * ca1), (sr12 * ca2)
941 if isnear0(x1) or isnear0(x2):
942 raise IntersectionError(_parallel_)
943 # handle domain error for equivalent longitudes,
944 # see also functions asin_safe and acos_safe at
945 # <https://www.EdWilliams.org/avform.htm#Math>
946 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1)
947 t21, t23 = acos1((sa1 - sa2 * cr12) / x2), radiansPI2(end2)
948 if sin(db) > 0:
949 t21 = PI2 - t21
950 else:
951 t12 = PI2 - t12
952 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
953 wrapPI(t21 - t23)) # angle 1-2-3)
954 if isnear0(sx1) and isnear0(sx2):
955 raise IntersectionError(_infinite_)
956 sx3 = sx1 * sx2
957# XXX if sx3 < 0:
958# XXX raise ValueError(_ambiguous_)
959 x3 = acos1(cr12 * sx3 - cx2 * cx1)
960 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
962 a, b = _destination2(a1, b1, r13, t13)
963 # like .ellipsoidalBaseDI,_intersect3, if this intersection
964 # is "before" the first point, use the antipodal intersection
965 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
966 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
968 else: # end point(s) or bearing(s)
969 _N_vector_ = _MODS.nvectorBase._N_vector_
971 x1, d1 = _int3d2(s1, end1, wrap, _1_, _N_vector_, hs)
972 x2, d2 = _int3d2(s2, end2, wrap, _2_, _N_vector_, hs)
973 x = x1.cross(x2)
974 if x.length < EPS: # [nearly] colinear or parallel lines
975 raise IntersectionError(_colinear_)
976 a, b = x.philam
977 # choose intersection similar to sphericalNvector
978 if not (_intdot(d1, a1, b1, a, b, wrap) *
979 _intdot(d2, a2, b2, a, b, wrap)) > 0:
980 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
982 h = fmean(hs) if height is None else Height(height)
983 return _LL3Tuple(degrees90(a), degrees180(b), h,
984 intersection, LatLon, LatLon_kwds)
987def intersection(start1, end1, start2, end2, height=None, wrap=False,
988 **LatLon_and_kwds):
989 '''Compute the intersection point of two lines, each defined by
990 two points or by a start point and a bearing from North.
992 @arg start1: Start point of the first line (L{LatLon}).
993 @arg end1: End point of the first line (L{LatLon}) or the bearing
994 at the first start point (compass C{degrees360}).
995 @arg start2: Start point of the second line (L{LatLon}).
996 @arg end2: End point of the second line (L{LatLon}) or the bearing
997 at the second start point (compass C{degrees360}).
998 @kwarg height: Optional height for the intersection point,
999 overriding the mean height (C{meter}).
1000 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{start2}}
1001 and both B{C{end*}} points (C{bool}).
1002 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=}L{LatLon} to use
1003 for the intersection point and optionally additional
1004 B{C{LatLon}} keyword arguments, ignored if C{B{LatLon}
1005 is None}.
1007 @return: The intersection point as a (B{C{LatLon}}) or if C{B{LatLon}
1008 is None} a L{LatLon3Tuple}C{(lat, lon, height)}. An alternate
1009 intersection point might be the L{antipode} to the returned result.
1011 @raise IntersectionError: Ambiguous or infinite intersection or colinear,
1012 parallel or otherwise non-intersecting lines.
1014 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}} or B{C{end2}}
1015 point not L{LatLon}.
1017 @raise ValueError: Invalid B{C{height}} or C{null} line.
1018 '''
1019 s1 = _T00.others(start1=start1)
1020 s2 = _T00.others(start2=start2)
1021 try:
1022 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap, **LatLon_and_kwds)
1023 except (TypeError, ValueError) as x:
1024 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1027def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1028 height=None, wrap=False, # was=True
1029 **LatLon_and_kwds):
1030 '''Compute the intersection points of two circles each defined by a
1031 center point and a radius.
1033 @arg center1: Center of the first circle (L{LatLon}).
1034 @arg rad1: Radius of the first circle (C{meter} or C{radians}, see
1035 B{C{radius}}).
1036 @arg center2: Center of the second circle (L{LatLon}).
1037 @arg rad2: Radius of the second circle (C{meter} or C{radians}, see
1038 B{C{radius}}).
1039 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1040 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1041 @kwarg eps: Required overlap (C{meter} or C{radians}, see B{C{radius}}).
1042 @kwarg height: Optional height for the intersection points (C{meter},
1043 conventionally) or C{None} for the I{"radical height"}
1044 at the I{radical line} between both centers.
1045 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{center2}}
1046 (C{bool}).
1047 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=}L{LatLon} to use for
1048 the intersection points and optionally additional B{C{LatLon}}
1049 keyword arguments, ignored if C{B{LatLon} is None}.
1051 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1052 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1053 lon, height)}. For abutting circles, both intersection
1054 points are the same instance, aka the I{radical center}.
1056 @raise IntersectionError: Concentric, antipodal, invalid or
1057 non-intersecting circles.
1059 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1061 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1062 B{C{eps}} or B{C{height}}.
1064 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1066 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1067 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1068 '''
1069 c1 = _T00.others(center1=center1)
1070 c2 = _T00.others(center2=center2)
1071 try:
1072 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1073 height=height, wrap=wrap,
1074 **LatLon_and_kwds)
1075 except (TypeError, ValueError) as x:
1076 raise _xError(x, center1=center1, rad1=rad1,
1077 center2=center2, rad2=rad2, wrap=wrap)
1080def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1081 height=None, too_d=None, wrap=False, # was=True
1082 LatLon=LatLon, **LatLon_kwds):
1083 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1084 # above, separated to allow callers to embellish any exceptions
1086 def _dest3(bearing, h):
1087 a, b = _destination2(a1, b1, r1, bearing)
1088 return _LL3Tuple(degrees90(a), degrees180(b), h,
1089 intersections2, LatLon, LatLon_kwds)
1091 a1, b1 = c1.philam
1092 a2, b2 = c2.philam
1093 if wrap:
1094 a2, b2 = _Wrap.philam(a2, b2)
1096 r1, r2, f = _rads3(rad1, rad2, radius)
1097 if f: # swapped radii, swap centers
1098 a1, a2 = a2, a1 # PYCHOK swap!
1099 b1, b2 = b2, b1 # PYCHOK swap!
1101 db, b2 = unrollPI(b1, b2, wrap=wrap)
1102 d = vincentys_(a2, a1, db) # radians
1103 if d < max(r1 - r2, EPS):
1104 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1106 r = eps if radius is None else (m2radians(
1107 eps, radius=radius) if eps else _0_0)
1108 if r < _0_0:
1109 raise _ValueError(eps=r)
1111 x = fsumf_(r1, r2, -d) # overlap
1112 if x > max(r, EPS):
1113 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1114 x = sd * sr1
1115 if isnear0(x):
1116 raise _ValueError(_invalid_)
1117 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1119 elif x < r: # PYCHOK no cover
1120 t = (d * radius) if too_d is None else too_d
1121 raise IntersectionError(_too_(_Fmt.distant(t)))
1123 if height is None: # "radical height"
1124 f = _radical2(d, r1, r2).ratio
1125 h = Height(favg(c1.height, c2.height, f=f))
1126 else:
1127 h = Height(height)
1129 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1130 if x < EPS4: # externally ...
1131 r = _dest3(b, h)
1132 elif x > _PI_EPS4: # internally ...
1133 r = _dest3(b + PI, h)
1134 else:
1135 return _dest3(b + x, h), _dest3(b - x, h)
1136 return r, r # ... abutting circles
1139@deprecated_function
1140def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1141 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1142 '''
1143 return ispolar(points, wrap=wrap)
1146def _LL3Tuple(lat, lon, height, where, LatLon, LatLon_kwds):
1147 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1148 '''
1149 n = where.__name__
1150 if LatLon is None:
1151 r = LatLon3Tuple(lat, lon, height, name=n)
1152 else:
1153 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1154 r = LatLon(lat, lon, **kwds)
1155 return r
1158def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):
1159 '''Compute the I{geographic} mean of several points.
1161 @arg points: Points to be averaged (L{LatLon}[]).
1162 @kwarg height: Optional height at mean point, overriding the mean height
1163 (C{meter}).
1164 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}} (C{bool}).
1165 @kwarg LatLon: Optional class to return the mean point (L{LatLon}) or C{None}.
1166 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments,
1167 ignored if C{B{LatLon} is None}.
1169 @return: The geographic mean and height (B{C{LatLon}}) or if C{B{LatLon}
1170 is None}, a L{LatLon3Tuple}C{(lat, lon, height)}.
1172 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1174 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1175 '''
1176 def _N_vs(ps, w):
1177 Ps = _T00.PointsIter(ps, wrap=w)
1178 for p in Ps.iterate(closed=False):
1179 yield p._N_vector
1181 m = _MODS.nvectorBase
1182 # geographic, vectorial mean
1183 n = m.sumOf(_N_vs(points, wrap), h=height, Vector=m.NvectorBase)
1184 lat, lon, h = n.latlonheight
1185 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1188@deprecated_function
1189def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1190 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1192 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1193 point (L{LatLon}) on the polygon and the C{distance}
1194 between the C{closest} and the given B{C{point}}. The
1195 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1196 lon)} if C{B{LatLon} is None} ...
1197 '''
1198 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1199 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1200 ll = LatLon2Tuple(ll.lat, ll.lon)
1201 return ll, d
1204def nearestOn3(point, points, closed=False, radius=R_M, wrap=False, adjust=True,
1205 limit=9, **LatLon_and_kwds):
1206 '''Locate the point on a path or polygon closest to a reference point.
1208 Distances are I{approximated} using function L{equirectangular4
1209 <pygeodesy.equirectangular4>}, subject to the supplied B{C{options}}.
1211 @arg point: The reference point (L{LatLon}).
1212 @arg points: The path or polygon points (L{LatLon}[]).
1213 @kwarg closed: Optionally, close the polygon (C{bool}).
1214 @kwarg radius: Mean earth radius (C{meter}).
1215 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1216 B{C{points}} (C{bool}).
1217 @kwarg adjust: See function L{equirectangular4<pygeodesy.equirectangular4>}
1218 (C{bool}).
1219 @kwarg limit: See function L{equirectangular4<pygeodesy.equirectangular4>}
1220 (C{degrees}), default C{9 degrees} is about C{1,000 Km} (for
1221 (mean spherical earth radius L{R_KM}).
1222 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=L{LatLon}} to return the
1223 closest point and optionally additional C{B{LatLon}} keyword
1224 arguments or specify C{B{LatLon}=None}.
1226 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1227 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1228 lon, height)} if C{B{LatLon} is None}. The C{distance} is
1229 the L{equirectangular4<pygeodesy.equirectangular4>} distance
1230 between the C{closest} and the given B{C{point}} converted to
1231 C{meter}, same units as B{C{radius}}. The C{angle} from the
1232 given B{C{point}} to the C{closest} is in compass C{degrees360},
1233 like function L{compassAngle<pygeodesy.compassAngle>}. The
1234 C{height} is the (interpolated) height at the C{closest} point.
1236 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1237 see function L{equirectangular4<pygeodesy.equirectangular4>}.
1239 @raise PointsError: Insufficient number of B{C{points}}.
1241 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1243 @raise ValueError: Invalid B{C{radius}}.
1245 @see: Functions L{equirectangular4<pygeodesy.equirectangular4>} and
1246 L{nearestOn5<pygeodesy.nearestOn5>}.
1247 '''
1248 t = _nearestOn5(point, points, closed=closed, wrap=wrap,
1249 adjust=adjust, limit=limit)
1250 d = degrees2m(t.distance, radius=radius)
1251 h = t.height
1252 n = nearestOn3.__name__
1254 LL, kwds = _xkwds_pop2(LatLon_and_kwds, LatLon=LatLon)
1255 r = LatLon3Tuple(t.lat, t.lon, h, name=n) if LL is None else \
1256 LL(t.lat, t.lon, **_xkwds(kwds, height=h, name=n))
1257 return NearestOn3Tuple(r, d, t.angle, name=n)
1260def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1261 '''Compute the perimeter of a (spherical) polygon or composite
1262 (with great circle arcs joining the points).
1264 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1265 or L{BooleanGH}).
1266 @kwarg closed: Optionally, close the polygon (C{bool}).
1267 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1268 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1269 B{C{points}} (C{bool}).
1271 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1272 or C{radians} if C{B{radius} is None}).
1274 @raise PointsError: Insufficient number of B{C{points}}.
1276 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1278 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1279 C{B{points}} a composite.
1281 @note: Distances are based on function L{vincentys_<pygeodesy.vincentys_>}.
1283 @see: Functions L{perimeterOf<pygeodesy.perimeterOf>},
1284 L{sphericalNvector.perimeterOf} and L{ellipsoidalKarney.perimeterOf}.
1285 '''
1286 def _rads(ps, c, w): # angular edge lengths in radians
1287 Ps = _T00.PointsIter(ps, loop=1, wrap=w)
1288 a1, b1 = Ps[0].philam
1289 for p in Ps.iterate(closed=c):
1290 a2, b2 = p.philam
1291 db, b2 = unrollPI(b1, b2, wrap=w and not (c and Ps.looped))
1292 yield vincentys_(a2, a1, db)
1293 a1, b1 = a2, b2
1295 if _MODS.booleans.isBoolean(points):
1296 if not closed:
1297 raise _ValueError(closed=closed, points=_composite_)
1298 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1299 else:
1300 r = fsum(_rads(points, closed, wrap))
1301 return _radians2m(r, radius)
1304def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1305 excess=excessAbc_,
1306 wrap=False):
1307 '''Compute the angles, sides, and area of a (spherical) triangle.
1309 @arg latA: First corner latitude (C{degrees}).
1310 @arg lonA: First corner longitude (C{degrees}).
1311 @arg latB: Second corner latitude (C{degrees}).
1312 @arg lonB: Second corner longitude (C{degrees}).
1313 @arg latC: Third corner latitude (C{degrees}).
1314 @arg lonC: Third corner longitude (C{degrees}).
1315 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1316 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1317 or C{None}.
1318 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1319 L{excessGirard_} or L{excessLHuilier_}).
1320 @kwarg wrap: If C{True}, wrap and L{unroll180<pygeodesy.unroll180>}
1321 longitudes (C{bool}).
1323 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1324 spherical angles C{A}, C{B} and C{C}, angular sides
1325 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1326 in I{square} C{meter} or same units as B{C{radius}}
1327 I{squared} or if C{B{radius}=0} or C{None}, a
1328 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1329 I{spherical deficit} C{D} and I{spherical excess}
1330 C{E} as the C{unit area}, all in C{radians}.
1331 '''
1332 t = triangle8_(Phid(latA=latA), Lamd(lonA=lonA),
1333 Phid(latB=latB), Lamd(lonB=lonB),
1334 Phid(latC=latC), Lamd(lonC=lonC),
1335 excess=excess, wrap=wrap)
1336 return _t7Tuple(t, radius)
1339def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_,
1340 wrap=False):
1341 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1342 excess} of a (spherical) triangle.
1344 @arg phiA: First corner latitude (C{radians}).
1345 @arg lamA: First corner longitude (C{radians}).
1346 @arg phiB: Second corner latitude (C{radians}).
1347 @arg lamB: Second corner longitude (C{radians}).
1348 @arg phiC: Third corner latitude (C{radians}).
1349 @arg lamC: Third corner longitude (C{radians}).
1350 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1351 L{excessGirard_} or L{excessLHuilier_}).
1352 @kwarg wrap: If C{True}, L{unrollPI<pygeodesy.unrollPI>} the
1353 longitudinal deltas (C{bool}).
1355 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1356 spherical angles C{A}, C{B} and C{C}, angular sides
1357 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1358 I{spherical excess} C{E}, all in C{radians}.
1359 '''
1360 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1361 d, _ = unrollPI(lamB, lamC, wrap=w)
1362 a = vincentys_(phiC, phiB, d)
1363 return a, (phiB, lamB, phiC, lamC, phiA, lamA) # rotate A, B, C
1365 def _A_r(a, sa, ca, sb, cb, sc, cc):
1366 s = sb * sc
1367 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1368 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2_'s
1370 # notation: side C{a} is oposite to corner C{A}, etc.
1371 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1372 b, r = _a_r(wrap, *r)
1373 c, _ = _a_r(wrap, *r)
1375 A, r = _A_r(a, *sincos2_(a, b, c))
1376 B, r = _A_r(b, *r)
1377 C, _ = _A_r(c, *r)
1379 D = fsumf_(PI2, -a, -b, -c) # deficit aka defect
1380 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else (
1381 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else
1382 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b))))
1384 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1387def _t7Tuple(t, radius):
1388 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1389 '''
1390 if radius: # not in (None, _0_0)
1391 r = radius if _isRadius(radius) else \
1392 _ellipsoidal_datum(radius).ellipsoid.Rmean
1393 A, B, C = map1(degrees, t.A, t.B, t.C)
1394 t = Triangle7Tuple(A, (r * t.a),
1395 B, (r * t.b),
1396 C, (r * t.c), t.E * r**2)
1397 return t
1400__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1401 areaOf, # functions
1402 intersecant2, intersection, intersections2, ispolar,
1403 isPoleEnclosedBy, # DEPRECATED, use ispolar
1404 meanOf,
1405 nearestOn2, nearestOn3,
1406 perimeterOf,
1407 sumOf, # XXX == vector3d.sumOf
1408 triangle7, triangle8_)
1410# **) MIT License
1411#
1412# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1413#
1414# Permission is hereby granted, free of charge, to any person obtaining a
1415# copy of this software and associated documentation files (the "Software"),
1416# to deal in the Software without restriction, including without limitation
1417# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1418# and/or sell copies of the Software, and to permit persons to whom the
1419# Software is furnished to do so, subject to the following conditions:
1420#
1421# The above copyright notice and this permission notice shall be included
1422# in all copies or substantial portions of the Software.
1423#
1424# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1425# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1426# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1427# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1428# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1429# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1430# OTHER DEALINGS IN THE SOFTWARE.