Coverage for pygeodesy/sphericalTrigonometry.py: 94%
390 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-10-04 14:05 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-10-04 14:05 -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-2016} published under the same MIT Licence**, see
14U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}.
15'''
16# make sure int/int division yields float quotient, see .basics
17from __future__ import division as _; del _ # PYCHOK semicolon
19from pygeodesy.basics import copysign0, isscalar, 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 _ValueError, IntersectionError, _xError, \
26 _xkwds, _xkwds_get, _xkwds_pop
27from pygeodesy.fmath import favg, fdot, fmean, hypot
28from pygeodesy.fsums import Fsum, fsum, fsumf_
29from pygeodesy.formy import antipode_, bearing_, _bearingTo2, excessAbc_, \
30 excessGirard_, excessLHuilier_, opposing_, _radical2, \
31 vincentys_
32from pygeodesy.interns import _1_, _2_, _coincident_, _composite_, _colinear_, \
33 _concentric_, _convex_, _end_, _infinite_, _invalid_,\
34 _line_, _near_, _not_, _null_, _parallel_, _point_, \
35 _SPACE_, _too_
36from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER
37# from pygeodesy.named import notImplemented # from .points
38# from pygeodesy.nvectorBase import NvectorBase, sumOf # _MODE
39from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, NearestOn3Tuple, \
40 Triangle7Tuple, Triangle8Tuple
41from pygeodesy.points import ispolar, nearestOn5 as _nearestOn5, \
42 Fmt as _Fmt, notImplemented # XXX shadowed
43from pygeodesy.props import deprecated_function, deprecated_method
44from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, _intersecant2, \
45 LatLonSphericalBase, _rads3, _r2m, _trilaterate5
46# from pygeodesy.streprs import Fmt as _Fmt # from .points XXX shadowed
47from pygeodesy.units import Bearing_, Height, Lam_, Phi_, Radius_, Scalar
48from pygeodesy.utily import acos1, asin1, atan1d, atan2d, degrees90, degrees180, \
49 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \
50 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI
51from pygeodesy.vector3d import sumOf, Vector3d
53from math import asin, atan2, cos, degrees, fabs, radians, sin
55__all__ = _ALL_LAZY.sphericalTrigonometry
56__version__ = '23.09.29'
58_PI_EPS4 = PI - EPS4
59if _PI_EPS4 >= PI:
60 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4)
63class Cartesian(CartesianSphericalBase):
64 '''Extended to convert geocentric, L{Cartesian} points to
65 spherical, geodetic L{LatLon}.
66 '''
68 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon
69 '''Convert this cartesian point to a C{spherical} geodetic point.
71 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
72 arguments. Use C{B{LatLon}=...} to override
73 this L{LatLon} class or specify C{B{LatLon}=None}.
75 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is C{None},
76 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
77 with C{C} and C{M} if available.
79 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
80 '''
81 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
82 return CartesianSphericalBase.toLatLon(self, **kwds)
85class LatLon(LatLonSphericalBase):
86 '''New point on spherical model earth model.
88 @example:
90 >>> p = LatLon(52.205, 0.119) # height=0
91 '''
93 def _ab1_ab2_db5(self, other, wrap):
94 '''(INTERNAL) Helper for several methods.
95 '''
96 a1, b1 = self.philam
97 a2, b2 = self.others(other, up=2).philam
98 if wrap:
99 a2, b2 = _Wrap.philam(a2, b2)
100 db, b2 = unrollPI(b1, b2, wrap=wrap)
101 else: # unrollPI shortcut
102 db = b2 - b1
103 return a1, b1, a2, b2, db
105 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
106 '''Compute the (signed) distance from the start to the closest
107 point on the great circle line defined by a start and an
108 end point.
110 That is, if a perpendicular is drawn from this point to the
111 great circle line, the along-track distance is the distance
112 from the start point to the point where the perpendicular
113 crosses the line.
115 @arg start: Start point of the great circle line (L{LatLon}).
116 @arg end: End point of the great circle line (L{LatLon}).
117 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
118 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
119 the B{C{start}} and B{C{end}} point (C{bool}).
121 @return: Distance along the great circle line (C{radians}
122 if C{B{radius} is None} or C{meter}, same units
123 as B{C{radius}}), positive if I{after} the
124 B{C{start}} toward the B{C{end}} point of the
125 line, I{negative} if before or C{0} if at the
126 B{C{start}} point.
128 @raise TypeError: Invalid B{C{start}} or B{C{end}} point.
130 @raise ValueError: Invalid B{C{radius}}.
132 @example:
134 >>> p = LatLon(53.2611, -0.7972)
136 >>> s = LatLon(53.3206, -1.7297)
137 >>> e = LatLon(53.1887, 0.1334)
138 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58
139 '''
140 r, x, b = self._a_x_b3(start, end, radius, wrap)
141 cx = cos(x)
142 return _0_0 if isnear0(cx) else \
143 _r2m(copysign0(acos1(cos(r) / cx), cos(b)), radius)
145 def _a_x_b3(self, start, end, radius, wrap):
146 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo.
147 '''
148 s = self.others(start=start)
149 e = self.others(end=end)
150 s, e, w = _unrollon3(self, s, e, wrap)
152 r = Radius_(radius)
153 r = s.distanceTo(self, r, wrap=w) / r
155 b = radians(s.initialBearingTo(self, wrap=w)
156 - s.initialBearingTo(e, wrap=w))
157 x = asin(sin(r) * sin(b))
158 return r, x, -b
160 @deprecated_method
161 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover
162 '''DEPRECATED, use method L{initialBearingTo}.
163 '''
164 return self.initialBearingTo(other, wrap=wrap, raiser=raiser)
166 def crossingParallels(self, other, lat, wrap=False):
167 '''Return the pair of meridians at which a great circle defined
168 by this and an other point crosses the given latitude.
170 @arg other: The other point defining great circle (L{LatLon}).
171 @arg lat: Latitude at the crossing (C{degrees}).
172 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
173 B{C{other}} point (C{bool}).
175 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
176 C{None} if the great circle doesn't reach B{C{lat}}.
177 '''
178 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
179 sa, ca, sa1, ca1, \
180 sa2, ca2, sdb, cdb = sincos2_(radians(lat), a1, a2, db)
181 sa1 *= ca2 * ca
183 x = sa1 * sdb
184 y = sa1 * cdb - ca1 * sa2 * ca
185 z = ca1 * sdb * ca2 * sa
187 h = hypot(x, y)
188 if h < EPS or fabs(z) > h: # PYCHOK no cover
189 return None # great circle doesn't reach latitude
191 m = atan2(-y, x) + b1 # longitude at max latitude
192 d = acos1(z / h) # delta longitude to intersections
193 return degrees180(m - d), degrees180(m + d)
195 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
196 '''Compute the (signed) distance from this point to
197 the great circle defined by a start and an end point.
199 @arg start: Start point of the great circle line (L{LatLon}).
200 @arg end: End point of the great circle line (L{LatLon}).
201 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
202 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
203 the B{C{start}} and B{C{end}} point (C{bool}).
205 @return: Distance to the great circle (C{radians} if
206 B{C{radius}} or C{meter}, same units as
207 B{C{radius}}), I{negative} if to the left or
208 I{positive} if to the right of the line.
210 @raise TypeError: If B{C{start}} or B{C{end}} is not L{LatLon}.
212 @raise ValueError: Invalid B{C{radius}}.
214 @example:
216 >>> p = LatLon(53.2611, -0.7972)
218 >>> s = LatLon(53.3206, -1.7297)
219 >>> e = LatLon(53.1887, 0.1334)
220 >>> d = p.crossTrackDistanceTo(s, e) # -307.5
221 '''
222 _, x, _ = self._a_x_b3(start, end, radius, wrap)
223 return _r2m(x, radius)
225 def destination(self, distance, bearing, radius=R_M, height=None):
226 '''Locate the destination from this point after having
227 travelled the given distance on the given initial bearing.
229 @arg distance: Distance travelled (C{meter}, same units as
230 B{C{radius}}).
231 @arg bearing: Bearing from this point (compass C{degrees360}).
232 @kwarg radius: Mean earth radius (C{meter}).
233 @kwarg height: Optional height at destination (C{meter}, same
234 units a B{C{radius}}).
236 @return: Destination point (L{LatLon}).
238 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
239 B{C{radius}} or B{C{height}}.
241 @example:
243 >>> p1 = LatLon(51.4778, -0.0015)
244 >>> p2 = p1.destination(7794, 300.7)
245 >>> p2.toStr() # '51.5135°N, 000.0983°W'
246 '''
247 a, b = self.philam
248 r, t = _angular(distance, radius, low=None), Bearing_(bearing)
250 a, b = _destination2(a, b, r, t)
251 h = self._heigHt(height)
252 return self.classof(degrees90(a), degrees180(b), height=h)
254 def distanceTo(self, other, radius=R_M, wrap=False):
255 '''Compute the (angular) distance from this to an other point.
257 @arg other: The other point (L{LatLon}).
258 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
259 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
260 the B{C{other}} point (C{bool}).
262 @return: Distance between this and the B{C{other}} point
263 (C{meter}, same units as B{C{radius}} or
264 C{radians} if B{C{radius}} is C{None}).
266 @raise TypeError: The B{C{other}} point is not L{LatLon}.
268 @raise ValueError: Invalid B{C{radius}}.
270 @example:
272 >>> p1 = LatLon(52.205, 0.119)
273 >>> p2 = LatLon(48.857, 2.351);
274 >>> d = p1.distanceTo(p2) # 404300
275 '''
276 a1, _, a2, _, db = self._ab1_ab2_db5(other, wrap)
277 return _r2m(vincentys_(a2, a1, db), radius)
279# @Property_RO
280# def Ecef(self):
281# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
282# '''
283# return _MODS.ecef.EcefKarney
285 def greatCircle(self, bearing, Vector=Vector3d, **Vector_kwds):
286 '''Compute the vector normal to great circle obtained by heading
287 on the given initial bearing from this point.
289 Direction of vector is such that initial bearing vector
290 b = c × n, where n is an n-vector representing this point.
292 @arg bearing: Bearing from this point (compass C{degrees360}).
293 @kwarg Vector: Vector class to return the great circle,
294 overriding the default L{Vector3d}.
295 @kwarg Vector_kwds: Optional, additional keyword argunents
296 for B{C{Vector}}.
298 @return: Vector representing great circle (C{Vector}).
300 @raise ValueError: Invalid B{C{bearing}}.
302 @example:
304 >>> p = LatLon(53.3206, -1.7297)
305 >>> g = p.greatCircle(96.0)
306 >>> g.toStr() # (-0.794, 0.129, 0.594)
307 '''
308 a, b = self.philam
309 sa, ca, sb, cb, st, ct = sincos2_(a, b, Bearing_(bearing))
311 return Vector(sb * ct - cb * sa * st,
312 -cb * ct - sb * sa * st,
313 ca * st, **Vector_kwds) # XXX .unit()?
315 def initialBearingTo(self, other, wrap=False, raiser=False):
316 '''Compute the initial bearing (forward azimuth) from this
317 to an other point.
319 @arg other: The other point (spherical L{LatLon}).
320 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
321 the B{C{other}} point (C{bool}).
322 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}),
323 use C{B{raiser}=True} for behavior like
324 C{sphericalNvector.LatLon.initialBearingTo}.
326 @return: Initial bearing (compass C{degrees360}).
328 @raise CrossError: If this and the B{C{other}} point coincide,
329 provided both B{C{raiser}} is C{True} and
330 L{pygeodesy.crosserrors} is C{True}.
332 @raise TypeError: The B{C{other}} point is not L{LatLon}.
334 @example:
336 >>> p1 = LatLon(52.205, 0.119)
337 >>> p2 = LatLon(48.857, 2.351)
338 >>> b = p1.initialBearingTo(p2) # 156.2
339 '''
340 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
341 # XXX behavior like sphericalNvector.LatLon.initialBearingTo
342 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(db)) < EPS:
343 raise CrossError(_point_, self, other=other, wrap=wrap, txt=_coincident_)
345 return degrees(bearing_(a1, b1, a2, b2, final=False))
347 def intermediateTo(self, other, fraction, height=None, wrap=False):
348 '''Locate the point at given fraction between (or along) this
349 and an other point.
351 @arg other: The other point (L{LatLon}).
352 @arg fraction: Fraction between both points (C{scalar},
353 0.0 at this and 1.0 at the other point).
354 @kwarg height: Optional height, overriding the intermediate
355 height (C{meter}).
356 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
357 B{C{other}} point (C{bool}).
359 @return: Intermediate point (L{LatLon}).
361 @raise TypeError: The B{C{other}} point is not L{LatLon}.
363 @raise ValueError: Invalid B{C{fraction}} or B{C{height}}.
365 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
367 @example:
369 >>> p1 = LatLon(52.205, 0.119)
370 >>> p2 = LatLon(48.857, 2.351)
371 >>> p = p1.intermediateTo(p2, 0.25) # 51.3721°N, 000.7073°E
372 '''
373 p = self
374 f = Scalar(fraction=fraction)
375 if not isnear0(f):
376 p = p.others(other)
377 if wrap:
378 p = _Wrap.point(p)
379 if not isnear1(f): # and not near0
380 a1, b1 = self.philam
381 a2, b2 = p.philam
382 db, b2 = unrollPI(b1, b2, wrap=wrap)
383 r = vincentys_(a2, a1, db)
384 sr = sin(r)
385 if isnon0(sr):
386 sa1, ca1, sa2, ca2, \
387 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2)
389 t = f * r
390 a = sin(r - t) # / sr superflous
391 b = sin( t) # / sr superflous
393 x = a * ca1 * cb1 + b * ca2 * cb2
394 y = a * ca1 * sb1 + b * ca2 * sb2
395 z = a * sa1 + b * sa2
397 a = atan1d(z, hypot(x, y))
398 b = atan2d(y, x)
400 else: # PYCHOK no cover
401 a = degrees90( favg(a1, a2, f=f)) # coincident
402 b = degrees180(favg(b1, b2, f=f))
404 h = self._havg(other, f=f, h=height)
405 p = self.classof(a, b, height=h)
406 return p
408 def intersection(self, end1, other, end2, height=None, wrap=False):
409 '''Compute the intersection point of two lines, each defined by
410 two points or a start point and bearing from North.
412 @arg end1: End point of this line (L{LatLon}) or the initial
413 bearing at this point (compass C{degrees360}).
414 @arg other: Start point of the other line (L{LatLon}).
415 @arg end2: End point of the other line (L{LatLon}) or the
416 initial bearing at the B{C{other}} point (compass
417 C{degrees360}).
418 @kwarg height: Optional height for intersection point,
419 overriding the mean height (C{meter}).
420 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
421 B{C{start2}} and both B{C{end*}} points (C{bool}).
423 @return: The intersection point (L{LatLon}). An alternate
424 intersection point might be the L{antipode} to
425 the returned result.
427 @raise IntersectionError: Ambiguous or infinite intersection
428 or colinear, parallel or otherwise
429 non-intersecting lines.
431 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}}
432 or B{C{end2}} not C{scalar} nor L{LatLon}.
434 @raise ValueError: Invalid B{C{height}} or C{null} line.
436 @example:
438 >>> p = LatLon(51.8853, 0.2545)
439 >>> s = LatLon(49.0034, 2.5735)
440 >>> i = p.intersection(108.547, s, 32.435) # '50.9078°N, 004.5084°E'
441 '''
442 try:
443 s2 = self.others(other)
444 return _intersect(self, end1, s2, end2, height=height, wrap=wrap,
445 LatLon=self.classof)
446 except (TypeError, ValueError) as x:
447 raise _xError(x, start1=self, end1=end1,
448 other=other, end2=end2, wrap=wrap)
450 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0,
451 height=None, wrap=True):
452 '''Compute the intersection points of two circles, each defined
453 by a center point and radius.
455 @arg rad1: Radius of the this circle (C{meter} or C{radians},
456 see B{C{radius}}).
457 @arg other: Center point of the other circle (L{LatLon}).
458 @arg rad2: Radius of the other circle (C{meter} or C{radians},
459 see B{C{radius}}).
460 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
461 B{C{rad2}} and B{C{eps}} are given in C{radians}).
462 @kwarg eps: Required overlap (C{meter} or C{radians}, see
463 B{C{radius}}).
464 @kwarg height: Optional height for the intersection points (C{meter},
465 conventionally) or C{None} for the I{"radical height"}
466 at the I{radical line} between both centers.
467 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
468 B{C{other}} point (C{bool}).
470 @return: 2-Tuple of the intersection points, each a L{LatLon}
471 instance. For abutting circles, both intersection
472 points are the same instance, aka the I{radical center}.
474 @raise IntersectionError: Concentric, antipodal, invalid or
475 non-intersecting circles.
477 @raise TypeError: If B{C{other}} is not L{LatLon}.
479 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
480 B{C{eps}} or B{C{height}}.
481 '''
482 try:
483 c2 = self.others(other)
484 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps,
485 height=height, wrap=wrap,
486 LatLon=self.classof)
487 except (TypeError, ValueError) as x:
488 raise _xError(x, center=self, rad1=rad1,
489 other=other, rad2=rad2, wrap=wrap)
491 @deprecated_method
492 def isEnclosedBy(self, points): # PYCHOK no cover
493 '''DEPRECATED, use method C{isenclosedBy}.'''
494 return self.isenclosedBy(points)
496 def isenclosedBy(self, points, wrap=False):
497 '''Check whether a (convex) polygon or composite encloses this point.
499 @arg points: The polygon points or composite (L{LatLon}[],
500 L{BooleanFHP} or L{BooleanGH}).
501 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
502 B{C{points}} (C{bool}).
504 @return: C{True} if this point is inside the polygon or
505 composite, C{False} otherwise.
507 @raise PointsError: Insufficient number of B{C{points}}.
509 @raise TypeError: Some B{C{points}} are not L{LatLon}.
511 @raise ValueError: Invalid B{C{points}}, non-convex polygon.
513 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
514 and L{pygeodesy.ispolar} especially if the B{C{points}} may
515 enclose a pole or wrap around the earth I{longitudinally}.
516 '''
517 if _MODS.booleans.isBoolean(points):
518 return points._encloses(self.lat, self.lon, wrap=wrap)
520 Ps = self.PointsIter(points, loop=2, dedup=True, wrap=wrap)
521 n0 = self._N_vector
523 v2 = Ps[0]._N_vector
524 p1 = Ps[1]
525 v1 = p1._N_vector
526 # check whether this point on same side of all
527 # polygon edges (to the left or right depending
528 # on the anti-/clockwise polygon direction)
529 gc1 = v2.cross(v1)
530 t0 = gc1.angleTo(n0) > PI_2
531 s0 = None
532 # get great-circle vector for each edge
533 for i, p2 in Ps.enumerate(closed=True):
534 if wrap and not Ps.looped:
535 p2 = _unrollon(p1, p2)
536 p1 = p2
537 v2 = p2._N_vector
538 gc = v1.cross(v2)
539 t = gc.angleTo(n0) > PI_2
540 if t != t0: # different sides of edge i
541 return False # outside
543 # check for convex polygon: angle between
544 # gc vectors, signed by direction of n0
545 # (otherwise the test above is not reliable)
546 s = signOf(gc1.angleTo(gc, vSign=n0))
547 if s != s0:
548 if s0 is None:
549 s0 = s
550 else:
551 t = _Fmt.SQUARE(points=i)
552 raise _ValueError(t, p2, wrap=wrap, txt=_not_(_convex_))
553 gc1, v1 = gc, v2
555 return True # inside
557 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
558 '''Find the midpoint between this and an other point.
560 @arg other: The other point (L{LatLon}).
561 @kwarg height: Optional height for midpoint, overriding
562 the mean height (C{meter}).
563 @kwarg fraction: Midpoint location from this point (C{scalar}),
564 may be negative or greater than 1.0.
565 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
566 B{C{other}} point (C{bool}).
568 @return: Midpoint (L{LatLon}).
570 @raise TypeError: The B{C{other}} point is not L{LatLon}.
572 @raise ValueError: Invalid B{C{height}}.
574 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
576 @example:
578 >>> p1 = LatLon(52.205, 0.119)
579 >>> p2 = LatLon(48.857, 2.351)
580 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
581 '''
582 if fraction is _0_5:
583 # see <https://MathForum.org/library/drmath/view/51822.html>
584 a1, b, a2, _, db = self._ab1_ab2_db5(other, wrap)
585 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db)
587 x = ca2 * cdb + ca1
588 y = ca2 * sdb
590 a = atan1d(sa1 + sa2, hypot(x, y))
591 b = degrees180(b + atan2(y, x))
593 h = self._havg(other, h=height)
594 r = self.classof(a, b, height=h)
595 else:
596 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
597 return r
599 def nearestOn(self, point1, point2, radius=R_M, **wrap_adjust_limit):
600 '''Locate the point between two points closest to this point.
602 Distances are approximated by function L{pygeodesy.equirectangular_},
603 subject to the supplied B{C{options}}.
605 @arg point1: Start point (L{LatLon}).
606 @arg point2: End point (L{LatLon}).
607 @kwarg radius: Mean earth radius (C{meter}).
608 @kwarg wrap_adjust_limit: Optional keyword arguments for functions
609 L{sphericalTrigonometry.nearestOn3} and
610 L{pygeodesy.equirectangular_},
612 @return: Closest point on the great circle line (L{LatLon}).
614 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
615 see function L{pygeodesy.equirectangular_}.
617 @raise NotImplementedError: Keyword argument C{B{within}=False}
618 is not (yet) supported.
620 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
622 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
624 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}
625 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
626 '''
627 # remove kwarg B{C{within}} if present
628 w = _xkwds_pop(wrap_adjust_limit, within=True)
629 if not w:
630 notImplemented(self, within=w)
632# # UNTESTED - handle C{B{within}=False} and C{B{within}=True}
633# wrap = _xkwds_get(options, wrap=False)
634# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap)
635# if fabs(a) < EPS or (within and a < EPS):
636# return point1
637# d = point1.distanceTo(point2, radius=radius, wrap=wrap)
638# if isnear0(d):
639# return point1 # or point2
640# elif fabs(d - a) < EPS or (a + EPS) > d:
641# return point2
642# f = a / d
643# if within:
644# if f > EPS1:
645# return point2
646# elif f < EPS:
647# return point1
648# return point1.intermediateTo(point2, f, wrap=wrap)
650 # without kwarg B{C{within}}, use backward compatible .nearestOn3
651 return self.nearestOn3([point1, point2], closed=False, radius=radius,
652 **wrap_adjust_limit)[0]
654 @deprecated_method
655 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover
656 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
658 @return: ... 2-Tuple C{(closest, distance)} of the closest
659 point (L{LatLon}) on the polygon and the distance
660 to that point from this point in C{meter}, same
661 units of B{C{radius}}.
662 '''
663 r = self.nearestOn3(points, closed=closed, radius=radius, **options)
664 return r.closest, r.distance
666 def nearestOn3(self, points, closed=False, radius=R_M, **wrap_adjust_limit):
667 '''Locate the point on a polygon closest to this point.
669 Distances are approximated by function L{pygeodesy.equirectangular_},
670 subject to the supplied B{C{options}}.
672 @arg points: The polygon points (L{LatLon}[]).
673 @kwarg closed: Optionally, close the polygon (C{bool}).
674 @kwarg radius: Mean earth radius (C{meter}).
675 @kwarg wrap_adjust_limit: Optional keyword arguments for function
676 L{sphericalTrigonometry.nearestOn3} and
677 L{pygeodesy.equirectangular_},
679 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
680 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_}
681 C{distance} between this and the C{closest} point converted to
682 C{meter}, same units as B{C{radius}}. The C{angle} from this
683 to the C{closest} point is in compass C{degrees360}, like
684 function L{pygeodesy.compassAngle}.
686 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
687 see function L{pygeodesy.equirectangular_}.
689 @raise PointsError: Insufficient number of B{C{points}}.
691 @raise TypeError: Some B{C{points}} are not C{LatLon}.
693 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
695 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_}
696 and L{pygeodesy.nearestOn5}.
697 '''
698 return nearestOn3(self, points, closed=closed, radius=radius,
699 LatLon=self.classof, **wrap_adjust_limit)
701 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
702 '''Convert this point to C{Karney}-based cartesian (ECEF)
703 coordinates.
705 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}}
706 and other keyword arguments, ignored
707 if C{B{Cartesian} is None}. Use
708 C{B{Cartesian}=...} to override
709 this L{Cartesian} class or specify
710 C{B{Cartesian}=None}.
712 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
713 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
714 C, M, datum)} with C{C} and C{M} if available.
716 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
717 '''
718 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
719 return LatLonSphericalBase.toCartesian(self, **kwds)
721 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
722 '''Compute the angles, sides and area of a spherical triangle.
724 @arg otherB: Second triangle point (C{LatLon}).
725 @arg otherC: Third triangle point (C{LatLon}).
726 @kwarg radius: Mean earth radius, ellipsoid or datum
727 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},
728 L{Datum} or L{a_f2Tuple}) or C{None}.
729 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
730 B{C{otherB}} and B{C{otherC}} points (C{bool}).
732 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if
733 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A,
734 a, B, b, C, c, D, E)}.
736 @see: Function L{triangle7} and U{Spherical trigonometry
737 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
738 '''
739 B = self.others(otherB=otherB)
740 C = self.others(otherC=otherC)
741 B, C, _ = _unrollon3(self, B, C, wrap)
743 r = self.philam + B.philam + C.philam
744 t = triangle8_(*r, wrap=wrap)
745 return self._xnamed(_t7Tuple(t, radius))
747 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
748 area=True, eps=EPS1, radius=R_M, wrap=False):
749 '''Trilaterate three points by I{area overlap} or I{perimeter
750 intersection} of three corresponding circles.
752 @arg distance1: Distance to this point (C{meter}, same units
753 as B{C{radius}}).
754 @arg point2: Second center point (C{LatLon}).
755 @arg distance2: Distance to point2 (C{meter}, same units as
756 B{C{radius}}).
757 @arg point3: Third center point (C{LatLon}).
758 @arg distance3: Distance to point3 (C{meter}, same units as
759 B{C{radius}}).
760 @kwarg area: If C{True} compute the area overlap, otherwise the
761 perimeter intersection of the circles (C{bool}).
762 @kwarg eps: The required I{minimal overlap} for C{B{area}=True}
763 or the I{intersection margin} for C{B{area}=False}
764 (C{meter}, same units as B{C{radius}}).
765 @kwarg radius: Mean earth radius (C{meter}, conventionally).
766 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
767 B{C{point2}} and B{C{point3}} (C{bool}).
769 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
770 with C{min} and C{max} in C{meter}, same units as B{C{eps}},
771 the corresponding trilaterated points C{minPoint} and
772 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number
773 of trilatered points found for the given B{C{eps}}.
775 If only a single trilaterated point is found, C{min I{is}
776 max}, C{minPoint I{is} maxPoint} and C{n = 1}.
778 For C{B{area}=True}, C{min} and C{max} are the smallest
779 respectively largest I{radial} overlap found.
781 For C{B{area}=False}, C{min} and C{max} represent the
782 nearest respectively farthest intersection margin.
784 If C{B{area}=True} and all 3 circles are concentric, C{n =
785 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}}
786 with the smallest B{C{distance#}} C{min} and C{max} the
787 largest B{C{distance#}}.
789 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
790 insufficient overlap for C{B{area}=True} or
791 no intersection or all (near-)concentric for
792 C{B{area}=False}.
794 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
796 @raise ValueError: Coincident B{C{point2}} or B{C{point3}} or invalid
797 B{C{distance1}}, B{C{distance2}}, B{C{distance3}}
798 or B{C{radius}}.
799 '''
800 return _trilaterate5(self, distance1,
801 self.others(point2=point2), distance2,
802 self.others(point3=point3), distance3,
803 area=area, radius=radius, eps=eps, wrap=wrap)
806_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
809def areaOf(points, radius=R_M, wrap=False): # was=True
810 '''Calculate the area of a (spherical) polygon or composite
811 (with the pointsjoined by great circle arcs).
813 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
814 or L{BooleanGH}).
815 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
816 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
817 or C{None}.
818 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{points}}
819 (C{bool}).
821 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
822 or C{radians} if B{C{radius}} is C{None}).
824 @raise PointsError: Insufficient number of B{C{points}}.
826 @raise TypeError: Some B{C{points}} are not L{LatLon}.
828 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
830 @note: The area is based on I{Karney}'s U{'Area of a spherical
831 polygon'<https://MathOverflow.net/questions/97711/
832 the-area-of-spherical-polygons>}, 3rd Answer.
834 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
835 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
836 L{ellipsoidalKarney.areaOf}.
838 @example:
840 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
841 >>> areaOf(b) # 8666058750.718977
843 >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1)
844 >>> areaOf(c) # 6.18e9
845 '''
846 if _MODS.booleans.isBoolean(points):
847 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
849 _at2, _t_2 = atan2, tan_2
850 _un, _w180 = unrollPI, wrap180
852 Ps = _T00.PointsIter(points, loop=1, wrap=wrap)
853 p1 = p2 = Ps[0]
854 a1, b1 = p1.philam
855 ta1, z1 = _t_2(a1), None
857 A = Fsum() # mean phi
858 R = Fsum() # see L{pygeodesy.excessKarney_}
859 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
860 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
861 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
862 D = Fsum()
863 for i, p2 in Ps.enumerate(closed=True):
864 a2, b2 = p2.philam
865 db, b2 = _un(b1, b2, wrap=wrap and not Ps.looped)
866 A += a2
867 ta2 = _t_2(a2)
868 tdb = _t_2(db, points=i)
869 R += _at2(tdb * (ta1 + ta2),
870 _1_0 + ta1 * ta2)
871 ta1, b1 = ta2, b2
873 if not p2.isequalTo(p1, eps=EPS):
874 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
875 if z1 is not None:
876 D += _w180(z - z1) # (z - z1 + 540) ...
877 D += _w180(z2 - z) # (z2 - z + 540) % 360 - 180
878 p1, z1 = p2, z2
880 R = abs(R * _2_0)
881 if abs(D) < _90_0: # ispolar(points)
882 R = abs(R - PI2)
883 if radius:
884 a = degrees(A.fover(len(A))) # mean lat
885 R *= _mean_radius(radius, a)**2
886 return float(R)
889def _destination2(a, b, r, t):
890 '''(INTERNAL) Destination lat- and longitude in C{radians}.
892 @arg a: Latitude (C{radians}).
893 @arg b: Longitude (C{radians}).
894 @arg r: Angular distance (C{radians}).
895 @arg t: Bearing (compass C{radians}).
897 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
898 '''
899 # see <https://www.EdWilliams.org/avform.htm#LL>
900 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
901 ca *= sr
903 a = asin1(ct * ca + cr * sa)
904 d = atan2(st * ca, cr - sa * sin(a))
905 # note, in EdWilliams.org/avform.htm W is + and E is -
906 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
909def _int3d2(s, end, wrap, _i_, Vector, hs):
910 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
911 # and similar logic in .ellipsoidalBaseDI._intersect3
912 a1, b1 = s.philam
914 if isscalar(end): # bearing, get pseudo-end point
915 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
916 else: # must be a point
917 s.others(end, name=_end_ + _i_)
918 hs.append(end.height)
919 a2, b2 = end.philam
920 if wrap:
921 a2, b2 = _Wrap.philam(a2, b2)
923 db, b2 = unrollPI(b1, b2, wrap=wrap)
924 if max(fabs(db), fabs(a2 - a1)) < EPS:
925 raise _ValueError(_SPACE_(_line_ + _i_, _null_))
926 # note, in EdWilliams.org/avform.htm W is + and E is -
927 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
928 -(b1 + b2) * _0_5)
929 cb21 *= sin(a1 - a2) # sa21
930 sb21 *= sin(a1 + a2) # sa12
931 x = Vector(sb12 * cb21 - cb12 * sb21,
932 cb12 * cb21 + sb12 * sb21,
933 cos(a1) * cos(a2) * sin(db)) # ll=start
934 return x.unit(), (db, (a2 - a1)) # negated d
937def _intdot(ds, a1, b1, a, b, wrap):
938 # compute dot product ds . (-b + b1, a - a1)
939 db, _ = unrollPI(b1, b, wrap=wrap)
940 return fdot(ds, db, a - a1)
943def intersecant2(center, circle, point, bearing, radius=R_M, exact=False,
944 height=None, wrap=False): # was=True
945 '''Compute the intersections of a circle and a line.
947 @arg center: Center of the circle (L{LatLon}).
948 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
949 or a point on the circle (L{LatLon}).
950 @arg point: A point in- or outside the circle (L{LatLon}).
951 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or
952 a second point on the line (L{LatLon}).
953 @kwarg radius: Mean earth radius (C{meter}, conventionally).
954 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
955 destination and distance, if C{False} use the basic
956 rhumb methods (C{bool}) or if C{None} use the I{great
957 circle} methods.
958 @kwarg height: Optional height for the intersection points (C{meter},
959 conventionally) or C{None}.
960 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}}
961 and the B{C{circle}} and B{C{bearing}} points (C{bool}).
963 @return: 2-Tuple of the intersection points (representing a chord),
964 each an instance of this class. For a tangent line, both
965 points are the same instance, the B{C{point}} or wrapped
966 or I{normalized}.
968 @raise IntersectionError: The circle and line do not intersect.
970 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
971 B{C{circle}} or B{C{bearing}} invalid.
973 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
974 B{C{exact}} or B{C{height}}.
975 '''
976 c = _T00.others(center=center)
977 p = _T00.others(point=point)
978 try:
979 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact,
980 height=height, wrap=wrap)
981 except (TypeError, ValueError) as x:
982 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing, exact=exact)
985def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
986 LatLon=None, **LatLon_kwds):
987 # (INTERNAL) Intersect two (spherical) lines, see L{intersection}
988 # above, separated to allow callers to embellish any exceptions
990 s1, s2 = start1, start2
991 if wrap:
992 s2 = _Wrap.point(s2)
993 hs = [s1.height, s2.height]
995 a1, b1 = s1.philam
996 a2, b2 = s2.philam
997 db, b2 = unrollPI(b1, b2, wrap=wrap)
998 r12 = vincentys_(a2, a1, db)
999 if fabs(r12) < EPS: # [nearly] coincident points
1000 a, b = favg(a1, a2), favg(b1, b2)
1002 # see <https://www.EdWilliams.org/avform.htm#Intersection>
1003 elif isscalar(end1) and isscalar(end2): # both bearings
1004 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
1006 x1, x2 = (sr12 * ca1), (sr12 * ca2)
1007 if isnear0(x1) or isnear0(x2):
1008 raise IntersectionError(_parallel_)
1009 # handle domain error for equivalent longitudes,
1010 # see also functions asin_safe and acos_safe at
1011 # <https://www.EdWilliams.org/avform.htm#Math>
1012 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1)
1013 t21, t23 = acos1((sa1 - sa2 * cr12) / x2), radiansPI2(end2)
1014 if sin(db) > 0:
1015 t21 = PI2 - t21
1016 else:
1017 t12 = PI2 - t12
1018 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
1019 wrapPI(t21 - t23)) # angle 1-2-3)
1020 if isnear0(sx1) and isnear0(sx2):
1021 raise IntersectionError(_infinite_)
1022 sx3 = sx1 * sx2
1023# XXX if sx3 < 0:
1024# XXX raise ValueError(_ambiguous_)
1025 x3 = acos1(cr12 * sx3 - cx2 * cx1)
1026 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
1028 a, b = _destination2(a1, b1, r13, t13)
1029 # like .ellipsoidalBaseDI,_intersect3, if this intersection
1030 # is "before" the first point, use the antipodal intersection
1031 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
1032 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1034 else: # end point(s) or bearing(s)
1035 _N_vector_ = _MODS.nvectorBase._N_vector_
1037 x1, d1 = _int3d2(s1, end1, wrap, _1_, _N_vector_, hs)
1038 x2, d2 = _int3d2(s2, end2, wrap, _2_, _N_vector_, hs)
1039 x = x1.cross(x2)
1040 if x.length < EPS: # [nearly] colinear or parallel lines
1041 raise IntersectionError(_colinear_)
1042 a, b = x.philam
1043 # choose intersection similar to sphericalNvector
1044 if not (_intdot(d1, a1, b1, a, b, wrap) *
1045 _intdot(d2, a2, b2, a, b, wrap)) > 0:
1046 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1048 h = fmean(hs) if height is None else Height(height)
1049 return _LL3Tuple(degrees90(a), degrees180(b), h,
1050 intersection, LatLon, LatLon_kwds)
1053def intersection(start1, end1, start2, end2, height=None, wrap=False,
1054 LatLon=LatLon, **LatLon_kwds):
1055 '''Compute the intersection point of two lines, each defined
1056 by two points or a start point and bearing from North.
1058 @arg start1: Start point of the first line (L{LatLon}).
1059 @arg end1: End point of the first line (L{LatLon}) or
1060 the initial bearing at the first start point
1061 (compass C{degrees360}).
1062 @arg start2: Start point of the second line (L{LatLon}).
1063 @arg end2: End point of the second line (L{LatLon}) or
1064 the initial bearing at the second start point
1065 (compass C{degrees360}).
1066 @kwarg height: Optional height for the intersection point,
1067 overriding the mean height (C{meter}).
1068 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1069 B{C{start2}} and both B{C{end*}} points (C{bool}).
1070 @kwarg LatLon: Optional class to return the intersection
1071 point (L{LatLon}) or C{None}.
1072 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1073 arguments, ignored if C{B{LatLon} is None}.
1075 @return: The intersection point as a (B{C{LatLon}}) or if
1076 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon,
1077 height)}. An alternate intersection point might
1078 be the L{antipode} to the returned result.
1080 @raise IntersectionError: Ambiguous or infinite intersection
1081 or colinear, parallel or otherwise
1082 non-intersecting lines.
1084 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}}
1085 or B{C{end2}} point not L{LatLon}.
1087 @raise ValueError: Invalid B{C{height}} or C{null} line.
1089 @example:
1091 >>> p = LatLon(51.8853, 0.2545)
1092 >>> s = LatLon(49.0034, 2.5735)
1093 >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E'
1094 '''
1095 s1 = _T00.others(start1=start1)
1096 s2 = _T00.others(start2=start2)
1097 try:
1098 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap,
1099 LatLon=LatLon, **LatLon_kwds)
1100 except (TypeError, ValueError) as x:
1101 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1104def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1105 height=None, wrap=False, # was=True
1106 LatLon=LatLon, **LatLon_kwds):
1107 '''Compute the intersection points of two circles each defined
1108 by a center point and a radius.
1110 @arg center1: Center of the first circle (L{LatLon}).
1111 @arg rad1: Radius of the first circle (C{meter} or C{radians},
1112 see B{C{radius}}).
1113 @arg center2: Center of the second circle (L{LatLon}).
1114 @arg rad2: Radius of the second circle (C{meter} or C{radians},
1115 see B{C{radius}}).
1116 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1117 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1118 @kwarg eps: Required overlap (C{meter} or C{radians}, see
1119 B{C{radius}}).
1120 @kwarg height: Optional height for the intersection points (C{meter},
1121 conventionally) or C{None} for the I{"radical height"}
1122 at the I{radical line} between both centers.
1123 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{center2}}
1124 (C{bool}).
1125 @kwarg LatLon: Optional class to return the intersection
1126 points (L{LatLon}) or C{None}.
1127 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1128 arguments, ignored if C{B{LatLon} is None}.
1130 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1131 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1132 lon, height)}. For abutting circles, both intersection
1133 points are the same instance, aka the I{radical center}.
1135 @raise IntersectionError: Concentric, antipodal, invalid or
1136 non-intersecting circles.
1138 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1140 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1141 B{C{eps}} or B{C{height}}.
1143 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1145 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1146 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1147 '''
1148 c1 = _T00.others(center1=center1)
1149 c2 = _T00.others(center2=center2)
1150 try:
1151 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1152 height=height, wrap=wrap,
1153 LatLon=LatLon, **LatLon_kwds)
1154 except (TypeError, ValueError) as x:
1155 raise _xError(x, center1=center1, rad1=rad1,
1156 center2=center2, rad2=rad2, wrap=wrap)
1159def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1160 height=None, too_d=None, wrap=False, # was=True
1161 LatLon=LatLon, **LatLon_kwds):
1162 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1163 # above, separated to allow callers to embellish any exceptions
1165 def _dest3(bearing, h):
1166 a, b = _destination2(a1, b1, r1, bearing)
1167 return _LL3Tuple(degrees90(a), degrees180(b), h,
1168 intersections2, LatLon, LatLon_kwds)
1170 a1, b1 = c1.philam
1171 a2, b2 = c2.philam
1172 if wrap:
1173 a2, b2 = _Wrap.philam(a2, b2)
1175 r1, r2, f = _rads3(rad1, rad2, radius)
1176 if f: # swapped radii, swap centers
1177 a1, a2 = a2, a1 # PYCHOK swap!
1178 b1, b2 = b2, b1 # PYCHOK swap!
1180 db, b2 = unrollPI(b1, b2, wrap=wrap)
1181 d = vincentys_(a2, a1, db) # radians
1182 if d < max(r1 - r2, EPS):
1183 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1185 r = eps if radius is None else (m2radians(
1186 eps, radius=radius) if eps else _0_0)
1187 if r < _0_0:
1188 raise _ValueError(eps=r)
1190 x = fsumf_(r1, r2, -d) # overlap
1191 if x > max(r, EPS):
1192 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1193 x = sd * sr1
1194 if isnear0(x):
1195 raise _ValueError(_invalid_)
1196 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1198 elif x < r: # PYCHOK no cover
1199 t = (d * radius) if too_d is None else too_d
1200 raise IntersectionError(_too_(_Fmt.distant(t)))
1202 if height is None: # "radical height"
1203 f = _radical2(d, r1, r2).ratio
1204 h = Height(favg(c1.height, c2.height, f=f))
1205 else:
1206 h = Height(height)
1208 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1209 if x < EPS4: # externally ...
1210 r = _dest3(b, h)
1211 elif x > _PI_EPS4: # internally ...
1212 r = _dest3(b + PI, h)
1213 else:
1214 return _dest3(b + x, h), _dest3(b - x, h)
1215 return r, r # ... abutting circles
1218@deprecated_function
1219def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1220 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1221 '''
1222 return ispolar(points, wrap=wrap)
1225def _LL3Tuple(lat, lon, height, where, LatLon, LatLon_kwds):
1226 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1227 '''
1228 n = where.__name__
1229 if LatLon is None:
1230 r = LatLon3Tuple(lat, lon, height, name=n)
1231 else:
1232 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1233 r = LatLon(lat, lon, **kwds)
1234 return r
1237def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):
1238 '''Compute the I{geographic} mean of several points.
1240 @arg points: Points to be averaged (L{LatLon}[]).
1241 @kwarg height: Optional height at mean point, overriding the mean
1242 height (C{meter}).
1243 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}}
1244 (C{bool}).
1245 @kwarg LatLon: Optional class to return the mean point (L{LatLon})
1246 or C{None}.
1247 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1248 arguments, ignored if C{B{LatLon} is None}.
1250 @return: The geographic mean and height (B{C{LatLon}}) or a
1251 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}}
1252 is C{None}.
1254 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1256 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1257 '''
1258 def _N_vs(ps, w):
1259 Ps = _T00.PointsIter(ps, wrap=w)
1260 for p in Ps.iterate(closed=False):
1261 yield p._N_vector
1263 m = _MODS.nvectorBase
1264 # geographic, vectorial mean
1265 n = m.sumOf(_N_vs(points, wrap), h=height, Vector=m.NvectorBase)
1266 lat, lon, h = n.latlonheight
1267 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1270@deprecated_function
1271def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1272 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1274 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1275 point (L{LatLon}) on the polygon and the C{distance}
1276 between the C{closest} and the given B{C{point}}. The
1277 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1278 lon)} if B{C{LatLon}} is C{None} ...
1279 '''
1280 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1281 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1282 ll = LatLon2Tuple(ll.lat, ll.lon)
1283 return ll, d
1286def nearestOn3(point, points, closed=False, radius=R_M, wrap=False, adjust=True,
1287 limit=9, **LatLon_and_kwds):
1288 '''Locate the point on a path or polygon closest to a reference point.
1290 Distances are I{approximated} using function L{pygeodesy.equirectangular_},
1291 subject to the supplied B{C{options}}.
1293 @arg point: The reference point (L{LatLon}).
1294 @arg points: The path or polygon points (L{LatLon}[]).
1295 @kwarg closed: Optionally, close the polygon (C{bool}).
1296 @kwarg radius: Mean earth radius (C{meter}).
1297 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1298 B{C{points}} (C{bool}).
1299 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}).
1300 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}),
1301 default C{9 degrees} is about C{1,000 Kmeter} (for mean
1302 spherical earth radius L{R_KM}).
1303 @kwarg LatLon: Optional class to return the closest point (L{LatLon})
1304 or C{None}.
1305 @kwarg options: Optional keyword arguments for function
1306 L{pygeodesy.equirectangular_}.
1308 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1309 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1310 lon, height)} if B{C{LatLon}} is C{None}. The C{distance}
1311 is the L{pygeodesy.equirectangular_} distance between the
1312 C{closest} and the given B{C{point}} converted to C{meter},
1313 same units as B{C{radius}}. The C{angle} from the given
1314 B{C{point}} to the C{closest} is in compass C{degrees360},
1315 like function L{pygeodesy.compassAngle}. The C{height} is
1316 the (interpolated) height at the C{closest} point.
1318 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1319 see function L{pygeodesy.equirectangular_}.
1321 @raise PointsError: Insufficient number of B{C{points}}.
1323 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1325 @raise ValueError: Invalid B{C{radius}}.
1327 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}.
1328 '''
1329 t = _nearestOn5(point, points, closed=closed, wrap=wrap,
1330 adjust=adjust, limit=limit)
1331 d = degrees2m(t.distance, radius=radius)
1332 h = t.height
1333 n = nearestOn3.__name__
1335 kwds = _xkwds(LatLon_and_kwds, height=h, name=n)
1336 LL = _xkwds_pop(kwds, LatLon=LatLon)
1337 r = LatLon3Tuple(t.lat, t.lon, h, name=n) if LL is None else \
1338 LL(t.lat, t.lon, **kwds)
1339 return NearestOn3Tuple(r, d, t.angle, name=n)
1342def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1343 '''Compute the perimeter of a (spherical) polygon or composite
1344 (with great circle arcs joining the points).
1346 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1347 or L{BooleanGH}).
1348 @kwarg closed: Optionally, close the polygon (C{bool}).
1349 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1350 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1351 B{C{points}} (C{bool}).
1353 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1354 or C{radians} if B{C{radius}} is C{None}).
1356 @raise PointsError: Insufficient number of B{C{points}}.
1358 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1360 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1361 C{B{points}} a composite.
1363 @note: Distances are based on function L{pygeodesy.vincentys_}.
1365 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf}
1366 and L{ellipsoidalKarney.perimeterOf}.
1367 '''
1368 def _rads(ps, c, w): # angular edge lengths in radians
1369 Ps = _T00.PointsIter(ps, loop=1, wrap=w)
1370 a1, b1 = Ps[0].philam
1371 for p in Ps.iterate(closed=c):
1372 a2, b2 = p.philam
1373 db, b2 = unrollPI(b1, b2, wrap=w and not (c and Ps.looped))
1374 yield vincentys_(a2, a1, db)
1375 a1, b1 = a2, b2
1377 if _MODS.booleans.isBoolean(points):
1378 if not closed:
1379 raise _ValueError(closed=closed, points=_composite_)
1380 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1381 else:
1382 r = fsum(_rads(points, closed, wrap), floats=True)
1383 return _r2m(r, radius)
1386def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1387 excess=excessAbc_,
1388 wrap=False):
1389 '''Compute the angles, sides, and area of a (spherical) triangle.
1391 @arg latA: First corner latitude (C{degrees}).
1392 @arg lonA: First corner longitude (C{degrees}).
1393 @arg latB: Second corner latitude (C{degrees}).
1394 @arg lonB: Second corner longitude (C{degrees}).
1395 @arg latC: Third corner latitude (C{degrees}).
1396 @arg lonC: Third corner longitude (C{degrees}).
1397 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1398 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1399 or C{None}.
1400 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1401 L{excessGirard_} or L{excessLHuilier_}).
1402 @kwarg wrap: If C{True}, wrap and L{pygeodesy.unroll180}
1403 longitudes (C{bool}).
1405 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1406 spherical angles C{A}, C{B} and C{C}, angular sides
1407 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1408 in I{square} C{meter} or same units as B{C{radius}}
1409 I{squared} or if C{B{radius}=0} or C{None}, a
1410 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in
1411 C{radians} with the I{spherical excess} C{E} as the
1412 C{unit area} in C{radians}.
1413 '''
1414 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA),
1415 Phi_(latB=latB), Lam_(lonB=lonB),
1416 Phi_(latC=latC), Lam_(lonC=lonC),
1417 excess=excess, wrap=wrap)
1418 return _t7Tuple(t, radius)
1421def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_,
1422 wrap=False):
1423 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1424 excess} of a (spherical) triangle.
1426 @arg phiA: First corner latitude (C{radians}).
1427 @arg lamA: First corner longitude (C{radians}).
1428 @arg phiB: Second corner latitude (C{radians}).
1429 @arg lamB: Second corner longitude (C{radians}).
1430 @arg phiC: Third corner latitude (C{radians}).
1431 @arg lamC: Third corner longitude (C{radians}).
1432 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1433 L{excessGirard_} or L{excessLHuilier_}).
1434 @kwarg wrap: If C{True}, L{pygeodesy.unrollPI} the
1435 longitudinal deltas (C{bool}).
1437 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1438 spherical angles C{A}, C{B} and C{C}, angular sides
1439 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1440 I{spherical excess} C{E}, all in C{radians}.
1441 '''
1442 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1443 d, _ = unrollPI(lamB, lamC, wrap=w)
1444 a = vincentys_(phiC, phiB, d)
1445 return a, (phiB, lamB, phiC, lamC, phiA, lamA) # rotate A, B, C
1447 def _A_r(a, sa, ca, sb, cb, sc, cc):
1448 s = sb * sc
1449 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1450 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2_'s
1452 # notation: side C{a} is oposite to corner C{A}, etc.
1453 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1454 b, r = _a_r(wrap, *r)
1455 c, _ = _a_r(wrap, *r)
1457 A, r = _A_r(a, *sincos2_(a, b, c))
1458 B, r = _A_r(b, *r)
1459 C, _ = _A_r(c, *r)
1461 D = fsumf_(PI2, -a, -b, -c) # deficit aka defect
1462 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else (
1463 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else
1464 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b))))
1466 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1469def _t7Tuple(t, radius):
1470 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1471 '''
1472 if radius: # not in (None, _0_0)
1473 r = radius if isscalar(radius) else \
1474 _ellipsoidal_datum(radius).ellipsoid.Rmean
1475 A, B, C = map1(degrees, t.A, t.B, t.C)
1476 t = Triangle7Tuple(A, (r * t.a),
1477 B, (r * t.b),
1478 C, (r * t.c), t.E * r**2)
1479 return t
1482__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1483 areaOf, # functions
1484 intersecant2, intersection, intersections2, ispolar,
1485 isPoleEnclosedBy, # DEPRECATED, use ispolar
1486 meanOf,
1487 nearestOn2, nearestOn3,
1488 perimeterOf,
1489 sumOf, # XXX == vector3d.sumOf
1490 triangle7, triangle8_)
1492# **) MIT License
1493#
1494# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1495#
1496# Permission is hereby granted, free of charge, to any person obtaining a
1497# copy of this software and associated documentation files (the "Software"),
1498# to deal in the Software without restriction, including without limitation
1499# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1500# and/or sell copies of the Software, and to permit persons to whom the
1501# Software is furnished to do so, subject to the following conditions:
1502#
1503# The above copyright notice and this permission notice shall be included
1504# in all copies or substantial portions of the Software.
1505#
1506# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1507# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1508# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1509# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1510# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1511# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1512# OTHER DEALINGS IN THE SOFTWARE.