Coverage for pygeodesy/sphericalTrigonometry.py: 93%
386 statements
« prev ^ index » next coverage.py v7.6.1, created at 2024-10-22 18:16 -0400
« prev ^ index » next coverage.py v7.6.1, created at 2024-10-22 18:16 -0400
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, 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.10.12'
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
186 the great circle defined by a start and 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 the given initial bearing.
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 on the given initial bearing from this point.
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 return Vector(sb * ct - cb * sa * st,
275 -cb * ct - sb * sa * st,
276 ca * st, **Vector_kwds) # XXX .unit()?
278 def initialBearingTo(self, other, wrap=False, raiser=False):
279 '''Compute the initial bearing (forward azimuth) from this
280 to an other point.
282 @arg other: The other point (spherical L{LatLon}).
283 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
284 the B{C{other}} point (C{bool}).
285 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}),
286 use C{B{raiser}=True} for behavior like
287 C{sphericalNvector.LatLon.initialBearingTo}.
289 @return: Initial bearing (compass C{degrees360}).
291 @raise CrossError: If this and the B{C{other}} point coincide
292 and if B{C{raiser}} and L{crosserrors
293 <pygeodesy.crosserrors>} are both C{True}.
295 @raise TypeError: The B{C{other}} point is not L{LatLon}.
296 '''
297 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
298 # XXX behavior like sphericalNvector.LatLon.initialBearingTo
299 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(db)) < EPS:
300 raise CrossError(_point_, self, other=other, wrap=wrap, txt=_coincident_)
302 return degrees(bearing_(a1, b1, a2, b2, final=False))
304 def intermediateTo(self, other, fraction, height=None, wrap=False):
305 '''Locate the point at given fraction between (or along) this
306 and an other point.
308 @arg other: The other point (L{LatLon}).
309 @arg fraction: Fraction between both points (C{scalar},
310 0.0 at this and 1.0 at the other point).
311 @kwarg height: Optional height, overriding the intermediate
312 height (C{meter}).
313 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
314 B{C{other}} point (C{bool}).
316 @return: Intermediate point (L{LatLon}).
318 @raise TypeError: The B{C{other}} point is not L{LatLon}.
320 @raise ValueError: Invalid B{C{fraction}} or B{C{height}}.
322 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
323 '''
324 p = self
325 f = Scalar(fraction=fraction)
326 if not isnear0(f):
327 p = p.others(other)
328 if wrap:
329 p = _Wrap.point(p)
330 if not isnear1(f): # and not near0
331 a1, b1 = self.philam
332 a2, b2 = p.philam
333 db, b2 = unrollPI(b1, b2, wrap=wrap)
334 r = vincentys_(a2, a1, db)
335 sr = sin(r)
336 if isnon0(sr):
337 sa1, ca1, sa2, ca2, \
338 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2)
340 t = f * r
341 a = sin(r - t) # / sr superflous
342 b = sin( t) # / sr superflous
344 x = a * ca1 * cb1 + b * ca2 * cb2
345 y = a * ca1 * sb1 + b * ca2 * sb2
346 z = a * sa1 + b * sa2
348 a = atan1d(z, hypot(x, y))
349 b = atan2d(y, x)
351 else: # PYCHOK no cover
352 a = degrees90( favg(a1, a2, f=f)) # coincident
353 b = degrees180(favg(b1, b2, f=f))
355 h = self._havg(other, f=f, h=height)
356 p = self.classof(a, b, height=h)
357 return p
359 def intersection(self, end1, other, end2, height=None, wrap=False):
360 '''Compute the intersection point of two lines, each defined by
361 two points or a start point and bearing from North.
363 @arg end1: End point of this line (L{LatLon}) or the initial
364 bearing at this point (compass C{degrees360}).
365 @arg other: Start point of the other line (L{LatLon}).
366 @arg end2: End point of the other line (L{LatLon}) or the
367 initial bearing at the B{C{other}} point (compass
368 C{degrees360}).
369 @kwarg height: Optional height for intersection point,
370 overriding the mean height (C{meter}).
371 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
372 B{C{start2}} and both B{C{end*}} points (C{bool}).
374 @return: The intersection point (L{LatLon}). An alternate
375 intersection point might be the L{antipode} to
376 the returned result.
378 @raise IntersectionError: Ambiguous or infinite intersection
379 or colinear, parallel or otherwise
380 non-intersecting lines.
382 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}}
383 or B{C{end2}} not C{scalar} nor L{LatLon}.
385 @raise ValueError: Invalid B{C{height}} or C{null} line.
386 '''
387 try:
388 s2 = self.others(other)
389 return _intersect(self, end1, s2, end2, height=height, wrap=wrap,
390 LatLon=self.classof)
391 except (TypeError, ValueError) as x:
392 raise _xError(x, start1=self, end1=end1,
393 other=other, end2=end2, wrap=wrap)
395 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0,
396 height=None, wrap=True):
397 '''Compute the intersection points of two circles, each defined
398 by a center point and radius.
400 @arg rad1: Radius of the this circle (C{meter} or C{radians},
401 see B{C{radius}}).
402 @arg other: Center point of the other circle (L{LatLon}).
403 @arg rad2: Radius of the other circle (C{meter} or C{radians},
404 see B{C{radius}}).
405 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
406 B{C{rad2}} and B{C{eps}} are given in C{radians}).
407 @kwarg eps: Required overlap (C{meter} or C{radians}, see
408 B{C{radius}}).
409 @kwarg height: Optional height for the intersection points (C{meter},
410 conventionally) or C{None} for the I{"radical height"}
411 at the I{radical line} between both centers.
412 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
413 B{C{other}} point (C{bool}).
415 @return: 2-Tuple of the intersection points, each a L{LatLon}
416 instance. For abutting circles, both intersection
417 points are the same instance, aka the I{radical center}.
419 @raise IntersectionError: Concentric, antipodal, invalid or
420 non-intersecting circles.
422 @raise TypeError: If B{C{other}} is not L{LatLon}.
424 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
425 B{C{eps}} or B{C{height}}.
426 '''
427 try:
428 c2 = self.others(other)
429 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps,
430 height=height, wrap=wrap,
431 LatLon=self.classof)
432 except (TypeError, ValueError) as x:
433 raise _xError(x, center=self, rad1=rad1,
434 other=other, rad2=rad2, wrap=wrap)
436 @deprecated_method
437 def isEnclosedBy(self, points): # PYCHOK no cover
438 '''DEPRECATED, use method C{isenclosedBy}.'''
439 return self.isenclosedBy(points)
441 def isenclosedBy(self, points, wrap=False):
442 '''Check whether a (convex) polygon or composite encloses this point.
444 @arg points: The polygon points or composite (L{LatLon}[],
445 L{BooleanFHP} or L{BooleanGH}).
446 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
447 B{C{points}} (C{bool}).
449 @return: C{True} if this point is inside the polygon or
450 composite, C{False} otherwise.
452 @raise PointsError: Insufficient number of B{C{points}}.
454 @raise TypeError: Some B{C{points}} are not L{LatLon}.
456 @raise ValueError: Invalid B{C{points}}, non-convex polygon.
458 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
459 and L{pygeodesy.ispolar} especially if the B{C{points}} may
460 enclose a pole or wrap around the earth I{longitudinally}.
461 '''
462 if _MODS.booleans.isBoolean(points):
463 return points._encloses(self.lat, self.lon, wrap=wrap)
465 Ps = self.PointsIter(points, loop=2, dedup=True, wrap=wrap)
466 n0 = self._N_vector
468 v2 = Ps[0]._N_vector
469 p1 = Ps[1]
470 v1 = p1._N_vector
471 # check whether this point on same side of all
472 # polygon edges (to the left or right depending
473 # on the anti-/clockwise polygon direction)
474 gc1 = v2.cross(v1)
475 t0 = gc1.angleTo(n0) > PI_2
476 s0 = None
477 # get great-circle vector for each edge
478 for i, p2 in Ps.enumerate(closed=True):
479 if wrap and not Ps.looped:
480 p2 = _unrollon(p1, p2)
481 p1 = p2
482 v2 = p2._N_vector
483 gc = v1.cross(v2)
484 t = gc.angleTo(n0) > PI_2
485 if t != t0: # different sides of edge i
486 return False # outside
488 # check for convex polygon: angle between
489 # gc vectors, signed by direction of n0
490 # (otherwise the test above is not reliable)
491 s = signOf(gc1.angleTo(gc, vSign=n0))
492 if s != s0:
493 if s0 is None:
494 s0 = s
495 else:
496 t = _Fmt.SQUARE(points=i)
497 raise _ValueError(t, p2, wrap=wrap, txt_not_=_convex_)
498 gc1, v1 = gc, v2
500 return True # inside
502 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
503 '''Find the midpoint between this and an other point.
505 @arg other: The other point (L{LatLon}).
506 @kwarg height: Optional height for midpoint, overriding
507 the mean height (C{meter}).
508 @kwarg fraction: Midpoint location from this point (C{scalar}),
509 may be negative or greater than 1.0.
510 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
511 B{C{other}} point (C{bool}).
513 @return: Midpoint (L{LatLon}).
515 @raise TypeError: The B{C{other}} point is not L{LatLon}.
517 @raise ValueError: Invalid B{C{height}}.
519 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
520 '''
521 if fraction is _0_5:
522 # see <https://MathForum.org/library/drmath/view/51822.html>
523 a1, b, a2, _, db = self._ab1_ab2_db5(other, wrap)
524 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db)
526 x = ca2 * cdb + ca1
527 y = ca2 * sdb
529 a = atan1d(sa1 + sa2, hypot(x, y))
530 b = degrees180(b + atan2(y, x))
532 h = self._havg(other, h=height)
533 r = self.classof(a, b, height=h)
534 else:
535 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
536 return r
538 def nearestOn(self, point1, point2, radius=R_M, **wrap_adjust_limit):
539 '''Locate the point between two points closest to this point.
541 Distances are approximated by function L{pygeodesy.equirectangular4},
542 subject to the supplied B{C{options}}.
544 @arg point1: Start point (L{LatLon}).
545 @arg point2: End point (L{LatLon}).
546 @kwarg radius: Mean earth radius (C{meter}).
547 @kwarg wrap_adjust_limit: Optional keyword arguments for functions
548 L{sphericalTrigonometry.nearestOn3} and
549 L{pygeodesy.equirectangular4},
551 @return: Closest point on the great circle line (L{LatLon}).
553 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
554 see function L{pygeodesy.equirectangular4}.
556 @raise NotImplementedError: Keyword argument C{B{within}=False}
557 is not (yet) supported.
559 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
561 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
563 @see: Functions L{pygeodesy.equirectangular4} and L{pygeodesy.nearestOn5}
564 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
565 '''
566 # remove kwarg B{C{within}} if present
567 w, kwds = _xkwds_pop2(wrap_adjust_limit, within=True)
568 if not w:
569 self._notImplemented(within=w)
571# # UNTESTED - handle C{B{within}=False} and C{B{within}=True}
572# wrap = _xkwds_get(options, wrap=False)
573# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap)
574# if fabs(a) < EPS or (within and a < EPS):
575# return point1
576# d = point1.distanceTo(point2, radius=radius, wrap=wrap)
577# if isnear0(d):
578# return point1 # or point2
579# elif fabs(d - a) < EPS or (a + EPS) > d:
580# return point2
581# f = a / d
582# if within:
583# if f > EPS1:
584# return point2
585# elif f < EPS:
586# return point1
587# return point1.intermediateTo(point2, f, wrap=wrap)
589 # without kwarg B{C{within}}, use backward compatible .nearestOn3
590 return self.nearestOn3([point1, point2], closed=False, radius=radius,
591 **kwds)[0]
593 @deprecated_method
594 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover
595 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
597 @return: ... 2-Tuple C{(closest, distance)} of the closest
598 point (L{LatLon}) on the polygon and the distance
599 to that point from this point in C{meter}, same
600 units of B{C{radius}}.
601 '''
602 r = self.nearestOn3(points, closed=closed, radius=radius, **options)
603 return r.closest, r.distance
605 def nearestOn3(self, points, closed=False, radius=R_M, **wrap_adjust_limit):
606 '''Locate the point on a polygon closest to this point.
608 Distances are approximated by function L{pygeodesy.equirectangular4},
609 subject to the supplied B{C{options}}.
611 @arg points: The polygon points (L{LatLon}[]).
612 @kwarg closed: Optionally, close the polygon (C{bool}).
613 @kwarg radius: Mean earth radius (C{meter}).
614 @kwarg wrap_adjust_limit: Optional keyword arguments for function
615 L{sphericalTrigonometry.nearestOn3} and
616 L{pygeodesy.equirectangular4},
618 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
619 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular4}
620 C{distance} between this and the C{closest} point converted to
621 C{meter}, same units as B{C{radius}}. The C{angle} from this
622 to the C{closest} point is in compass C{degrees360}, like
623 function L{pygeodesy.compassAngle}.
625 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
626 see function L{pygeodesy.equirectangular4}.
628 @raise PointsError: Insufficient number of B{C{points}}.
630 @raise TypeError: Some B{C{points}} are not C{LatLon}.
632 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
634 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular4}
635 and L{pygeodesy.nearestOn5}.
636 '''
637 return nearestOn3(self, points, closed=closed, radius=radius,
638 LatLon=self.classof, **wrap_adjust_limit)
640 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
641 '''Convert this point to C{Karney}-based cartesian (ECEF) coordinates.
643 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}} and other
644 keyword arguments, ignored if C{B{Cartesian} is
645 None}. Use C{B{Cartesian}=...} to override this
646 L{Cartesian} class or specify C{B{Cartesian}=None}.
648 @return: The cartesian point (L{Cartesian}) or if C{B{Cartesian} is None},
649 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with C{C}
650 and C{M} if available.
652 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
653 '''
654 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
655 return LatLonSphericalBase.toCartesian(self, **kwds)
657 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
658 '''Compute the angles, sides and area of a spherical triangle.
660 @arg otherB: Second triangle point (C{LatLon}).
661 @arg otherC: Third triangle point (C{LatLon}).
662 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, L{Ellipsoid},
663 L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) or C{None}.
664 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll points B{C{otherB}}
665 and B{C{otherC}} (C{bool}).
667 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if B{C{radius} is
668 None}, a L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)}.
670 @see: Function L{triangle7} and U{Spherical trigonometry
671 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
672 '''
673 B = self.others(otherB=otherB)
674 C = self.others(otherC=otherC)
675 B, C, _ = _unrollon3(self, B, C, wrap)
677 r = self.philam + B.philam + C.philam
678 t = triangle8_(*r, wrap=wrap)
679 return self._xnamed(_t7Tuple(t, radius))
681 def triangulate(self, bearing1, other, bearing2, **height_wrap):
682 '''Locate a point given this, an other point and the (initial) bearing
683 at this and at the other point.
685 @arg bearing1: Bearing at this point (compass C{degrees360}).
686 @arg other: The other point (C{LatLon}).
687 @arg bearing2: Bearing at the other point (compass C{degrees360}).
688 @kwarg height_wrap_tol: Optional keyword arguments C{B{height}=None},
689 C{B{wrap}=False}, see method L{intersection}.
691 @return: Triangulated point (C{LatLon}).
693 @see: Method L{intersection} for further details.
694 '''
695 if _isDegrees(bearing1) and _isDegrees(bearing2):
696 return self.intersection(bearing1, other, bearing2, **height_wrap)
697 raise _TypeError(bearing1=bearing1, bearing2=bearing2, **height_wrap)
699 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
700 area=True, eps=EPS1, radius=R_M, wrap=False):
701 '''Trilaterate three points by I{area overlap} or I{perimeter intersection}
702 of three corresponding circles.
704 @arg distance1: Distance to this point (C{meter}, same units as B{C{radius}}).
705 @arg point2: Second center point (C{LatLon}).
706 @arg distance2: Distance to point2 (C{meter}, same units as B{C{radius}}).
707 @arg point3: Third center point (C{LatLon}).
708 @arg distance3: Distance to point3 (C{meter}, same units as B{C{radius}}).
709 @kwarg area: If C{True}, compute the area overlap, otherwise the perimeter
710 intersection of the circles (C{bool}).
711 @kwarg eps: The required I{minimal overlap} for C{B{area}=True} or the
712 I{intersection margin} if C{B{area}=False} (C{meter}, same
713 units as B{C{radius}}).
714 @kwarg radius: Mean earth radius (C{meter}, conventionally).
715 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{point2}} and
716 B{C{point3}} (C{bool}).
718 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)} with
719 C{min} and C{max} in C{meter}, same units as B{C{eps}}, the
720 corresponding trilaterated points C{minPoint} and C{maxPoint}
721 as I{spherical} C{LatLon} and C{n}, the number of trilatered
722 points found for the given B{C{eps}}.
724 If only a single trilaterated point is found, C{min I{is} max},
725 C{minPoint I{is} maxPoint} and C{n = 1}.
727 For C{B{area}=True}, C{min} and C{max} are the smallest respectively
728 largest I{radial} overlap found.
730 For C{B{area}=False}, C{min} and C{max} represent the nearest
731 respectively farthest intersection margin.
733 If C{B{area}=True} and all 3 circles are concentric, C{n=0} and
734 C{minPoint} and C{maxPoint} are both the B{C{point#}} with the
735 smallest B{C{distance#}} C{min} and C{max} the largest B{C{distance#}}.
737 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
738 insufficient overlap for C{B{area}=True} or
739 no intersection or all (near-)concentric if
740 C{B{area}=False}.
742 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
744 @raise ValueError: Coincident B{C{point2}} or B{C{point3}} or invalid
745 B{C{distance1}}, B{C{distance2}}, B{C{distance3}}
746 or B{C{radius}}.
747 '''
748 return _trilaterate5(self, distance1,
749 self.others(point2=point2), distance2,
750 self.others(point3=point3), distance3,
751 area=area, radius=radius, eps=eps, wrap=wrap)
754_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
757def areaOf(points, radius=R_M, wrap=False): # was=True
758 '''Calculate the area of a (spherical) polygon or composite
759 (with the pointsjoined by great circle arcs).
761 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
762 or L{BooleanGH}).
763 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
764 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
765 or C{None}.
766 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{points}}
767 (C{bool}).
769 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
770 or C{radians} if C{B{radius} is None}).
772 @raise PointsError: Insufficient number of B{C{points}}.
774 @raise TypeError: Some B{C{points}} are not L{LatLon}.
776 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
778 @note: The area is based on I{Karney}'s U{'Area of a spherical
779 polygon'<https://MathOverflow.net/questions/97711/
780 the-area-of-spherical-polygons>}, 3rd Answer.
782 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
783 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
784 L{ellipsoidalKarney.areaOf}.
785 '''
786 if _MODS.booleans.isBoolean(points):
787 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
789 _at2, _t_2 = atan2, tan_2
790 _un, _w180 = unrollPI, wrap180
792 Ps = _T00.PointsIter(points, loop=1, wrap=wrap)
793 p1 = p2 = Ps[0]
794 a1, b1 = p1.philam
795 ta1, z1 = _t_2(a1), None
797 A = Fsum() # mean phi
798 R = Fsum() # see L{pygeodesy.excessKarney_}
799 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
800 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
801 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
802 D = Fsum()
803 for i, p2 in Ps.enumerate(closed=True):
804 a2, b2 = p2.philam
805 db, b2 = _un(b1, b2, wrap=wrap and not Ps.looped)
806 A += a2
807 ta2 = _t_2(a2)
808 tdb = _t_2(db, points=i)
809 R += _at2(tdb * (ta1 + ta2),
810 _1_0 + ta1 * ta2)
811 ta1, b1 = ta2, b2
813 if not p2.isequalTo(p1, eps=EPS):
814 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
815 if z1 is not None:
816 D += _w180(z - z1) # (z - z1 + 540) ...
817 D += _w180(z2 - z) # (z2 - z + 540) % 360 - 180
818 p1, z1 = p2, z2
820 R = abs(R * _2_0)
821 if abs(D) < _90_0: # ispolar(points)
822 R = abs(R - PI2)
823 if radius:
824 a = degrees(A.fover(len(A))) # mean lat
825 R *= _mean_radius(radius, a)**2
826 return float(R)
829def _destination2(a, b, r, t):
830 '''(INTERNAL) Destination lat- and longitude in C{radians}.
832 @arg a: Latitude (C{radians}).
833 @arg b: Longitude (C{radians}).
834 @arg r: Angular distance (C{radians}).
835 @arg t: Bearing (compass C{radians}).
837 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
838 '''
839 # see <https://www.EdWilliams.org/avform.htm#LL>
840 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
841 ca *= sr
843 a = asin1(ct * ca + cr * sa)
844 d = atan2(st * ca, cr - sa * sin(a))
845 # note, in EdWilliams.org/avform.htm W is + and E is -
846 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
849def _int3d2(s, end, wrap, _i_, Vector, hs):
850 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
851 # and similar logic in .ellipsoidalBaseDI._intersect3
852 a1, b1 = s.philam
854 if _isDegrees(end): # bearing, get pseudo-end point
855 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
856 else: # must be a point
857 s.others(end, name=_end_ + _i_)
858 hs.append(end.height)
859 a2, b2 = end.philam
860 if wrap:
861 a2, b2 = _Wrap.philam(a2, b2)
863 db, b2 = unrollPI(b1, b2, wrap=wrap)
864 if max(fabs(db), fabs(a2 - a1)) < EPS:
865 raise _ValueError(_SPACE_(_line_ + _i_, _null_))
866 # note, in EdWilliams.org/avform.htm W is + and E is -
867 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
868 -(b1 + b2) * _0_5)
869 cb21 *= sin(a1 - a2) # sa21
870 sb21 *= sin(a1 + a2) # sa12
871 x = Vector(sb12 * cb21 - cb12 * sb21,
872 cb12 * cb21 + sb12 * sb21,
873 cos(a1) * cos(a2) * sin(db)) # ll=start
874 return x.unit(), (db, (a2 - a1)) # negated d
877def _intdot(ds, a1, b1, a, b, wrap):
878 # compute dot product ds . (-b + b1, a - a1)
879 db, _ = unrollPI(b1, b, wrap=wrap)
880 return fdot(ds, db, a - a1)
883def intersecant2(center, circle, point, other, **radius_exact_height_wrap):
884 '''Compute the intersections of a circle and a (great circle) line given as
885 two points or as a point and bearing.
887 @arg center: Center of the circle (L{LatLon}).
888 @arg circle: Radius of the circle (C{meter}, same units as the earth
889 B{C{radius}}) or a point on the circle (L{LatLon}).
890 @arg point: A point on the (great circle) line (L{LatLon}).
891 @arg other: An other point on the (great circle) line (L{LatLon}) or
892 the bearing at the B{C{point}} (compass C{degrees360}).
893 @kwarg radius_exact_height_wrap: Optional keyword arguments, see method
894 L{intersecant2<pygeodesy.sphericalBase.LatLonSphericalBase.
895 intersecant2>} for further details.
897 @return: 2-Tuple of the intersection points (representing a chord), each
898 an instance of the B{C{point}} class. Both points are the same
899 instance if the (great circle) line is tangent to the circle.
901 @raise IntersectionError: The circle and line do not intersect.
903 @raise TypeError: If B{C{center}}, B{C{point}}, B{C{circle}} or B{C{other}}
904 not L{LatLon}.
906 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},
907 B{C{exact}}, B{C{height}} or B{C{napieradius}}.
908 '''
909 c = _T00.others(center=center)
910 p = _T00.others(point=point)
911 try:
912 return _intersecant2(c, circle, p, other, **radius_exact_height_wrap)
913 except (TypeError, ValueError) as x:
914 raise _xError(x, center=center, circle=circle, point=point, other=other,
915 **radius_exact_height_wrap)
918def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
919 LatLon=None, **LatLon_kwds):
920 # (INTERNAL) Intersect two (spherical) lines, see L{intersection}
921 # above, separated to allow callers to embellish any exceptions
923 s1, s2 = start1, start2
924 if wrap:
925 s2 = _Wrap.point(s2)
926 hs = [s1.height, s2.height]
928 a1, b1 = s1.philam
929 a2, b2 = s2.philam
930 db, b2 = unrollPI(b1, b2, wrap=wrap)
931 r12 = vincentys_(a2, a1, db)
932 if fabs(r12) < EPS: # [nearly] coincident points
933 a, b = favg(a1, a2), favg(b1, b2)
935 # see <https://www.EdWilliams.org/avform.htm#Intersection>
936 elif _isDegrees(end1) and _isDegrees(end2): # both bearings
937 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
939 x1, x2 = (sr12 * ca1), (sr12 * ca2)
940 if isnear0(x1) or isnear0(x2):
941 raise IntersectionError(_parallel_)
942 # handle domain error for equivalent longitudes,
943 # see also functions asin_safe and acos_safe at
944 # <https://www.EdWilliams.org/avform.htm#Math>
945 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1)
946 t21, t23 = acos1((sa1 - sa2 * cr12) / x2), radiansPI2(end2)
947 if sin(db) > 0:
948 t21 = PI2 - t21
949 else:
950 t12 = PI2 - t12
951 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
952 wrapPI(t21 - t23)) # angle 1-2-3)
953 if isnear0(sx1) and isnear0(sx2):
954 raise IntersectionError(_infinite_)
955 sx3 = sx1 * sx2
956# XXX if sx3 < 0:
957# XXX raise ValueError(_ambiguous_)
958 x3 = acos1(cr12 * sx3 - cx2 * cx1)
959 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
961 a, b = _destination2(a1, b1, r13, t13)
962 # like .ellipsoidalBaseDI,_intersect3, if this intersection
963 # is "before" the first point, use the antipodal intersection
964 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
965 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
967 else: # end point(s) or bearing(s)
968 _N_vector_ = _MODS.nvectorBase._N_vector_
970 x1, d1 = _int3d2(s1, end1, wrap, _1_, _N_vector_, hs)
971 x2, d2 = _int3d2(s2, end2, wrap, _2_, _N_vector_, hs)
972 x = x1.cross(x2)
973 if x.length < EPS: # [nearly] colinear or parallel lines
974 raise IntersectionError(_colinear_)
975 a, b = x.philam
976 # choose intersection similar to sphericalNvector
977 if not (_intdot(d1, a1, b1, a, b, wrap) *
978 _intdot(d2, a2, b2, a, b, wrap)) > 0:
979 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
981 h = fmean(hs) if height is None else Height(height)
982 return _LL3Tuple(degrees90(a), degrees180(b), h,
983 intersection, LatLon, LatLon_kwds)
986def intersection(start1, end1, start2, end2, height=None, wrap=False,
987 LatLon=LatLon, **LatLon_kwds):
988 '''Compute the intersection point of two lines, each defined
989 by two points or a start point and bearing from North.
991 @arg start1: Start point of the first line (L{LatLon}).
992 @arg end1: End point of the first line (L{LatLon}) or
993 the initial bearing at the first start point
994 (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
997 the initial bearing at the second start point
998 (compass C{degrees360}).
999 @kwarg height: Optional height for the intersection point,
1000 overriding the mean height (C{meter}).
1001 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1002 B{C{start2}} and both B{C{end*}} points (C{bool}).
1003 @kwarg LatLon: Optional class to return the intersection
1004 point (L{LatLon}) or C{None}.
1005 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1006 arguments, ignored if C{B{LatLon} is None}.
1008 @return: The intersection point as a (B{C{LatLon}}) or if
1009 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon,
1010 height)}. An alternate intersection point might
1011 be the L{antipode} to the returned result.
1013 @raise IntersectionError: Ambiguous or infinite intersection
1014 or colinear, parallel or otherwise
1015 non-intersecting lines.
1017 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}}
1018 or B{C{end2}} point not L{LatLon}.
1020 @raise ValueError: Invalid B{C{height}} or C{null} line.
1021 '''
1022 s1 = _T00.others(start1=start1)
1023 s2 = _T00.others(start2=start2)
1024 try:
1025 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap,
1026 LatLon=LatLon, **LatLon_kwds)
1027 except (TypeError, ValueError) as x:
1028 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1031def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1032 height=None, wrap=False, # was=True
1033 LatLon=LatLon, **LatLon_kwds):
1034 '''Compute the intersection points of two circles each defined
1035 by a center point and a radius.
1037 @arg center1: Center of the first circle (L{LatLon}).
1038 @arg rad1: Radius of the first circle (C{meter} or C{radians},
1039 see B{C{radius}}).
1040 @arg center2: Center of the second circle (L{LatLon}).
1041 @arg rad2: Radius of the second circle (C{meter} or C{radians},
1042 see B{C{radius}}).
1043 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1044 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1045 @kwarg eps: Required overlap (C{meter} or C{radians}, see
1046 B{C{radius}}).
1047 @kwarg height: Optional height for the intersection points (C{meter},
1048 conventionally) or C{None} for the I{"radical height"}
1049 at the I{radical line} between both centers.
1050 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{center2}}
1051 (C{bool}).
1052 @kwarg LatLon: Optional class to return the intersection
1053 points (L{LatLon}) or C{None}.
1054 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1055 arguments, ignored if C{B{LatLon} is None}.
1057 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1058 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1059 lon, height)}. For abutting circles, both intersection
1060 points are the same instance, aka the I{radical center}.
1062 @raise IntersectionError: Concentric, antipodal, invalid or
1063 non-intersecting circles.
1065 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1067 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1068 B{C{eps}} or B{C{height}}.
1070 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1072 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1073 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1074 '''
1075 c1 = _T00.others(center1=center1)
1076 c2 = _T00.others(center2=center2)
1077 try:
1078 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1079 height=height, wrap=wrap,
1080 LatLon=LatLon, **LatLon_kwds)
1081 except (TypeError, ValueError) as x:
1082 raise _xError(x, center1=center1, rad1=rad1,
1083 center2=center2, rad2=rad2, wrap=wrap)
1086def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1087 height=None, too_d=None, wrap=False, # was=True
1088 LatLon=LatLon, **LatLon_kwds):
1089 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1090 # above, separated to allow callers to embellish any exceptions
1092 def _dest3(bearing, h):
1093 a, b = _destination2(a1, b1, r1, bearing)
1094 return _LL3Tuple(degrees90(a), degrees180(b), h,
1095 intersections2, LatLon, LatLon_kwds)
1097 a1, b1 = c1.philam
1098 a2, b2 = c2.philam
1099 if wrap:
1100 a2, b2 = _Wrap.philam(a2, b2)
1102 r1, r2, f = _rads3(rad1, rad2, radius)
1103 if f: # swapped radii, swap centers
1104 a1, a2 = a2, a1 # PYCHOK swap!
1105 b1, b2 = b2, b1 # PYCHOK swap!
1107 db, b2 = unrollPI(b1, b2, wrap=wrap)
1108 d = vincentys_(a2, a1, db) # radians
1109 if d < max(r1 - r2, EPS):
1110 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1112 r = eps if radius is None else (m2radians(
1113 eps, radius=radius) if eps else _0_0)
1114 if r < _0_0:
1115 raise _ValueError(eps=r)
1117 x = fsumf_(r1, r2, -d) # overlap
1118 if x > max(r, EPS):
1119 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1120 x = sd * sr1
1121 if isnear0(x):
1122 raise _ValueError(_invalid_)
1123 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1125 elif x < r: # PYCHOK no cover
1126 t = (d * radius) if too_d is None else too_d
1127 raise IntersectionError(_too_(_Fmt.distant(t)))
1129 if height is None: # "radical height"
1130 f = _radical2(d, r1, r2).ratio
1131 h = Height(favg(c1.height, c2.height, f=f))
1132 else:
1133 h = Height(height)
1135 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1136 if x < EPS4: # externally ...
1137 r = _dest3(b, h)
1138 elif x > _PI_EPS4: # internally ...
1139 r = _dest3(b + PI, h)
1140 else:
1141 return _dest3(b + x, h), _dest3(b - x, h)
1142 return r, r # ... abutting circles
1145@deprecated_function
1146def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1147 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1148 '''
1149 return ispolar(points, wrap=wrap)
1152def _LL3Tuple(lat, lon, height, where, LatLon, LatLon_kwds):
1153 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1154 '''
1155 n = where.__name__
1156 if LatLon is None:
1157 r = LatLon3Tuple(lat, lon, height, name=n)
1158 else:
1159 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1160 r = LatLon(lat, lon, **kwds)
1161 return r
1164def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):
1165 '''Compute the I{geographic} mean of several points.
1167 @arg points: Points to be averaged (L{LatLon}[]).
1168 @kwarg height: Optional height at mean point, overriding the mean height
1169 (C{meter}).
1170 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}} (C{bool}).
1171 @kwarg LatLon: Optional class to return the mean point (L{LatLon}) or C{None}.
1172 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments,
1173 ignored if C{B{LatLon} is None}.
1175 @return: The geographic mean and height (B{C{LatLon}}) or if C{B{LatLon}
1176 is None}, a L{LatLon3Tuple}C{(lat, lon, height)}.
1178 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1180 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1181 '''
1182 def _N_vs(ps, w):
1183 Ps = _T00.PointsIter(ps, wrap=w)
1184 for p in Ps.iterate(closed=False):
1185 yield p._N_vector
1187 m = _MODS.nvectorBase
1188 # geographic, vectorial mean
1189 n = m.sumOf(_N_vs(points, wrap), h=height, Vector=m.NvectorBase)
1190 lat, lon, h = n.latlonheight
1191 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1194@deprecated_function
1195def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1196 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1198 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1199 point (L{LatLon}) on the polygon and the C{distance}
1200 between the C{closest} and the given B{C{point}}. The
1201 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1202 lon)} if C{B{LatLon} is None} ...
1203 '''
1204 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1205 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1206 ll = LatLon2Tuple(ll.lat, ll.lon)
1207 return ll, d
1210def nearestOn3(point, points, closed=False, radius=R_M, wrap=False, adjust=True,
1211 limit=9, **LatLon_and_kwds):
1212 '''Locate the point on a path or polygon closest to a reference point.
1214 Distances are I{approximated} using function L{equirectangular4
1215 <pygeodesy.equirectangular4>}, subject to the supplied B{C{options}}.
1217 @arg point: The reference point (L{LatLon}).
1218 @arg points: The path or polygon points (L{LatLon}[]).
1219 @kwarg closed: Optionally, close the polygon (C{bool}).
1220 @kwarg radius: Mean earth radius (C{meter}).
1221 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1222 B{C{points}} (C{bool}).
1223 @kwarg adjust: See function L{equirectangular4<pygeodesy.equirectangular4>}
1224 (C{bool}).
1225 @kwarg limit: See function L{equirectangular4<pygeodesy.equirectangular4>}
1226 (C{degrees}), default C{9 degrees} is about C{1,000 Km} (for
1227 (mean spherical earth radius L{R_KM}).
1228 @kwarg LatLon_and_kwds: Optional class C{B{LatLon}=L{LatLon}} to return the
1229 closest point and optional, additional C{B{LatLon}} keyword
1230 arguments or specify C{B{LatLon}=None}.
1232 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1233 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1234 lon, height)} if C{B{LatLon} is None}. The C{distance} is
1235 the L{equirectangular4<pygeodesy.equirectangular4>} distance
1236 between the C{closest} and the given B{C{point}} converted to
1237 C{meter}, same units as B{C{radius}}. The C{angle} from the
1238 given B{C{point}} to the C{closest} is in compass C{degrees360},
1239 like function L{compassAngle<pygeodesy.compassAngle>}. The
1240 C{height} is the (interpolated) height at the C{closest} point.
1242 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1243 see function L{equirectangular4<pygeodesy.equirectangular4>}.
1245 @raise PointsError: Insufficient number of B{C{points}}.
1247 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1249 @raise ValueError: Invalid B{C{radius}}.
1251 @see: Functions L{equirectangular4<pygeodesy.equirectangular4>} and
1252 L{nearestOn5<pygeodesy.nearestOn5>}.
1253 '''
1254 t = _nearestOn5(point, points, closed=closed, wrap=wrap,
1255 adjust=adjust, limit=limit)
1256 d = degrees2m(t.distance, radius=radius)
1257 h = t.height
1258 n = nearestOn3.__name__
1260 LL, kwds = _xkwds_pop2(LatLon_and_kwds, LatLon=LatLon)
1261 r = LatLon3Tuple(t.lat, t.lon, h, name=n) if LL is None else \
1262 LL(t.lat, t.lon, **_xkwds(kwds, height=h, name=n))
1263 return NearestOn3Tuple(r, d, t.angle, name=n)
1266def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1267 '''Compute the perimeter of a (spherical) polygon or composite
1268 (with great circle arcs joining the points).
1270 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1271 or L{BooleanGH}).
1272 @kwarg closed: Optionally, close the polygon (C{bool}).
1273 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1274 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1275 B{C{points}} (C{bool}).
1277 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1278 or C{radians} if C{B{radius} is None}).
1280 @raise PointsError: Insufficient number of B{C{points}}.
1282 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1284 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1285 C{B{points}} a composite.
1287 @note: Distances are based on function L{vincentys_<pygeodesy.vincentys_>}.
1289 @see: Functions L{perimeterOf<pygeodesy.perimeterOf>},
1290 L{sphericalNvector.perimeterOf} and L{ellipsoidalKarney.perimeterOf}.
1291 '''
1292 def _rads(ps, c, w): # angular edge lengths in radians
1293 Ps = _T00.PointsIter(ps, loop=1, wrap=w)
1294 a1, b1 = Ps[0].philam
1295 for p in Ps.iterate(closed=c):
1296 a2, b2 = p.philam
1297 db, b2 = unrollPI(b1, b2, wrap=w and not (c and Ps.looped))
1298 yield vincentys_(a2, a1, db)
1299 a1, b1 = a2, b2
1301 if _MODS.booleans.isBoolean(points):
1302 if not closed:
1303 raise _ValueError(closed=closed, points=_composite_)
1304 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1305 else:
1306 r = fsum(_rads(points, closed, wrap))
1307 return _radians2m(r, radius)
1310def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1311 excess=excessAbc_,
1312 wrap=False):
1313 '''Compute the angles, sides, and area of a (spherical) triangle.
1315 @arg latA: First corner latitude (C{degrees}).
1316 @arg lonA: First corner longitude (C{degrees}).
1317 @arg latB: Second corner latitude (C{degrees}).
1318 @arg lonB: Second corner longitude (C{degrees}).
1319 @arg latC: Third corner latitude (C{degrees}).
1320 @arg lonC: Third corner longitude (C{degrees}).
1321 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1322 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1323 or C{None}.
1324 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1325 L{excessGirard_} or L{excessLHuilier_}).
1326 @kwarg wrap: If C{True}, wrap and L{unroll180<pygeodesy.unroll180>}
1327 longitudes (C{bool}).
1329 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1330 spherical angles C{A}, C{B} and C{C}, angular sides
1331 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1332 in I{square} C{meter} or same units as B{C{radius}}
1333 I{squared} or if C{B{radius}=0} or C{None}, a
1334 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in
1335 C{radians} with the I{spherical excess} C{E} as the
1336 C{unit area} in C{radians}.
1337 '''
1338 t = triangle8_(Phid(latA=latA), Lamd(lonA=lonA),
1339 Phid(latB=latB), Lamd(lonB=lonB),
1340 Phid(latC=latC), Lamd(lonC=lonC),
1341 excess=excess, wrap=wrap)
1342 return _t7Tuple(t, radius)
1345def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_,
1346 wrap=False):
1347 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1348 excess} of a (spherical) triangle.
1350 @arg phiA: First corner latitude (C{radians}).
1351 @arg lamA: First corner longitude (C{radians}).
1352 @arg phiB: Second corner latitude (C{radians}).
1353 @arg lamB: Second corner longitude (C{radians}).
1354 @arg phiC: Third corner latitude (C{radians}).
1355 @arg lamC: Third corner longitude (C{radians}).
1356 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1357 L{excessGirard_} or L{excessLHuilier_}).
1358 @kwarg wrap: If C{True}, L{unrollPI<pygeodesy.unrollPI>} the
1359 longitudinal deltas (C{bool}).
1361 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1362 spherical angles C{A}, C{B} and C{C}, angular sides
1363 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1364 I{spherical excess} C{E}, all in C{radians}.
1365 '''
1366 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1367 d, _ = unrollPI(lamB, lamC, wrap=w)
1368 a = vincentys_(phiC, phiB, d)
1369 return a, (phiB, lamB, phiC, lamC, phiA, lamA) # rotate A, B, C
1371 def _A_r(a, sa, ca, sb, cb, sc, cc):
1372 s = sb * sc
1373 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1374 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2_'s
1376 # notation: side C{a} is oposite to corner C{A}, etc.
1377 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1378 b, r = _a_r(wrap, *r)
1379 c, _ = _a_r(wrap, *r)
1381 A, r = _A_r(a, *sincos2_(a, b, c))
1382 B, r = _A_r(b, *r)
1383 C, _ = _A_r(c, *r)
1385 D = fsumf_(PI2, -a, -b, -c) # deficit aka defect
1386 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else (
1387 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else
1388 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b))))
1390 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1393def _t7Tuple(t, radius):
1394 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1395 '''
1396 if radius: # not in (None, _0_0)
1397 r = radius if _isRadius(radius) else \
1398 _ellipsoidal_datum(radius).ellipsoid.Rmean
1399 A, B, C = map1(degrees, t.A, t.B, t.C)
1400 t = Triangle7Tuple(A, (r * t.a),
1401 B, (r * t.b),
1402 C, (r * t.c), t.E * r**2)
1403 return t
1406__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1407 areaOf, # functions
1408 intersecant2, intersection, intersections2, ispolar,
1409 isPoleEnclosedBy, # DEPRECATED, use ispolar
1410 meanOf,
1411 nearestOn2, nearestOn3,
1412 perimeterOf,
1413 sumOf, # XXX == vector3d.sumOf
1414 triangle7, triangle8_)
1416# **) MIT License
1417#
1418# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1419#
1420# Permission is hereby granted, free of charge, to any person obtaining a
1421# copy of this software and associated documentation files (the "Software"),
1422# to deal in the Software without restriction, including without limitation
1423# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1424# and/or sell copies of the Software, and to permit persons to whom the
1425# Software is furnished to do so, subject to the following conditions:
1426#
1427# The above copyright notice and this permission notice shall be included
1428# in all copies or substantial portions of the Software.
1429#
1430# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1431# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1432# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1433# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1434# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1435# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1436# OTHER DEALINGS IN THE SOFTWARE.