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