Coverage for pygeodesy/sphericalTrigonometry.py: 94%
390 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-12-02 13:46 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2023-12-02 13:46 -0500
2# -*- coding: utf-8 -*-
4u'''Spherical, C{trigonometry}-based geodesy.
6Trigonometric classes geodetic (lat-/longitude) L{LatLon} and
7geocentric (ECEF) L{Cartesian} and functions L{areaOf}, L{intersection},
8L{intersections2}, L{isPoleEnclosedBy}, L{meanOf}, L{nearestOn3} and
9L{perimeterOf}, I{all spherical}.
11Pure Python implementation of geodetic (lat-/longitude) methods using
12spherical trigonometry, transcoded from JavaScript originals by
13I{(C) Chris Veness 2011-2016} published under the same MIT Licence**, see
14U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}.
15'''
16# make sure int/int division yields float quotient, see .basics
17from __future__ import division as _; del _ # PYCHOK semicolon
19from pygeodesy.basics import copysign0, 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 _m2radians, CartesianSphericalBase, \
45 _intersecant2, LatLonSphericalBase, \
46 _rads3, _radians2m, _trilaterate5
47# from pygeodesy.streprs import Fmt as _Fmt # from .points XXX shadowed
48from pygeodesy.units import Bearing_, Height, Lam_, Phi_, Radius_, Scalar
49from pygeodesy.utily import acos1, asin1, atan1d, atan2d, degrees90, degrees180, \
50 degrees2m, m2radians, radiansPI2, sincos2_, tan_2, \
51 unrollPI, _unrollon, _unrollon3, _Wrap, wrap180, wrapPI
52from pygeodesy.vector3d import sumOf, Vector3d
54from math import asin, atan2, cos, degrees, fabs, radians, sin
56__all__ = _ALL_LAZY.sphericalTrigonometry
57__version__ = '23.10.24'
59_PI_EPS4 = PI - EPS4
60if _PI_EPS4 >= PI:
61 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4)
64class Cartesian(CartesianSphericalBase):
65 '''Extended to convert geocentric, L{Cartesian} points to
66 spherical, geodetic L{LatLon}.
67 '''
69 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon
70 '''Convert this cartesian point to a C{spherical} geodetic point.
72 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
73 arguments. Use C{B{LatLon}=...} to override
74 this L{LatLon} class or specify C{B{LatLon}=None}.
76 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is C{None},
77 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
78 with C{C} and C{M} if available.
80 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
81 '''
82 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
83 return CartesianSphericalBase.toLatLon(self, **kwds)
86class LatLon(LatLonSphericalBase):
87 '''New point on spherical model earth model.
89 @example:
91 >>> p = LatLon(52.205, 0.119) # height=0
92 '''
94 def _ab1_ab2_db5(self, other, wrap):
95 '''(INTERNAL) Helper for several methods.
96 '''
97 a1, b1 = self.philam
98 a2, b2 = self.others(other, up=2).philam
99 if wrap:
100 a2, b2 = _Wrap.philam(a2, b2)
101 db, b2 = unrollPI(b1, b2, wrap=wrap)
102 else: # unrollPI shortcut
103 db = b2 - b1
104 return a1, b1, a2, b2, db
106 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
107 '''Compute the (signed) distance from the start to the closest
108 point on the great circle line defined by a start and an
109 end point.
111 That is, if a perpendicular is drawn from this point to the
112 great circle line, the along-track distance is the distance
113 from the start point to the point where the perpendicular
114 crosses the line.
116 @arg start: Start point of the great circle line (L{LatLon}).
117 @arg end: End point of the great circle line (L{LatLon}).
118 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
119 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
120 the B{C{start}} and B{C{end}} point (C{bool}).
122 @return: Distance along the great circle line (C{radians}
123 if C{B{radius} is None} or C{meter}, same units
124 as B{C{radius}}), positive if I{after} the
125 B{C{start}} toward the B{C{end}} point of the
126 line, I{negative} if before or C{0} if at the
127 B{C{start}} point.
129 @raise TypeError: Invalid B{C{start}} or B{C{end}} point.
131 @raise ValueError: Invalid B{C{radius}}.
133 @example:
135 >>> p = LatLon(53.2611, -0.7972)
137 >>> s = LatLon(53.3206, -1.7297)
138 >>> e = LatLon(53.1887, 0.1334)
139 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58
140 '''
141 r, x, b = self._a_x_b3(start, end, radius, wrap)
142 cx = cos(x)
143 return _0_0 if isnear0(cx) else \
144 _radians2m(copysign0(acos1(cos(r) / cx), cos(b)), radius)
146 def _a_x_b3(self, start, end, radius, wrap):
147 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo.
148 '''
149 s = self.others(start=start)
150 e = self.others(end=end)
151 s, e, w = _unrollon3(self, s, e, wrap)
153 r = Radius_(radius)
154 r = s.distanceTo(self, r, wrap=w) / r
156 b = radians(s.initialBearingTo(self, wrap=w)
157 - s.initialBearingTo(e, wrap=w))
158 x = asin(sin(r) * sin(b))
159 return r, x, -b
161 @deprecated_method
162 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover
163 '''DEPRECATED, use method L{initialBearingTo}.
164 '''
165 return self.initialBearingTo(other, wrap=wrap, raiser=raiser)
167 def crossingParallels(self, other, lat, wrap=False):
168 '''Return the pair of meridians at which a great circle defined
169 by this and an other point crosses the given latitude.
171 @arg other: The other point defining great circle (L{LatLon}).
172 @arg lat: Latitude at the crossing (C{degrees}).
173 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
174 B{C{other}} point (C{bool}).
176 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
177 C{None} if the great circle doesn't reach B{C{lat}}.
178 '''
179 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
180 sa, ca, sa1, ca1, \
181 sa2, ca2, sdb, cdb = sincos2_(radians(lat), a1, a2, db)
182 sa1 *= ca2 * ca
184 x = sa1 * sdb
185 y = sa1 * cdb - ca1 * sa2 * ca
186 z = ca1 * sdb * ca2 * sa
188 h = hypot(x, y)
189 if h < EPS or fabs(z) > h: # PYCHOK no cover
190 return None # great circle doesn't reach latitude
192 m = atan2(-y, x) + b1 # longitude at max latitude
193 d = acos1(z / h) # delta longitude to intersections
194 return degrees180(m - d), degrees180(m + d)
196 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
197 '''Compute the (signed) distance from this point to
198 the great circle defined by a start and an end point.
200 @arg start: Start point of the great circle line (L{LatLon}).
201 @arg end: End point of the great circle line (L{LatLon}).
202 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
203 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
204 the B{C{start}} and B{C{end}} point (C{bool}).
206 @return: Distance to the great circle (C{radians} if
207 B{C{radius}} or C{meter}, same units as
208 B{C{radius}}), I{negative} if to the left or
209 I{positive} if to the right of the line.
211 @raise TypeError: If B{C{start}} or B{C{end}} is not L{LatLon}.
213 @raise ValueError: Invalid B{C{radius}}.
215 @example:
217 >>> p = LatLon(53.2611, -0.7972)
219 >>> s = LatLon(53.3206, -1.7297)
220 >>> e = LatLon(53.1887, 0.1334)
221 >>> d = p.crossTrackDistanceTo(s, e) # -307.5
222 '''
223 _, x, _ = self._a_x_b3(start, end, radius, wrap)
224 return _radians2m(x, radius)
226 def destination(self, distance, bearing, radius=R_M, height=None):
227 '''Locate the destination from this point after having
228 travelled the given distance on the given initial bearing.
230 @arg distance: Distance travelled (C{meter}, same units as
231 B{C{radius}}).
232 @arg bearing: Bearing from this point (compass C{degrees360}).
233 @kwarg radius: Mean earth radius (C{meter}).
234 @kwarg height: Optional height at destination (C{meter}, same
235 units a B{C{radius}}).
237 @return: Destination point (L{LatLon}).
239 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
240 B{C{radius}} or B{C{height}}.
242 @example:
244 >>> p1 = LatLon(51.4778, -0.0015)
245 >>> p2 = p1.destination(7794, 300.7)
246 >>> p2.toStr() # '51.5135°N, 000.0983°W'
247 '''
248 a, b = self.philam
249 r, t = _m2radians(distance, radius, low=None), Bearing_(bearing)
251 a, b = _destination2(a, b, r, t)
252 h = self._heigHt(height)
253 return self.classof(degrees90(a), degrees180(b), height=h)
255 def distanceTo(self, other, radius=R_M, wrap=False):
256 '''Compute the (angular) distance from this to an other point.
258 @arg other: The other point (L{LatLon}).
259 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
260 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
261 the B{C{other}} point (C{bool}).
263 @return: Distance between this and the B{C{other}} point
264 (C{meter}, same units as B{C{radius}} or
265 C{radians} if B{C{radius}} is C{None}).
267 @raise TypeError: The B{C{other}} point is not L{LatLon}.
269 @raise ValueError: Invalid B{C{radius}}.
271 @example:
273 >>> p1 = LatLon(52.205, 0.119)
274 >>> p2 = LatLon(48.857, 2.351);
275 >>> d = p1.distanceTo(p2) # 404300
276 '''
277 a1, _, a2, _, db = self._ab1_ab2_db5(other, wrap)
278 return _radians2m(vincentys_(a2, a1, db), radius)
280# @Property_RO
281# def Ecef(self):
282# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
283# '''
284# return _MODS.ecef.EcefKarney
286 def greatCircle(self, bearing, Vector=Vector3d, **Vector_kwds):
287 '''Compute the vector normal to great circle obtained by heading
288 on the given initial bearing from this point.
290 Direction of vector is such that initial bearing vector
291 b = c × n, where n is an n-vector representing this point.
293 @arg bearing: Bearing from this point (compass C{degrees360}).
294 @kwarg Vector: Vector class to return the great circle,
295 overriding the default L{Vector3d}.
296 @kwarg Vector_kwds: Optional, additional keyword argunents
297 for B{C{Vector}}.
299 @return: Vector representing great circle (C{Vector}).
301 @raise ValueError: Invalid B{C{bearing}}.
303 @example:
305 >>> p = LatLon(53.3206, -1.7297)
306 >>> g = p.greatCircle(96.0)
307 >>> g.toStr() # (-0.794, 0.129, 0.594)
308 '''
309 a, b = self.philam
310 sa, ca, sb, cb, st, ct = sincos2_(a, b, Bearing_(bearing))
312 return Vector(sb * ct - cb * sa * st,
313 -cb * ct - sb * sa * st,
314 ca * st, **Vector_kwds) # XXX .unit()?
316 def initialBearingTo(self, other, wrap=False, raiser=False):
317 '''Compute the initial bearing (forward azimuth) from this
318 to an other point.
320 @arg other: The other point (spherical L{LatLon}).
321 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
322 the B{C{other}} point (C{bool}).
323 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}),
324 use C{B{raiser}=True} for behavior like
325 C{sphericalNvector.LatLon.initialBearingTo}.
327 @return: Initial bearing (compass C{degrees360}).
329 @raise CrossError: If this and the B{C{other}} point coincide,
330 provided both B{C{raiser}} is C{True} and
331 L{pygeodesy.crosserrors} is C{True}.
333 @raise TypeError: The B{C{other}} point is not L{LatLon}.
335 @example:
337 >>> p1 = LatLon(52.205, 0.119)
338 >>> p2 = LatLon(48.857, 2.351)
339 >>> b = p1.initialBearingTo(p2) # 156.2
340 '''
341 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
342 # XXX behavior like sphericalNvector.LatLon.initialBearingTo
343 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(db)) < EPS:
344 raise CrossError(_point_, self, other=other, wrap=wrap, txt=_coincident_)
346 return degrees(bearing_(a1, b1, a2, b2, final=False))
348 def intermediateTo(self, other, fraction, height=None, wrap=False):
349 '''Locate the point at given fraction between (or along) this
350 and an other point.
352 @arg other: The other point (L{LatLon}).
353 @arg fraction: Fraction between both points (C{scalar},
354 0.0 at this and 1.0 at the other point).
355 @kwarg height: Optional height, overriding the intermediate
356 height (C{meter}).
357 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
358 B{C{other}} point (C{bool}).
360 @return: Intermediate point (L{LatLon}).
362 @raise TypeError: The B{C{other}} point is not L{LatLon}.
364 @raise ValueError: Invalid B{C{fraction}} or B{C{height}}.
366 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
368 @example:
370 >>> p1 = LatLon(52.205, 0.119)
371 >>> p2 = LatLon(48.857, 2.351)
372 >>> p = p1.intermediateTo(p2, 0.25) # 51.3721°N, 000.7073°E
373 '''
374 p = self
375 f = Scalar(fraction=fraction)
376 if not isnear0(f):
377 p = p.others(other)
378 if wrap:
379 p = _Wrap.point(p)
380 if not isnear1(f): # and not near0
381 a1, b1 = self.philam
382 a2, b2 = p.philam
383 db, b2 = unrollPI(b1, b2, wrap=wrap)
384 r = vincentys_(a2, a1, db)
385 sr = sin(r)
386 if isnon0(sr):
387 sa1, ca1, sa2, ca2, \
388 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2)
390 t = f * r
391 a = sin(r - t) # / sr superflous
392 b = sin( t) # / sr superflous
394 x = a * ca1 * cb1 + b * ca2 * cb2
395 y = a * ca1 * sb1 + b * ca2 * sb2
396 z = a * sa1 + b * sa2
398 a = atan1d(z, hypot(x, y))
399 b = atan2d(y, x)
401 else: # PYCHOK no cover
402 a = degrees90( favg(a1, a2, f=f)) # coincident
403 b = degrees180(favg(b1, b2, f=f))
405 h = self._havg(other, f=f, h=height)
406 p = self.classof(a, b, height=h)
407 return p
409 def intersection(self, end1, other, end2, height=None, wrap=False):
410 '''Compute the intersection point of two lines, each defined by
411 two points or a start point and bearing from North.
413 @arg end1: End point of this line (L{LatLon}) or the initial
414 bearing at this point (compass C{degrees360}).
415 @arg other: Start point of the other line (L{LatLon}).
416 @arg end2: End point of the other line (L{LatLon}) or the
417 initial bearing at the B{C{other}} point (compass
418 C{degrees360}).
419 @kwarg height: Optional height for intersection point,
420 overriding the mean height (C{meter}).
421 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
422 B{C{start2}} and both B{C{end*}} points (C{bool}).
424 @return: The intersection point (L{LatLon}). An alternate
425 intersection point might be the L{antipode} to
426 the returned result.
428 @raise IntersectionError: Ambiguous or infinite intersection
429 or colinear, parallel or otherwise
430 non-intersecting lines.
432 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}}
433 or B{C{end2}} not C{scalar} nor L{LatLon}.
435 @raise ValueError: Invalid B{C{height}} or C{null} line.
437 @example:
439 >>> p = LatLon(51.8853, 0.2545)
440 >>> s = LatLon(49.0034, 2.5735)
441 >>> i = p.intersection(108.547, s, 32.435) # '50.9078°N, 004.5084°E'
442 '''
443 try:
444 s2 = self.others(other)
445 return _intersect(self, end1, s2, end2, height=height, wrap=wrap,
446 LatLon=self.classof)
447 except (TypeError, ValueError) as x:
448 raise _xError(x, start1=self, end1=end1,
449 other=other, end2=end2, wrap=wrap)
451 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0,
452 height=None, wrap=True):
453 '''Compute the intersection points of two circles, each defined
454 by a center point and radius.
456 @arg rad1: Radius of the this circle (C{meter} or C{radians},
457 see B{C{radius}}).
458 @arg other: Center point of the other circle (L{LatLon}).
459 @arg rad2: Radius of the other circle (C{meter} or C{radians},
460 see B{C{radius}}).
461 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
462 B{C{rad2}} and B{C{eps}} are given in C{radians}).
463 @kwarg eps: Required overlap (C{meter} or C{radians}, see
464 B{C{radius}}).
465 @kwarg height: Optional height for the intersection points (C{meter},
466 conventionally) or C{None} for the I{"radical height"}
467 at the I{radical line} between both centers.
468 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
469 B{C{other}} point (C{bool}).
471 @return: 2-Tuple of the intersection points, each a L{LatLon}
472 instance. For abutting circles, both intersection
473 points are the same instance, aka the I{radical center}.
475 @raise IntersectionError: Concentric, antipodal, invalid or
476 non-intersecting circles.
478 @raise TypeError: If B{C{other}} is not L{LatLon}.
480 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
481 B{C{eps}} or B{C{height}}.
482 '''
483 try:
484 c2 = self.others(other)
485 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps,
486 height=height, wrap=wrap,
487 LatLon=self.classof)
488 except (TypeError, ValueError) as x:
489 raise _xError(x, center=self, rad1=rad1,
490 other=other, rad2=rad2, wrap=wrap)
492 @deprecated_method
493 def isEnclosedBy(self, points): # PYCHOK no cover
494 '''DEPRECATED, use method C{isenclosedBy}.'''
495 return self.isenclosedBy(points)
497 def isenclosedBy(self, points, wrap=False):
498 '''Check whether a (convex) polygon or composite encloses this point.
500 @arg points: The polygon points or composite (L{LatLon}[],
501 L{BooleanFHP} or L{BooleanGH}).
502 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
503 B{C{points}} (C{bool}).
505 @return: C{True} if this point is inside the polygon or
506 composite, C{False} otherwise.
508 @raise PointsError: Insufficient number of B{C{points}}.
510 @raise TypeError: Some B{C{points}} are not L{LatLon}.
512 @raise ValueError: Invalid B{C{points}}, non-convex polygon.
514 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
515 and L{pygeodesy.ispolar} especially if the B{C{points}} may
516 enclose a pole or wrap around the earth I{longitudinally}.
517 '''
518 if _MODS.booleans.isBoolean(points):
519 return points._encloses(self.lat, self.lon, wrap=wrap)
521 Ps = self.PointsIter(points, loop=2, dedup=True, wrap=wrap)
522 n0 = self._N_vector
524 v2 = Ps[0]._N_vector
525 p1 = Ps[1]
526 v1 = p1._N_vector
527 # check whether this point on same side of all
528 # polygon edges (to the left or right depending
529 # on the anti-/clockwise polygon direction)
530 gc1 = v2.cross(v1)
531 t0 = gc1.angleTo(n0) > PI_2
532 s0 = None
533 # get great-circle vector for each edge
534 for i, p2 in Ps.enumerate(closed=True):
535 if wrap and not Ps.looped:
536 p2 = _unrollon(p1, p2)
537 p1 = p2
538 v2 = p2._N_vector
539 gc = v1.cross(v2)
540 t = gc.angleTo(n0) > PI_2
541 if t != t0: # different sides of edge i
542 return False # outside
544 # check for convex polygon: angle between
545 # gc vectors, signed by direction of n0
546 # (otherwise the test above is not reliable)
547 s = signOf(gc1.angleTo(gc, vSign=n0))
548 if s != s0:
549 if s0 is None:
550 s0 = s
551 else:
552 t = _Fmt.SQUARE(points=i)
553 raise _ValueError(t, p2, wrap=wrap, txt=_not_(_convex_))
554 gc1, v1 = gc, v2
556 return True # inside
558 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
559 '''Find the midpoint between this and an other point.
561 @arg other: The other point (L{LatLon}).
562 @kwarg height: Optional height for midpoint, overriding
563 the mean height (C{meter}).
564 @kwarg fraction: Midpoint location from this point (C{scalar}),
565 may be negative or greater than 1.0.
566 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
567 B{C{other}} point (C{bool}).
569 @return: Midpoint (L{LatLon}).
571 @raise TypeError: The B{C{other}} point is not L{LatLon}.
573 @raise ValueError: Invalid B{C{height}}.
575 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
577 @example:
579 >>> p1 = LatLon(52.205, 0.119)
580 >>> p2 = LatLon(48.857, 2.351)
581 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
582 '''
583 if fraction is _0_5:
584 # see <https://MathForum.org/library/drmath/view/51822.html>
585 a1, b, a2, _, db = self._ab1_ab2_db5(other, wrap)
586 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db)
588 x = ca2 * cdb + ca1
589 y = ca2 * sdb
591 a = atan1d(sa1 + sa2, hypot(x, y))
592 b = degrees180(b + atan2(y, x))
594 h = self._havg(other, h=height)
595 r = self.classof(a, b, height=h)
596 else:
597 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
598 return r
600 def nearestOn(self, point1, point2, radius=R_M, **wrap_adjust_limit):
601 '''Locate the point between two points closest to this point.
603 Distances are approximated by function L{pygeodesy.equirectangular_},
604 subject to the supplied B{C{options}}.
606 @arg point1: Start point (L{LatLon}).
607 @arg point2: End point (L{LatLon}).
608 @kwarg radius: Mean earth radius (C{meter}).
609 @kwarg wrap_adjust_limit: Optional keyword arguments for functions
610 L{sphericalTrigonometry.nearestOn3} and
611 L{pygeodesy.equirectangular_},
613 @return: Closest point on the great circle line (L{LatLon}).
615 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
616 see function L{pygeodesy.equirectangular_}.
618 @raise NotImplementedError: Keyword argument C{B{within}=False}
619 is not (yet) supported.
621 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
623 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
625 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}
626 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
627 '''
628 # remove kwarg B{C{within}} if present
629 w = _xkwds_pop(wrap_adjust_limit, within=True)
630 if not w:
631 notImplemented(self, within=w)
633# # UNTESTED - handle C{B{within}=False} and C{B{within}=True}
634# wrap = _xkwds_get(options, wrap=False)
635# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap)
636# if fabs(a) < EPS or (within and a < EPS):
637# return point1
638# d = point1.distanceTo(point2, radius=radius, wrap=wrap)
639# if isnear0(d):
640# return point1 # or point2
641# elif fabs(d - a) < EPS or (a + EPS) > d:
642# return point2
643# f = a / d
644# if within:
645# if f > EPS1:
646# return point2
647# elif f < EPS:
648# return point1
649# return point1.intermediateTo(point2, f, wrap=wrap)
651 # without kwarg B{C{within}}, use backward compatible .nearestOn3
652 return self.nearestOn3([point1, point2], closed=False, radius=radius,
653 **wrap_adjust_limit)[0]
655 @deprecated_method
656 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover
657 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
659 @return: ... 2-Tuple C{(closest, distance)} of the closest
660 point (L{LatLon}) on the polygon and the distance
661 to that point from this point in C{meter}, same
662 units of B{C{radius}}.
663 '''
664 r = self.nearestOn3(points, closed=closed, radius=radius, **options)
665 return r.closest, r.distance
667 def nearestOn3(self, points, closed=False, radius=R_M, **wrap_adjust_limit):
668 '''Locate the point on a polygon closest to this point.
670 Distances are approximated by function L{pygeodesy.equirectangular_},
671 subject to the supplied B{C{options}}.
673 @arg points: The polygon points (L{LatLon}[]).
674 @kwarg closed: Optionally, close the polygon (C{bool}).
675 @kwarg radius: Mean earth radius (C{meter}).
676 @kwarg wrap_adjust_limit: Optional keyword arguments for function
677 L{sphericalTrigonometry.nearestOn3} and
678 L{pygeodesy.equirectangular_},
680 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
681 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_}
682 C{distance} between this and the C{closest} point converted to
683 C{meter}, same units as B{C{radius}}. The C{angle} from this
684 to the C{closest} point is in compass C{degrees360}, like
685 function L{pygeodesy.compassAngle}.
687 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
688 see function L{pygeodesy.equirectangular_}.
690 @raise PointsError: Insufficient number of B{C{points}}.
692 @raise TypeError: Some B{C{points}} are not C{LatLon}.
694 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
696 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_}
697 and L{pygeodesy.nearestOn5}.
698 '''
699 return nearestOn3(self, points, closed=closed, radius=radius,
700 LatLon=self.classof, **wrap_adjust_limit)
702 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
703 '''Convert this point to C{Karney}-based cartesian (ECEF)
704 coordinates.
706 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}}
707 and other keyword arguments, ignored
708 if C{B{Cartesian} is None}. Use
709 C{B{Cartesian}=...} to override
710 this L{Cartesian} class or specify
711 C{B{Cartesian}=None}.
713 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
714 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
715 C, M, datum)} with C{C} and C{M} if available.
717 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
718 '''
719 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
720 return LatLonSphericalBase.toCartesian(self, **kwds)
722 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
723 '''Compute the angles, sides and area of a spherical triangle.
725 @arg otherB: Second triangle point (C{LatLon}).
726 @arg otherC: Third triangle point (C{LatLon}).
727 @kwarg radius: Mean earth radius, ellipsoid or datum
728 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},
729 L{Datum} or L{a_f2Tuple}) or C{None}.
730 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
731 B{C{otherB}} and B{C{otherC}} points (C{bool}).
733 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if
734 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A,
735 a, B, b, C, c, D, E)}.
737 @see: Function L{triangle7} and U{Spherical trigonometry
738 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
739 '''
740 B = self.others(otherB=otherB)
741 C = self.others(otherC=otherC)
742 B, C, _ = _unrollon3(self, B, C, wrap)
744 r = self.philam + B.philam + C.philam
745 t = triangle8_(*r, wrap=wrap)
746 return self._xnamed(_t7Tuple(t, radius))
748 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
749 area=True, eps=EPS1, radius=R_M, wrap=False):
750 '''Trilaterate three points by I{area overlap} or I{perimeter
751 intersection} of three corresponding circles.
753 @arg distance1: Distance to this point (C{meter}, same units
754 as B{C{radius}}).
755 @arg point2: Second center point (C{LatLon}).
756 @arg distance2: Distance to point2 (C{meter}, same units as
757 B{C{radius}}).
758 @arg point3: Third center point (C{LatLon}).
759 @arg distance3: Distance to point3 (C{meter}, same units as
760 B{C{radius}}).
761 @kwarg area: If C{True} compute the area overlap, otherwise the
762 perimeter intersection of the circles (C{bool}).
763 @kwarg eps: The required I{minimal overlap} for C{B{area}=True}
764 or the I{intersection margin} for C{B{area}=False}
765 (C{meter}, same units as B{C{radius}}).
766 @kwarg radius: Mean earth radius (C{meter}, conventionally).
767 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
768 B{C{point2}} and B{C{point3}} (C{bool}).
770 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
771 with C{min} and C{max} in C{meter}, same units as B{C{eps}},
772 the corresponding trilaterated points C{minPoint} and
773 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number
774 of trilatered points found for the given B{C{eps}}.
776 If only a single trilaterated point is found, C{min I{is}
777 max}, C{minPoint I{is} maxPoint} and C{n = 1}.
779 For C{B{area}=True}, C{min} and C{max} are the smallest
780 respectively largest I{radial} overlap found.
782 For C{B{area}=False}, C{min} and C{max} represent the
783 nearest respectively farthest intersection margin.
785 If C{B{area}=True} and all 3 circles are concentric, C{n =
786 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}}
787 with the smallest B{C{distance#}} C{min} and C{max} the
788 largest B{C{distance#}}.
790 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
791 insufficient overlap for C{B{area}=True} or
792 no intersection or all (near-)concentric for
793 C{B{area}=False}.
795 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
797 @raise ValueError: Coincident B{C{point2}} or B{C{point3}} or invalid
798 B{C{distance1}}, B{C{distance2}}, B{C{distance3}}
799 or B{C{radius}}.
800 '''
801 return _trilaterate5(self, distance1,
802 self.others(point2=point2), distance2,
803 self.others(point3=point3), distance3,
804 area=area, radius=radius, eps=eps, wrap=wrap)
807_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
810def areaOf(points, radius=R_M, wrap=False): # was=True
811 '''Calculate the area of a (spherical) polygon or composite
812 (with the pointsjoined by great circle arcs).
814 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
815 or L{BooleanGH}).
816 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
817 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
818 or C{None}.
819 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{points}}
820 (C{bool}).
822 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
823 or C{radians} if B{C{radius}} is C{None}).
825 @raise PointsError: Insufficient number of B{C{points}}.
827 @raise TypeError: Some B{C{points}} are not L{LatLon}.
829 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
831 @note: The area is based on I{Karney}'s U{'Area of a spherical
832 polygon'<https://MathOverflow.net/questions/97711/
833 the-area-of-spherical-polygons>}, 3rd Answer.
835 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
836 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
837 L{ellipsoidalKarney.areaOf}.
839 @example:
841 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
842 >>> areaOf(b) # 8666058750.718977
844 >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1)
845 >>> areaOf(c) # 6.18e9
846 '''
847 if _MODS.booleans.isBoolean(points):
848 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
850 _at2, _t_2 = atan2, tan_2
851 _un, _w180 = unrollPI, wrap180
853 Ps = _T00.PointsIter(points, loop=1, wrap=wrap)
854 p1 = p2 = Ps[0]
855 a1, b1 = p1.philam
856 ta1, z1 = _t_2(a1), None
858 A = Fsum() # mean phi
859 R = Fsum() # see L{pygeodesy.excessKarney_}
860 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
861 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
862 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
863 D = Fsum()
864 for i, p2 in Ps.enumerate(closed=True):
865 a2, b2 = p2.philam
866 db, b2 = _un(b1, b2, wrap=wrap and not Ps.looped)
867 A += a2
868 ta2 = _t_2(a2)
869 tdb = _t_2(db, points=i)
870 R += _at2(tdb * (ta1 + ta2),
871 _1_0 + ta1 * ta2)
872 ta1, b1 = ta2, b2
874 if not p2.isequalTo(p1, eps=EPS):
875 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
876 if z1 is not None:
877 D += _w180(z - z1) # (z - z1 + 540) ...
878 D += _w180(z2 - z) # (z2 - z + 540) % 360 - 180
879 p1, z1 = p2, z2
881 R = abs(R * _2_0)
882 if abs(D) < _90_0: # ispolar(points)
883 R = abs(R - PI2)
884 if radius:
885 a = degrees(A.fover(len(A))) # mean lat
886 R *= _mean_radius(radius, a)**2
887 return float(R)
890def _destination2(a, b, r, t):
891 '''(INTERNAL) Destination lat- and longitude in C{radians}.
893 @arg a: Latitude (C{radians}).
894 @arg b: Longitude (C{radians}).
895 @arg r: Angular distance (C{radians}).
896 @arg t: Bearing (compass C{radians}).
898 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
899 '''
900 # see <https://www.EdWilliams.org/avform.htm#LL>
901 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
902 ca *= sr
904 a = asin1(ct * ca + cr * sa)
905 d = atan2(st * ca, cr - sa * sin(a))
906 # note, in EdWilliams.org/avform.htm W is + and E is -
907 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
910def _int3d2(s, end, wrap, _i_, Vector, hs):
911 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
912 # and similar logic in .ellipsoidalBaseDI._intersect3
913 a1, b1 = s.philam
915 if isscalar(end): # bearing, get pseudo-end point
916 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
917 else: # must be a point
918 s.others(end, name=_end_ + _i_)
919 hs.append(end.height)
920 a2, b2 = end.philam
921 if wrap:
922 a2, b2 = _Wrap.philam(a2, b2)
924 db, b2 = unrollPI(b1, b2, wrap=wrap)
925 if max(fabs(db), fabs(a2 - a1)) < EPS:
926 raise _ValueError(_SPACE_(_line_ + _i_, _null_))
927 # note, in EdWilliams.org/avform.htm W is + and E is -
928 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
929 -(b1 + b2) * _0_5)
930 cb21 *= sin(a1 - a2) # sa21
931 sb21 *= sin(a1 + a2) # sa12
932 x = Vector(sb12 * cb21 - cb12 * sb21,
933 cb12 * cb21 + sb12 * sb21,
934 cos(a1) * cos(a2) * sin(db)) # ll=start
935 return x.unit(), (db, (a2 - a1)) # negated d
938def _intdot(ds, a1, b1, a, b, wrap):
939 # compute dot product ds . (-b + b1, a - a1)
940 db, _ = unrollPI(b1, b, wrap=wrap)
941 return fdot(ds, db, a - a1)
944def intersecant2(center, circle, point, other, **radius_exact_height_wrap):
945 '''Compute the intersections of a circle and a (great circle) line given as
946 two points or as a point and bearing.
948 @arg center: Center of the circle (L{LatLon}).
949 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
950 or a point on the circle (L{LatLon}).
951 @arg point: A point on the (great circle) line (L{LatLon}).
952 @arg other: An other point on the (great circle) line (L{LatLon}) or
953 the bearing at the B{C{point}} (compass C{degrees360}).
954 @kwarg radius_exact_height_wrap: Optional keyword arguments, see
955 method L{LatLon.intersecant2} for further details.
957 @return: 2-Tuple of the intersection points (representing a chord), each
958 an instance of the B{C{point}} class. Both points are the same
959 instance if the (great circle) line is tangent to the circle.
961 @raise IntersectionError: The circle and line do not intersect.
963 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
964 B{C{circle}} or B{C{other}} invalid.
966 @raise UnitError: Invalid B{C{circle}}, B{C{other}}, B{C{radius}},
967 B{C{exact}}, B{C{height}} or B{C{napieradius}}.
968 '''
969 c = _T00.others(center=center)
970 p = _T00.others(point=point)
971 try:
972 return _intersecant2(c, circle, p, other, **radius_exact_height_wrap)
973 except (TypeError, ValueError) as x:
974 raise _xError(x, center=center, circle=circle, point=point, other=other,
975 **radius_exact_height_wrap)
978def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
979 LatLon=None, **LatLon_kwds):
980 # (INTERNAL) Intersect two (spherical) lines, see L{intersection}
981 # above, separated to allow callers to embellish any exceptions
983 s1, s2 = start1, start2
984 if wrap:
985 s2 = _Wrap.point(s2)
986 hs = [s1.height, s2.height]
988 a1, b1 = s1.philam
989 a2, b2 = s2.philam
990 db, b2 = unrollPI(b1, b2, wrap=wrap)
991 r12 = vincentys_(a2, a1, db)
992 if fabs(r12) < EPS: # [nearly] coincident points
993 a, b = favg(a1, a2), favg(b1, b2)
995 # see <https://www.EdWilliams.org/avform.htm#Intersection>
996 elif isscalar(end1) and isscalar(end2): # both bearings
997 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
999 x1, x2 = (sr12 * ca1), (sr12 * ca2)
1000 if isnear0(x1) or isnear0(x2):
1001 raise IntersectionError(_parallel_)
1002 # handle domain error for equivalent longitudes,
1003 # see also functions asin_safe and acos_safe at
1004 # <https://www.EdWilliams.org/avform.htm#Math>
1005 t12, t13 = acos1((sa2 - sa1 * cr12) / x1), radiansPI2(end1)
1006 t21, t23 = acos1((sa1 - sa2 * cr12) / x2), radiansPI2(end2)
1007 if sin(db) > 0:
1008 t21 = PI2 - t21
1009 else:
1010 t12 = PI2 - t12
1011 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
1012 wrapPI(t21 - t23)) # angle 1-2-3)
1013 if isnear0(sx1) and isnear0(sx2):
1014 raise IntersectionError(_infinite_)
1015 sx3 = sx1 * sx2
1016# XXX if sx3 < 0:
1017# XXX raise ValueError(_ambiguous_)
1018 x3 = acos1(cr12 * sx3 - cx2 * cx1)
1019 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
1021 a, b = _destination2(a1, b1, r13, t13)
1022 # like .ellipsoidalBaseDI,_intersect3, if this intersection
1023 # is "before" the first point, use the antipodal intersection
1024 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
1025 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1027 else: # end point(s) or bearing(s)
1028 _N_vector_ = _MODS.nvectorBase._N_vector_
1030 x1, d1 = _int3d2(s1, end1, wrap, _1_, _N_vector_, hs)
1031 x2, d2 = _int3d2(s2, end2, wrap, _2_, _N_vector_, hs)
1032 x = x1.cross(x2)
1033 if x.length < EPS: # [nearly] colinear or parallel lines
1034 raise IntersectionError(_colinear_)
1035 a, b = x.philam
1036 # choose intersection similar to sphericalNvector
1037 if not (_intdot(d1, a1, b1, a, b, wrap) *
1038 _intdot(d2, a2, b2, a, b, wrap)) > 0:
1039 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1041 h = fmean(hs) if height is None else Height(height)
1042 return _LL3Tuple(degrees90(a), degrees180(b), h,
1043 intersection, LatLon, LatLon_kwds)
1046def intersection(start1, end1, start2, end2, height=None, wrap=False,
1047 LatLon=LatLon, **LatLon_kwds):
1048 '''Compute the intersection point of two lines, each defined
1049 by two points or a start point and bearing from North.
1051 @arg start1: Start point of the first line (L{LatLon}).
1052 @arg end1: End point of the first line (L{LatLon}) or
1053 the initial bearing at the first start point
1054 (compass C{degrees360}).
1055 @arg start2: Start point of the second line (L{LatLon}).
1056 @arg end2: End point of the second line (L{LatLon}) or
1057 the initial bearing at the second start point
1058 (compass C{degrees360}).
1059 @kwarg height: Optional height for the intersection point,
1060 overriding the mean height (C{meter}).
1061 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1062 B{C{start2}} and both B{C{end*}} points (C{bool}).
1063 @kwarg LatLon: Optional class to return the intersection
1064 point (L{LatLon}) or C{None}.
1065 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1066 arguments, ignored if C{B{LatLon} is None}.
1068 @return: The intersection point as a (B{C{LatLon}}) or if
1069 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon,
1070 height)}. An alternate intersection point might
1071 be the L{antipode} to the returned result.
1073 @raise IntersectionError: Ambiguous or infinite intersection
1074 or colinear, parallel or otherwise
1075 non-intersecting lines.
1077 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}}
1078 or B{C{end2}} point not L{LatLon}.
1080 @raise ValueError: Invalid B{C{height}} or C{null} line.
1082 @example:
1084 >>> p = LatLon(51.8853, 0.2545)
1085 >>> s = LatLon(49.0034, 2.5735)
1086 >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E'
1087 '''
1088 s1 = _T00.others(start1=start1)
1089 s2 = _T00.others(start2=start2)
1090 try:
1091 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap,
1092 LatLon=LatLon, **LatLon_kwds)
1093 except (TypeError, ValueError) as x:
1094 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1097def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1098 height=None, wrap=False, # was=True
1099 LatLon=LatLon, **LatLon_kwds):
1100 '''Compute the intersection points of two circles each defined
1101 by a center point and a radius.
1103 @arg center1: Center of the first circle (L{LatLon}).
1104 @arg rad1: Radius of the first circle (C{meter} or C{radians},
1105 see B{C{radius}}).
1106 @arg center2: Center of the second circle (L{LatLon}).
1107 @arg rad2: Radius of the second circle (C{meter} or C{radians},
1108 see B{C{radius}}).
1109 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1110 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1111 @kwarg eps: Required overlap (C{meter} or C{radians}, see
1112 B{C{radius}}).
1113 @kwarg height: Optional height for the intersection points (C{meter},
1114 conventionally) or C{None} for the I{"radical height"}
1115 at the I{radical line} between both centers.
1116 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{center2}}
1117 (C{bool}).
1118 @kwarg LatLon: Optional class to return the intersection
1119 points (L{LatLon}) or C{None}.
1120 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1121 arguments, ignored if C{B{LatLon} is None}.
1123 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1124 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1125 lon, height)}. For abutting circles, both intersection
1126 points are the same instance, aka the I{radical center}.
1128 @raise IntersectionError: Concentric, antipodal, invalid or
1129 non-intersecting circles.
1131 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1133 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1134 B{C{eps}} or B{C{height}}.
1136 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1138 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1139 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1140 '''
1141 c1 = _T00.others(center1=center1)
1142 c2 = _T00.others(center2=center2)
1143 try:
1144 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1145 height=height, wrap=wrap,
1146 LatLon=LatLon, **LatLon_kwds)
1147 except (TypeError, ValueError) as x:
1148 raise _xError(x, center1=center1, rad1=rad1,
1149 center2=center2, rad2=rad2, wrap=wrap)
1152def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1153 height=None, too_d=None, wrap=False, # was=True
1154 LatLon=LatLon, **LatLon_kwds):
1155 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1156 # above, separated to allow callers to embellish any exceptions
1158 def _dest3(bearing, h):
1159 a, b = _destination2(a1, b1, r1, bearing)
1160 return _LL3Tuple(degrees90(a), degrees180(b), h,
1161 intersections2, LatLon, LatLon_kwds)
1163 a1, b1 = c1.philam
1164 a2, b2 = c2.philam
1165 if wrap:
1166 a2, b2 = _Wrap.philam(a2, b2)
1168 r1, r2, f = _rads3(rad1, rad2, radius)
1169 if f: # swapped radii, swap centers
1170 a1, a2 = a2, a1 # PYCHOK swap!
1171 b1, b2 = b2, b1 # PYCHOK swap!
1173 db, b2 = unrollPI(b1, b2, wrap=wrap)
1174 d = vincentys_(a2, a1, db) # radians
1175 if d < max(r1 - r2, EPS):
1176 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1178 r = eps if radius is None else (m2radians(
1179 eps, radius=radius) if eps else _0_0)
1180 if r < _0_0:
1181 raise _ValueError(eps=r)
1183 x = fsumf_(r1, r2, -d) # overlap
1184 if x > max(r, EPS):
1185 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1186 x = sd * sr1
1187 if isnear0(x):
1188 raise _ValueError(_invalid_)
1189 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1191 elif x < r: # PYCHOK no cover
1192 t = (d * radius) if too_d is None else too_d
1193 raise IntersectionError(_too_(_Fmt.distant(t)))
1195 if height is None: # "radical height"
1196 f = _radical2(d, r1, r2).ratio
1197 h = Height(favg(c1.height, c2.height, f=f))
1198 else:
1199 h = Height(height)
1201 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1202 if x < EPS4: # externally ...
1203 r = _dest3(b, h)
1204 elif x > _PI_EPS4: # internally ...
1205 r = _dest3(b + PI, h)
1206 else:
1207 return _dest3(b + x, h), _dest3(b - x, h)
1208 return r, r # ... abutting circles
1211@deprecated_function
1212def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1213 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1214 '''
1215 return ispolar(points, wrap=wrap)
1218def _LL3Tuple(lat, lon, height, where, LatLon, LatLon_kwds):
1219 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1220 '''
1221 n = where.__name__
1222 if LatLon is None:
1223 r = LatLon3Tuple(lat, lon, height, name=n)
1224 else:
1225 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1226 r = LatLon(lat, lon, **kwds)
1227 return r
1230def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):
1231 '''Compute the I{geographic} mean of several points.
1233 @arg points: Points to be averaged (L{LatLon}[]).
1234 @kwarg height: Optional height at mean point, overriding the mean
1235 height (C{meter}).
1236 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}}
1237 (C{bool}).
1238 @kwarg LatLon: Optional class to return the mean point (L{LatLon})
1239 or C{None}.
1240 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1241 arguments, ignored if C{B{LatLon} is None}.
1243 @return: The geographic mean and height (B{C{LatLon}}) or a
1244 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}}
1245 is C{None}.
1247 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1249 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1250 '''
1251 def _N_vs(ps, w):
1252 Ps = _T00.PointsIter(ps, wrap=w)
1253 for p in Ps.iterate(closed=False):
1254 yield p._N_vector
1256 m = _MODS.nvectorBase
1257 # geographic, vectorial mean
1258 n = m.sumOf(_N_vs(points, wrap), h=height, Vector=m.NvectorBase)
1259 lat, lon, h = n.latlonheight
1260 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1263@deprecated_function
1264def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1265 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1267 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1268 point (L{LatLon}) on the polygon and the C{distance}
1269 between the C{closest} and the given B{C{point}}. The
1270 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1271 lon)} if B{C{LatLon}} is C{None} ...
1272 '''
1273 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1274 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1275 ll = LatLon2Tuple(ll.lat, ll.lon)
1276 return ll, d
1279def nearestOn3(point, points, closed=False, radius=R_M, wrap=False, adjust=True,
1280 limit=9, **LatLon_and_kwds):
1281 '''Locate the point on a path or polygon closest to a reference point.
1283 Distances are I{approximated} using function L{pygeodesy.equirectangular_},
1284 subject to the supplied B{C{options}}.
1286 @arg point: The reference point (L{LatLon}).
1287 @arg points: The path or polygon points (L{LatLon}[]).
1288 @kwarg closed: Optionally, close the polygon (C{bool}).
1289 @kwarg radius: Mean earth radius (C{meter}).
1290 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1291 B{C{points}} (C{bool}).
1292 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}).
1293 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}),
1294 default C{9 degrees} is about C{1,000 Kmeter} (for mean
1295 spherical earth radius L{R_KM}).
1296 @kwarg LatLon: Optional class to return the closest point (L{LatLon})
1297 or C{None}.
1298 @kwarg options: Optional keyword arguments for function
1299 L{pygeodesy.equirectangular_}.
1301 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1302 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1303 lon, height)} if B{C{LatLon}} is C{None}. The C{distance}
1304 is the L{pygeodesy.equirectangular_} distance between the
1305 C{closest} and the given B{C{point}} converted to C{meter},
1306 same units as B{C{radius}}. The C{angle} from the given
1307 B{C{point}} to the C{closest} is in compass C{degrees360},
1308 like function L{pygeodesy.compassAngle}. The C{height} is
1309 the (interpolated) height at the C{closest} point.
1311 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1312 see function L{pygeodesy.equirectangular_}.
1314 @raise PointsError: Insufficient number of B{C{points}}.
1316 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1318 @raise ValueError: Invalid B{C{radius}}.
1320 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}.
1321 '''
1322 t = _nearestOn5(point, points, closed=closed, wrap=wrap,
1323 adjust=adjust, limit=limit)
1324 d = degrees2m(t.distance, radius=radius)
1325 h = t.height
1326 n = nearestOn3.__name__
1328 kwds = _xkwds(LatLon_and_kwds, height=h, name=n)
1329 LL = _xkwds_pop(kwds, LatLon=LatLon)
1330 r = LatLon3Tuple(t.lat, t.lon, h, name=n) if LL is None else \
1331 LL(t.lat, t.lon, **kwds)
1332 return NearestOn3Tuple(r, d, t.angle, name=n)
1335def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1336 '''Compute the perimeter of a (spherical) polygon or composite
1337 (with great circle arcs joining the points).
1339 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1340 or L{BooleanGH}).
1341 @kwarg closed: Optionally, close the polygon (C{bool}).
1342 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1343 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1344 B{C{points}} (C{bool}).
1346 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1347 or C{radians} if B{C{radius}} is C{None}).
1349 @raise PointsError: Insufficient number of B{C{points}}.
1351 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1353 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1354 C{B{points}} a composite.
1356 @note: Distances are based on function L{pygeodesy.vincentys_}.
1358 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf}
1359 and L{ellipsoidalKarney.perimeterOf}.
1360 '''
1361 def _rads(ps, c, w): # angular edge lengths in radians
1362 Ps = _T00.PointsIter(ps, loop=1, wrap=w)
1363 a1, b1 = Ps[0].philam
1364 for p in Ps.iterate(closed=c):
1365 a2, b2 = p.philam
1366 db, b2 = unrollPI(b1, b2, wrap=w and not (c and Ps.looped))
1367 yield vincentys_(a2, a1, db)
1368 a1, b1 = a2, b2
1370 if _MODS.booleans.isBoolean(points):
1371 if not closed:
1372 raise _ValueError(closed=closed, points=_composite_)
1373 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1374 else:
1375 r = fsum(_rads(points, closed, wrap), floats=True)
1376 return _radians2m(r, radius)
1379def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1380 excess=excessAbc_,
1381 wrap=False):
1382 '''Compute the angles, sides, and area of a (spherical) triangle.
1384 @arg latA: First corner latitude (C{degrees}).
1385 @arg lonA: First corner longitude (C{degrees}).
1386 @arg latB: Second corner latitude (C{degrees}).
1387 @arg lonB: Second corner longitude (C{degrees}).
1388 @arg latC: Third corner latitude (C{degrees}).
1389 @arg lonC: Third corner longitude (C{degrees}).
1390 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1391 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1392 or C{None}.
1393 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1394 L{excessGirard_} or L{excessLHuilier_}).
1395 @kwarg wrap: If C{True}, wrap and L{pygeodesy.unroll180}
1396 longitudes (C{bool}).
1398 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1399 spherical angles C{A}, C{B} and C{C}, angular sides
1400 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1401 in I{square} C{meter} or same units as B{C{radius}}
1402 I{squared} or if C{B{radius}=0} or C{None}, a
1403 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in
1404 C{radians} with the I{spherical excess} C{E} as the
1405 C{unit area} in C{radians}.
1406 '''
1407 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA),
1408 Phi_(latB=latB), Lam_(lonB=lonB),
1409 Phi_(latC=latC), Lam_(lonC=lonC),
1410 excess=excess, wrap=wrap)
1411 return _t7Tuple(t, radius)
1414def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_,
1415 wrap=False):
1416 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1417 excess} of a (spherical) triangle.
1419 @arg phiA: First corner latitude (C{radians}).
1420 @arg lamA: First corner longitude (C{radians}).
1421 @arg phiB: Second corner latitude (C{radians}).
1422 @arg lamB: Second corner longitude (C{radians}).
1423 @arg phiC: Third corner latitude (C{radians}).
1424 @arg lamC: Third corner longitude (C{radians}).
1425 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1426 L{excessGirard_} or L{excessLHuilier_}).
1427 @kwarg wrap: If C{True}, L{pygeodesy.unrollPI} the
1428 longitudinal deltas (C{bool}).
1430 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1431 spherical angles C{A}, C{B} and C{C}, angular sides
1432 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1433 I{spherical excess} C{E}, all in C{radians}.
1434 '''
1435 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1436 d, _ = unrollPI(lamB, lamC, wrap=w)
1437 a = vincentys_(phiC, phiB, d)
1438 return a, (phiB, lamB, phiC, lamC, phiA, lamA) # rotate A, B, C
1440 def _A_r(a, sa, ca, sb, cb, sc, cc):
1441 s = sb * sc
1442 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1443 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2_'s
1445 # notation: side C{a} is oposite to corner C{A}, etc.
1446 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1447 b, r = _a_r(wrap, *r)
1448 c, _ = _a_r(wrap, *r)
1450 A, r = _A_r(a, *sincos2_(a, b, c))
1451 B, r = _A_r(b, *r)
1452 C, _ = _A_r(c, *r)
1454 D = fsumf_(PI2, -a, -b, -c) # deficit aka defect
1455 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else (
1456 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else
1457 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b))))
1459 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1462def _t7Tuple(t, radius):
1463 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1464 '''
1465 if radius: # not in (None, _0_0)
1466 r = radius if isscalar(radius) else \
1467 _ellipsoidal_datum(radius).ellipsoid.Rmean
1468 A, B, C = map1(degrees, t.A, t.B, t.C)
1469 t = Triangle7Tuple(A, (r * t.a),
1470 B, (r * t.b),
1471 C, (r * t.c), t.E * r**2)
1472 return t
1475__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1476 areaOf, # functions
1477 intersecant2, intersection, intersections2, ispolar,
1478 isPoleEnclosedBy, # DEPRECATED, use ispolar
1479 meanOf,
1480 nearestOn2, nearestOn3,
1481 perimeterOf,
1482 sumOf, # XXX == vector3d.sumOf
1483 triangle7, triangle8_)
1485# **) MIT License
1486#
1487# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1488#
1489# Permission is hereby granted, free of charge, to any person obtaining a
1490# copy of this software and associated documentation files (the "Software"),
1491# to deal in the Software without restriction, including without limitation
1492# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1493# and/or sell copies of the Software, and to permit persons to whom the
1494# Software is furnished to do so, subject to the following conditions:
1495#
1496# The above copyright notice and this permission notice shall be included
1497# in all copies or substantial portions of the Software.
1498#
1499# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1500# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1501# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1502# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1503# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1504# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1505# OTHER DEALINGS IN THE SOFTWARE.