Coverage for pygeodesy/sphericalTrigonometry.py: 94%
391 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -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, _r2m, _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, degrees90, degrees180, degrees2m, \
50 m2radians, radiansPI2, sincos2_, tan_2, _unrollon, \
51 unrollPI, _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.07.01'
59_parallel_ = 'parallel'
61_PI_EPS4 = PI - EPS4
62if _PI_EPS4 >= PI:
63 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4)
66class Cartesian(CartesianSphericalBase):
67 '''Extended to convert geocentric, L{Cartesian} points to
68 spherical, geodetic L{LatLon}.
69 '''
71 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon
72 '''Convert this cartesian point to a C{spherical} geodetic point.
74 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
75 arguments. Use C{B{LatLon}=...} to override
76 this L{LatLon} class or specify C{B{LatLon}=None}.
78 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is C{None},
79 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
80 with C{C} and C{M} if available.
82 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
83 '''
84 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
85 return CartesianSphericalBase.toLatLon(self, **kwds)
88class LatLon(LatLonSphericalBase):
89 '''New point on spherical model earth model.
91 @example:
93 >>> p = LatLon(52.205, 0.119) # height=0
94 '''
96 def _ab1_ab2_db5(self, other, wrap):
97 '''(INTERNAL) Helper for several methods.
98 '''
99 a1, b1 = self.philam
100 a2, b2 = self.others(other, up=2).philam
101 if wrap:
102 a2, b2 = _Wrap.philam(a2, b2)
103 db, b2 = unrollPI(b1, b2, wrap=wrap)
104 else: # unrollPI shortcut
105 db = b2 - b1
106 return a1, b1, a2, b2, db
108 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
109 '''Compute the (signed) distance from the start to the closest
110 point on the great circle line defined by a start and an
111 end point.
113 That is, if a perpendicular is drawn from this point to the
114 great circle line, the along-track distance is the distance
115 from the start point to the point where the perpendicular
116 crosses the line.
118 @arg start: Start point of the great circle line (L{LatLon}).
119 @arg end: End point of the great circle line (L{LatLon}).
120 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
121 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
122 the B{C{start}} and B{C{end}} point (C{bool}).
124 @return: Distance along the great circle line (C{radians}
125 if C{B{radius} is None} or C{meter}, same units
126 as B{C{radius}}), positive if I{after} the
127 B{C{start}} toward the B{C{end}} point of the
128 line, I{negative} if before or C{0} if at the
129 B{C{start}} point.
131 @raise TypeError: Invalid B{C{start}} or B{C{end}} point.
133 @raise ValueError: Invalid B{C{radius}}.
135 @example:
137 >>> p = LatLon(53.2611, -0.7972)
139 >>> s = LatLon(53.3206, -1.7297)
140 >>> e = LatLon(53.1887, 0.1334)
141 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58
142 '''
143 r, x, b = self._a_x_b3(start, end, radius, wrap)
144 cx = cos(x)
145 return _0_0 if isnear0(cx) else \
146 _r2m(copysign0(acos1(cos(r) / cx), cos(b)), radius)
148 def _a_x_b3(self, start, end, radius, wrap):
149 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo.
150 '''
151 s = self.others(start=start)
152 e = self.others(end=end)
153 s, e, w = _unrollon3(self, s, e, wrap)
155 r = Radius_(radius)
156 r = s.distanceTo(self, r, wrap=w) / r
158 b = radians(s.initialBearingTo(self, wrap=w)
159 - s.initialBearingTo(e, wrap=w))
160 x = asin(sin(r) * sin(b))
161 return r, x, -b
163 @deprecated_method
164 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover
165 '''DEPRECATED, use method L{initialBearingTo}.
166 '''
167 return self.initialBearingTo(other, wrap=wrap, raiser=raiser)
169 def crossingParallels(self, other, lat, wrap=False):
170 '''Return the pair of meridians at which a great circle defined
171 by this and an other point crosses the given latitude.
173 @arg other: The other point defining great circle (L{LatLon}).
174 @arg lat: Latitude at the crossing (C{degrees}).
175 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
176 B{C{other}} point (C{bool}).
178 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
179 C{None} if the great circle doesn't reach B{C{lat}}.
180 '''
181 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
182 sa, ca, sa1, ca1, \
183 sa2, ca2, sdb, cdb = sincos2_(radians(lat), a1, a2, db)
184 sa1 *= ca2 * ca
186 x = sa1 * sdb
187 y = sa1 * cdb - ca1 * sa2 * ca
188 z = ca1 * sdb * ca2 * sa
190 h = hypot(x, y)
191 if h < EPS or fabs(z) > h: # PYCHOK no cover
192 return None # great circle doesn't reach latitude
194 m = atan2(-y, x) + b1 # longitude at max latitude
195 d = acos1(z / h) # delta longitude to intersections
196 return degrees180(m - d), degrees180(m + d)
198 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
199 '''Compute the (signed) distance from this point to
200 the great circle defined by a start and an end point.
202 @arg start: Start point of the great circle line (L{LatLon}).
203 @arg end: End point of the great circle line (L{LatLon}).
204 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
205 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
206 the B{C{start}} and B{C{end}} point (C{bool}).
208 @return: Distance to the great circle (C{radians} if
209 B{C{radius}} or C{meter}, same units as
210 B{C{radius}}), I{negative} if to the left or
211 I{positive} if to the right of the line.
213 @raise TypeError: If B{C{start}} or B{C{end}} is not L{LatLon}.
215 @raise ValueError: Invalid B{C{radius}}.
217 @example:
219 >>> p = LatLon(53.2611, -0.7972)
221 >>> s = LatLon(53.3206, -1.7297)
222 >>> e = LatLon(53.1887, 0.1334)
223 >>> d = p.crossTrackDistanceTo(s, e) # -307.5
224 '''
225 _, x, _ = self._a_x_b3(start, end, radius, wrap)
226 return _r2m(x, radius)
228 def destination(self, distance, bearing, radius=R_M, height=None):
229 '''Locate the destination from this point after having
230 travelled the given distance on the given initial bearing.
232 @arg distance: Distance travelled (C{meter}, same units as
233 B{C{radius}}).
234 @arg bearing: Bearing from this point (compass C{degrees360}).
235 @kwarg radius: Mean earth radius (C{meter}).
236 @kwarg height: Optional height at destination (C{meter}, same
237 units a B{C{radius}}).
239 @return: Destination point (L{LatLon}).
241 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
242 B{C{radius}} or B{C{height}}.
244 @example:
246 >>> p1 = LatLon(51.4778, -0.0015)
247 >>> p2 = p1.destination(7794, 300.7)
248 >>> p2.toStr() # '51.5135°N, 000.0983°W'
249 '''
250 a, b = self.philam
251 r, t = _angular(distance, radius, low=None), Bearing_(bearing)
253 a, b = _destination2(a, b, r, t)
254 h = self._heigHt(height)
255 return self.classof(degrees90(a), degrees180(b), height=h)
257 def distanceTo(self, other, radius=R_M, wrap=False):
258 '''Compute the (angular) distance from this to an other point.
260 @arg other: The other point (L{LatLon}).
261 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
262 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
263 the B{C{other}} point (C{bool}).
265 @return: Distance between this and the B{C{other}} point
266 (C{meter}, same units as B{C{radius}} or
267 C{radians} if B{C{radius}} is C{None}).
269 @raise TypeError: The B{C{other}} point is not L{LatLon}.
271 @raise ValueError: Invalid B{C{radius}}.
273 @example:
275 >>> p1 = LatLon(52.205, 0.119)
276 >>> p2 = LatLon(48.857, 2.351);
277 >>> d = p1.distanceTo(p2) # 404300
278 '''
279 a1, _, a2, _, db = self._ab1_ab2_db5(other, wrap)
280 return _r2m(vincentys_(a2, a1, db), radius)
282# @Property_RO
283# def Ecef(self):
284# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
285# '''
286# return _MODS.ecef.EcefKarney
288 def greatCircle(self, bearing, Vector=Vector3d, **Vector_kwds):
289 '''Compute the vector normal to great circle obtained by heading
290 on the given initial bearing from this point.
292 Direction of vector is such that initial bearing vector
293 b = c × n, where n is an n-vector representing this point.
295 @arg bearing: Bearing from this point (compass C{degrees360}).
296 @kwarg Vector: Vector class to return the great circle,
297 overriding the default L{Vector3d}.
298 @kwarg Vector_kwds: Optional, additional keyword argunents
299 for B{C{Vector}}.
301 @return: Vector representing great circle (C{Vector}).
303 @raise ValueError: Invalid B{C{bearing}}.
305 @example:
307 >>> p = LatLon(53.3206, -1.7297)
308 >>> g = p.greatCircle(96.0)
309 >>> g.toStr() # (-0.794, 0.129, 0.594)
310 '''
311 a, b = self.philam
312 sa, ca, sb, cb, st, ct = sincos2_(a, b, Bearing_(bearing))
314 return Vector(sb * ct - cb * sa * st,
315 -cb * ct - sb * sa * st,
316 ca * st, **Vector_kwds) # XXX .unit()?
318 def initialBearingTo(self, other, wrap=False, raiser=False):
319 '''Compute the initial bearing (forward azimuth) from this
320 to an other point.
322 @arg other: The other point (spherical L{LatLon}).
323 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
324 the B{C{other}} point (C{bool}).
325 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}),
326 use C{B{raiser}=True} for behavior like
327 C{sphericalNvector.LatLon.initialBearingTo}.
329 @return: Initial bearing (compass C{degrees360}).
331 @raise CrossError: If this and the B{C{other}} point coincide,
332 provided both B{C{raiser}} is C{True} and
333 L{pygeodesy.crosserrors} is C{True}.
335 @raise TypeError: The B{C{other}} point is not L{LatLon}.
337 @example:
339 >>> p1 = LatLon(52.205, 0.119)
340 >>> p2 = LatLon(48.857, 2.351)
341 >>> b = p1.initialBearingTo(p2) # 156.2
342 '''
343 a1, b1, a2, b2, db = self._ab1_ab2_db5(other, wrap)
344 # XXX behavior like sphericalNvector.LatLon.initialBearingTo
345 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(db)) < EPS:
346 raise CrossError(_point_, self, other=other, wrap=wrap, txt=_coincident_)
348 return degrees(bearing_(a1, b1, a2, b2, final=False))
350 def intermediateTo(self, other, fraction, height=None, wrap=False):
351 '''Locate the point at given fraction between (or along) this
352 and an other point.
354 @arg other: The other point (L{LatLon}).
355 @arg fraction: Fraction between both points (C{scalar},
356 0.0 at this and 1.0 at the other point).
357 @kwarg height: Optional height, overriding the intermediate
358 height (C{meter}).
359 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
360 B{C{other}} point (C{bool}).
362 @return: Intermediate point (L{LatLon}).
364 @raise TypeError: The B{C{other}} point is not L{LatLon}.
366 @raise ValueError: Invalid B{C{fraction}} or B{C{height}}.
368 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
370 @example:
372 >>> p1 = LatLon(52.205, 0.119)
373 >>> p2 = LatLon(48.857, 2.351)
374 >>> p = p1.intermediateTo(p2, 0.25) # 51.3721°N, 000.7073°E
375 '''
376 p = self
377 f = Scalar(fraction=fraction)
378 if not isnear0(f):
379 p = p.others(other)
380 if wrap:
381 p = _Wrap.point(p)
382 if not isnear1(f): # and not near0
383 a1, b1 = self.philam
384 a2, b2 = p.philam
385 db, b2 = unrollPI(b1, b2, wrap=wrap)
386 r = vincentys_(a2, a1, db)
387 sr = sin(r)
388 if isnon0(sr):
389 sa1, ca1, sa2, ca2, \
390 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2)
392 t = f * r
393 a = sin(r - t) # / sr superflous
394 b = sin( t) # / sr superflous
396 x = a * ca1 * cb1 + b * ca2 * cb2
397 y = a * ca1 * sb1 + b * ca2 * sb2
398 z = a * sa1 + b * sa2
400 a = atan2(z, hypot(x, y))
401 b = atan2(y, x)
403 else: # PYCHOK no cover
404 a = favg(a1, a2, f=f) # coincident
405 b = favg(b1, b2, f=f)
407 h = self._havg(other, f=f, h=height)
408 p = self.classof(degrees90(a), degrees180(b), height=h)
409 return p
411 def intersection(self, end1, other, end2, height=None, wrap=False):
412 '''Compute the intersection point of two lines, each defined by
413 two points or a start point and bearing from North.
415 @arg end1: End point of this line (L{LatLon}) or the initial
416 bearing at this point (compass C{degrees360}).
417 @arg other: Start point of the other line (L{LatLon}).
418 @arg end2: End point of the other line (L{LatLon}) or the
419 initial bearing at the B{C{other}} point (compass
420 C{degrees360}).
421 @kwarg height: Optional height for intersection point,
422 overriding the mean height (C{meter}).
423 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
424 B{C{start2}} and both B{C{end*}} points (C{bool}).
426 @return: The intersection point (L{LatLon}). An alternate
427 intersection point might be the L{antipode} to
428 the returned result.
430 @raise IntersectionError: Ambiguous or infinite intersection
431 or colinear, parallel or otherwise
432 non-intersecting lines.
434 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}}
435 or B{C{end2}} not C{scalar} nor L{LatLon}.
437 @raise ValueError: Invalid B{C{height}} or C{null} line.
439 @example:
441 >>> p = LatLon(51.8853, 0.2545)
442 >>> s = LatLon(49.0034, 2.5735)
443 >>> i = p.intersection(108.547, s, 32.435) # '50.9078°N, 004.5084°E'
444 '''
445 try:
446 s2 = self.others(other)
447 return _intersect(self, end1, s2, end2, height=height, wrap=wrap,
448 LatLon=self.classof)
449 except (TypeError, ValueError) as x:
450 raise _xError(x, start1=self, end1=end1,
451 other=other, end2=end2, wrap=wrap)
453 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0,
454 height=None, wrap=True):
455 '''Compute the intersection points of two circles, each defined
456 by a center point and radius.
458 @arg rad1: Radius of the this circle (C{meter} or C{radians},
459 see B{C{radius}}).
460 @arg other: Center point of the other circle (L{LatLon}).
461 @arg rad2: Radius of the other circle (C{meter} or C{radians},
462 see B{C{radius}}).
463 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
464 B{C{rad2}} and B{C{eps}} are given in C{radians}).
465 @kwarg eps: Required overlap (C{meter} or C{radians}, see
466 B{C{radius}}).
467 @kwarg height: Optional height for the intersection points (C{meter},
468 conventionally) or C{None} for the I{"radical height"}
469 at the I{radical line} between both centers.
470 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
471 B{C{other}} point (C{bool}).
473 @return: 2-Tuple of the intersection points, each a L{LatLon}
474 instance. For abutting circles, both intersection
475 points are the same instance, aka the I{radical center}.
477 @raise IntersectionError: Concentric, antipodal, invalid or
478 non-intersecting circles.
480 @raise TypeError: If B{C{other}} is not L{LatLon}.
482 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
483 B{C{eps}} or B{C{height}}.
484 '''
485 try:
486 c2 = self.others(other)
487 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps,
488 height=height, wrap=wrap,
489 LatLon=self.classof)
490 except (TypeError, ValueError) as x:
491 raise _xError(x, center=self, rad1=rad1,
492 other=other, rad2=rad2, wrap=wrap)
494 @deprecated_method
495 def isEnclosedBy(self, points): # PYCHOK no cover
496 '''DEPRECATED, use method C{isenclosedBy}.'''
497 return self.isenclosedBy(points)
499 def isenclosedBy(self, points, wrap=False):
500 '''Check whether a (convex) polygon or composite encloses this point.
502 @arg points: The polygon points or composite (L{LatLon}[],
503 L{BooleanFHP} or L{BooleanGH}).
504 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
505 B{C{points}} (C{bool}).
507 @return: C{True} if this point is inside the polygon or
508 composite, C{False} otherwise.
510 @raise PointsError: Insufficient number of B{C{points}}.
512 @raise TypeError: Some B{C{points}} are not L{LatLon}.
514 @raise ValueError: Invalid B{C{points}}, non-convex polygon.
516 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
517 and L{pygeodesy.ispolar} especially if the B{C{points}} may
518 enclose a pole or wrap around the earth I{longitudinally}.
519 '''
520 if _MODS.booleans.isBoolean(points):
521 return points._encloses(self.lat, self.lon, wrap=wrap)
523 Ps = self.PointsIter(points, loop=2, dedup=True, wrap=wrap)
524 n0 = self._N_vector
526 v2 = Ps[0]._N_vector
527 p1 = Ps[1]
528 v1 = p1._N_vector
529 # check whether this point on same side of all
530 # polygon edges (to the left or right depending
531 # on the anti-/clockwise polygon direction)
532 gc1 = v2.cross(v1)
533 t0 = gc1.angleTo(n0) > PI_2
534 s0 = None
535 # get great-circle vector for each edge
536 for i, p2 in Ps.enumerate(closed=True):
537 if wrap and not Ps.looped:
538 p2 = _unrollon(p1, p2)
539 p1 = p2
540 v2 = p2._N_vector
541 gc = v1.cross(v2)
542 t = gc.angleTo(n0) > PI_2
543 if t != t0: # different sides of edge i
544 return False # outside
546 # check for convex polygon: angle between
547 # gc vectors, signed by direction of n0
548 # (otherwise the test above is not reliable)
549 s = signOf(gc1.angleTo(gc, vSign=n0))
550 if s != s0:
551 if s0 is None:
552 s0 = s
553 else:
554 t = _Fmt.SQUARE(points=i)
555 raise _ValueError(t, p2, wrap=wrap, txt=_not_(_convex_))
556 gc1, v1 = gc, v2
558 return True # inside
560 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
561 '''Find the midpoint between this and an other point.
563 @arg other: The other point (L{LatLon}).
564 @kwarg height: Optional height for midpoint, overriding
565 the mean height (C{meter}).
566 @kwarg fraction: Midpoint location from this point (C{scalar}),
567 may be negative or greater than 1.0.
568 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
569 B{C{other}} point (C{bool}).
571 @return: Midpoint (L{LatLon}).
573 @raise TypeError: The B{C{other}} point is not L{LatLon}.
575 @raise ValueError: Invalid B{C{height}}.
577 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
579 @example:
581 >>> p1 = LatLon(52.205, 0.119)
582 >>> p2 = LatLon(48.857, 2.351)
583 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
584 '''
585 if fraction is _0_5:
586 # see <https://MathForum.org/library/drmath/view/51822.html>
587 a1, b, a2, _, db = self._ab1_ab2_db5(other, wrap)
588 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db)
590 x = ca2 * cdb + ca1
591 y = ca2 * sdb
593 a = atan2(sa1 + sa2, hypot(x, y))
594 b += atan2(y, x)
596 h = self._havg(other, h=height)
597 r = self.classof(degrees90(a), degrees180(b), height=h)
598 else:
599 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
600 return r
602 def nearestOn(self, point1, point2, radius=R_M, **wrap_adjust_limit):
603 '''Locate the point between two points closest to this point.
605 Distances are approximated by function L{pygeodesy.equirectangular_},
606 subject to the supplied B{C{options}}.
608 @arg point1: Start point (L{LatLon}).
609 @arg point2: End point (L{LatLon}).
610 @kwarg radius: Mean earth radius (C{meter}).
611 @kwarg wrap_adjust_limit: Optional keyword arguments for functions
612 L{sphericalTrigonometry.nearestOn3} and
613 L{pygeodesy.equirectangular_},
615 @return: Closest point on the great circle line (L{LatLon}).
617 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
618 see function L{pygeodesy.equirectangular_}.
620 @raise NotImplementedError: Keyword argument C{B{within}=False}
621 is not (yet) supported.
623 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
625 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
627 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}
628 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
629 '''
630 # remove kwarg B{C{within}} if present
631 w = _xkwds_pop(wrap_adjust_limit, within=True)
632 if not w:
633 notImplemented(self, within=w)
635# # UNTESTED - handle C{B{within}=False} and C{B{within}=True}
636# wrap = _xkwds_get(options, wrap=False)
637# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap)
638# if fabs(a) < EPS or (within and a < EPS):
639# return point1
640# d = point1.distanceTo(point2, radius=radius, wrap=wrap)
641# if isnear0(d):
642# return point1 # or point2
643# elif fabs(d - a) < EPS or (a + EPS) > d:
644# return point2
645# f = a / d
646# if within:
647# if f > EPS1:
648# return point2
649# elif f < EPS:
650# return point1
651# return point1.intermediateTo(point2, f, wrap=wrap)
653 # without kwarg B{C{within}}, use backward compatible .nearestOn3
654 return self.nearestOn3([point1, point2], closed=False, radius=radius,
655 **wrap_adjust_limit)[0]
657 @deprecated_method
658 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover
659 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
661 @return: ... 2-Tuple C{(closest, distance)} of the closest
662 point (L{LatLon}) on the polygon and the distance
663 to that point from this point in C{meter}, same
664 units of B{C{radius}}.
665 '''
666 r = self.nearestOn3(points, closed=closed, radius=radius, **options)
667 return r.closest, r.distance
669 def nearestOn3(self, points, closed=False, radius=R_M, **wrap_adjust_limit):
670 '''Locate the point on a polygon closest to this point.
672 Distances are approximated by function L{pygeodesy.equirectangular_},
673 subject to the supplied B{C{options}}.
675 @arg points: The polygon points (L{LatLon}[]).
676 @kwarg closed: Optionally, close the polygon (C{bool}).
677 @kwarg radius: Mean earth radius (C{meter}).
678 @kwarg wrap_adjust_limit: Optional keyword arguments for function
679 L{sphericalTrigonometry.nearestOn3} and
680 L{pygeodesy.equirectangular_},
682 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
683 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_}
684 C{distance} between this and the C{closest} point converted to
685 C{meter}, same units as B{C{radius}}. The C{angle} from this
686 to the C{closest} point is in compass C{degrees360}, like
687 function L{pygeodesy.compassAngle}.
689 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
690 see function L{pygeodesy.equirectangular_}.
692 @raise PointsError: Insufficient number of B{C{points}}.
694 @raise TypeError: Some B{C{points}} are not C{LatLon}.
696 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
698 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_}
699 and L{pygeodesy.nearestOn5}.
700 '''
701 return nearestOn3(self, points, closed=closed, radius=radius,
702 LatLon=self.classof, **wrap_adjust_limit)
704 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
705 '''Convert this point to C{Karney}-based cartesian (ECEF)
706 coordinates.
708 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}}
709 and other keyword arguments, ignored
710 if C{B{Cartesian} is None}. Use
711 C{B{Cartesian}=...} to override
712 this L{Cartesian} class or specify
713 C{B{Cartesian}=None}.
715 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
716 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
717 C, M, datum)} with C{C} and C{M} if available.
719 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
720 '''
721 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
722 return LatLonSphericalBase.toCartesian(self, **kwds)
724 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
725 '''Compute the angles, sides and area of a spherical triangle.
727 @arg otherB: Second triangle point (C{LatLon}).
728 @arg otherC: Third triangle point (C{LatLon}).
729 @kwarg radius: Mean earth radius, ellipsoid or datum
730 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},
731 L{Datum} or L{a_f2Tuple}) or C{None}.
732 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
733 B{C{otherB}} and B{C{otherC}} points (C{bool}).
735 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if
736 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A,
737 a, B, b, C, c, D, E)}.
739 @see: Function L{triangle7} and U{Spherical trigonometry
740 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
741 '''
742 B = self.others(otherB=otherB)
743 C = self.others(otherC=otherC)
744 B, C, _ = _unrollon3(self, B, C, wrap)
746 r = self.philam + B.philam + C.philam
747 t = triangle8_(*r, wrap=wrap)
748 return self._xnamed(_t7Tuple(t, radius))
750 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
751 area=True, eps=EPS1, radius=R_M, wrap=False):
752 '''Trilaterate three points by area overlap or perimeter intersection
753 of three corresponding circles.
755 @arg distance1: Distance to this point (C{meter}, same units
756 as B{C{radius}}).
757 @arg point2: Second center point (C{LatLon}).
758 @arg distance2: Distance to point2 (C{meter}, same units as
759 B{C{radius}}).
760 @arg point3: Third center point (C{LatLon}).
761 @arg distance3: Distance to point3 (C{meter}, same units as
762 B{C{radius}}).
763 @kwarg area: If C{True} compute the area overlap, otherwise the
764 perimeter intersection of the circles (C{bool}).
765 @kwarg eps: The required I{minimal overlap} for C{B{area}=True}
766 or the I{intersection margin} for C{B{area}=False}
767 (C{meter}, same units as B{C{radius}}).
768 @kwarg radius: Mean earth radius (C{meter}, conventionally).
769 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
770 B{C{point2}} and B{C{point3}} (C{bool}).
772 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
773 with C{min} and C{max} in C{meter}, same units as B{C{eps}},
774 the corresponding trilaterated points C{minPoint} and
775 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number
776 of trilatered points found for the given B{C{eps}}.
778 If only a single trilaterated point is found, C{min I{is}
779 max}, C{minPoint I{is} maxPoint} and C{n = 1}.
781 For C{B{area}=True}, C{min} and C{max} are the smallest
782 respectively largest I{radial} overlap found.
784 For C{B{area}=False}, C{min} and C{max} represent the
785 nearest respectively farthest intersection margin.
787 If C{B{area}=True} and all 3 circles are concentric, C{n =
788 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}}
789 with the smallest B{C{distance#}} C{min} and C{max} the
790 largest B{C{distance#}}.
792 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
793 insufficient overlap for C{B{area}=True} or
794 no intersection or all (near-)concentric for
795 C{B{area}=False}.
797 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
799 @raise ValueError: Coincident B{C{point2}} or B{C{point3}} or invalid
800 B{C{distance1}}, B{C{distance2}}, B{C{distance3}}
801 or B{C{radius}}.
802 '''
803 return _trilaterate5(self, distance1,
804 self.others(point2=point2), distance2,
805 self.others(point3=point3), distance3,
806 area=area, radius=radius, eps=eps, wrap=wrap)
809_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
812def areaOf(points, radius=R_M, wrap=False): # was=True
813 '''Calculate the area of a (spherical) polygon or composite
814 (with the pointsjoined by great circle arcs).
816 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
817 or L{BooleanGH}).
818 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
819 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
820 or C{None}.
821 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{points}}
822 (C{bool}).
824 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
825 or C{radians} if B{C{radius}} is C{None}).
827 @raise PointsError: Insufficient number of B{C{points}}.
829 @raise TypeError: Some B{C{points}} are not L{LatLon}.
831 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
833 @note: The area is based on I{Karney}'s U{'Area of a spherical
834 polygon'<https://MathOverflow.net/questions/97711/
835 the-area-of-spherical-polygons>}, 3rd Answer.
837 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
838 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
839 L{ellipsoidalKarney.areaOf}.
841 @example:
843 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
844 >>> areaOf(b) # 8666058750.718977
846 >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1)
847 >>> areaOf(c) # 6.18e9
848 '''
849 if _MODS.booleans.isBoolean(points):
850 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
852 _at2, _t_2 = atan2, tan_2
853 _un, _w180 = unrollPI, wrap180
855 Ps = _T00.PointsIter(points, loop=1, wrap=wrap)
856 p1 = p2 = Ps[0]
857 a1, b1 = p1.philam
858 ta1, z1 = _t_2(a1), None
860 A = Fsum() # mean phi
861 R = Fsum() # see L{pygeodesy.excessKarney_}
862 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
863 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
864 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
865 D = Fsum()
866 for i, p2 in Ps.enumerate(closed=True):
867 a2, b2 = p2.philam
868 db, b2 = _un(b1, b2, wrap=wrap and not Ps.looped)
869 A += a2
870 ta2 = _t_2(a2)
871 tdb = _t_2(db, points=i)
872 R += _at2(tdb * (ta1 + ta2),
873 _1_0 + ta1 * ta2)
874 ta1, b1 = ta2, b2
876 if not p2.isequalTo(p1, eps=EPS):
877 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
878 if z1 is not None:
879 D += _w180(z - z1) # (z - z1 + 540) ...
880 D += _w180(z2 - z) # (z2 - z + 540) % 360 - 180
881 p1, z1 = p2, z2
883 R = abs(R * _2_0)
884 if abs(D) < _90_0: # ispolar(points)
885 R = abs(R - PI2)
886 if radius:
887 a = degrees(A.fover(len(A))) # mean lat
888 R *= _mean_radius(radius, a)**2
889 return float(R)
892def _destination2(a, b, r, t):
893 '''(INTERNAL) Destination lat- and longitude in C{radians}.
895 @arg a: Latitude (C{radians}).
896 @arg b: Longitude (C{radians}).
897 @arg r: Angular distance (C{radians}).
898 @arg t: Bearing (compass C{radians}).
900 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
901 '''
902 # see <https://www.EdWilliams.org/avform.htm#LL>
903 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
904 ca *= sr
906 a = asin1(ct * ca + cr * sa)
907 d = atan2(st * ca, cr - sa * sin(a))
908 # note, in EdWilliams.org/avform.htm W is + and E is -
909 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
912def _int3d2(s, end, wrap, _i_, Vector, hs):
913 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
914 # and similar logic in .ellipsoidalBaseDI._intersect3
915 a1, b1 = s.philam
917 if isscalar(end): # bearing, get pseudo-end point
918 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
919 else: # must be a point
920 s.others(end, name=_end_ + _i_)
921 hs.append(end.height)
922 a2, b2 = end.philam
923 if wrap:
924 a2, b2 = _Wrap.philam(a2, b2)
926 db, b2 = unrollPI(b1, b2, wrap=wrap)
927 if max(fabs(db), fabs(a2 - a1)) < EPS:
928 raise _ValueError(_SPACE_(_line_ + _i_, _null_))
929 # note, in EdWilliams.org/avform.htm W is + and E is -
930 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
931 -(b1 + b2) * _0_5)
932 cb21 *= sin(a1 - a2) # sa21
933 sb21 *= sin(a1 + a2) # sa12
934 x = Vector(sb12 * cb21 - cb12 * sb21,
935 cb12 * cb21 + sb12 * sb21,
936 cos(a1) * cos(a2) * sin(db)) # ll=start
937 return x.unit(), (db, (a2 - a1)) # negated d
940def _intdot(ds, a1, b1, a, b, wrap):
941 # compute dot product ds . (-b + b1, a - a1)
942 db, _ = unrollPI(b1, b, wrap=wrap)
943 return fdot(ds, db, a - a1)
946def intersecant2(center, circle, point, bearing, radius=R_M, exact=False,
947 height=None, wrap=False): # was=True
948 '''Compute the intersections of a circle and a line.
950 @arg center: Center of the circle (L{LatLon}).
951 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
952 or a point on the circle (L{LatLon}).
953 @arg point: A point in- or outside the circle (L{LatLon}).
954 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or
955 a second point on the line (L{LatLon}).
956 @kwarg radius: Mean earth radius (C{meter}, conventionally).
957 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
958 destination and distance, if C{False} use the basic
959 rhumb methods (C{bool}) or if C{None} use the I{great
960 circle} methods.
961 @kwarg height: Optional height for the intersection points (C{meter},
962 conventionally) or C{None}.
963 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the B{C{point}}
964 and the B{C{circle}} and B{C{bearing}} points (C{bool}).
966 @return: 2-Tuple of the intersection points (representing a chord),
967 each an instance of this class. For a tangent line, both
968 points are the same instance, the B{C{point}} or wrapped
969 or I{normalized}.
971 @raise IntersectionError: The circle and line do not intersect.
973 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
974 B{C{circle}} or B{C{bearing}} invalid.
976 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
977 B{C{exact}} or B{C{height}}.
978 '''
979 c = _T00.others(center=center)
980 p = _T00.others(point=point)
981 try:
982 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact,
983 height=height, wrap=wrap)
984 except (TypeError, ValueError) as x:
985 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing, exact=exact)
988def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
989 LatLon=None, **LatLon_kwds):
990 # (INTERNAL) Intersect two (spherical) lines, see L{intersection}
991 # above, separated to allow callers to embellish any exceptions
993 s1, s2 = start1, start2
994 if wrap:
995 s2 = _Wrap.point(s2)
996 hs = [s1.height, s2.height]
998 a1, b1 = s1.philam
999 a2, b2 = s2.philam
1000 db, b2 = unrollPI(b1, b2, wrap=wrap)
1001 r12 = vincentys_(a2, a1, db)
1002 if fabs(r12) < EPS: # [nearly] coincident points
1003 a, b = favg(a1, a2), favg(b1, b2)
1005 # see <https://www.EdWilliams.org/avform.htm#Intersection>
1006 elif isscalar(end1) and isscalar(end2): # both bearings
1007 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
1009 x1, x2 = (sr12 * ca1), (sr12 * ca2)
1010 if isnear0(x1) or isnear0(x2):
1011 raise IntersectionError(_parallel_)
1012 # handle domain error for equivalent longitudes,
1013 # see also functions asin_safe and acos_safe at
1014 # <https://www.EdWilliams.org/avform.htm#Math>
1015 t1, t2 = acos1((sa2 - sa1 * cr12) / x1), \
1016 acos1((sa1 - sa2 * cr12) / x2)
1017 if sin(db) > 0:
1018 t12, t21 = t1, PI2 - t2
1019 else:
1020 t12, t21 = PI2 - t1, t2
1021 t13, t23 = radiansPI2(end1), radiansPI2(end2)
1022 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
1023 wrapPI(t21 - t23)) # angle 1-2-3)
1024 if isnear0(sx1) and isnear0(sx2):
1025 raise IntersectionError(_infinite_)
1026 sx3 = sx1 * sx2
1027# XXX if sx3 < 0:
1028# XXX raise ValueError(_ambiguous_)
1029 x3 = acos1(cr12 * sx3 - cx2 * cx1)
1030 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
1032 a, b = _destination2(a1, b1, r13, t13)
1033 # like .ellipsoidalBaseDI,_intersect3, if this intersection
1034 # is "before" the first point, use the antipodal intersection
1035 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
1036 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1038 else: # end point(s) or bearing(s)
1039 _N_vector_ = _MODS.nvectorBase._N_vector_
1041 x1, d1 = _int3d2(s1, end1, wrap, _1_, _N_vector_, hs)
1042 x2, d2 = _int3d2(s2, end2, wrap, _2_, _N_vector_, hs)
1043 x = x1.cross(x2)
1044 if x.length < EPS: # [nearly] colinear or parallel lines
1045 raise IntersectionError(_colinear_)
1046 a, b = x.philam
1047 # choose intersection similar to sphericalNvector
1048 if not (_intdot(d1, a1, b1, a, b, wrap) *
1049 _intdot(d2, a2, b2, a, b, wrap)) > 0:
1050 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1052 h = fmean(hs) if height is None else Height(height)
1053 return _LL3Tuple(degrees90(a), degrees180(b), h,
1054 intersection, LatLon, LatLon_kwds)
1057def intersection(start1, end1, start2, end2, height=None, wrap=False,
1058 LatLon=LatLon, **LatLon_kwds):
1059 '''Compute the intersection point of two lines, each defined
1060 by two points or a start point and bearing from North.
1062 @arg start1: Start point of the first line (L{LatLon}).
1063 @arg end1: End point of the first line (L{LatLon}) or
1064 the initial bearing at the first start point
1065 (compass C{degrees360}).
1066 @arg start2: Start point of the second line (L{LatLon}).
1067 @arg end2: End point of the second line (L{LatLon}) or
1068 the initial bearing at the second start point
1069 (compass C{degrees360}).
1070 @kwarg height: Optional height for the intersection point,
1071 overriding the mean height (C{meter}).
1072 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll
1073 B{C{start2}} and both B{C{end*}} points (C{bool}).
1074 @kwarg LatLon: Optional class to return the intersection
1075 point (L{LatLon}) or C{None}.
1076 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1077 arguments, ignored if C{B{LatLon} is None}.
1079 @return: The intersection point as a (B{C{LatLon}}) or if
1080 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon,
1081 height)}. An alternate intersection point might
1082 be the L{antipode} to the returned result.
1084 @raise IntersectionError: Ambiguous or infinite intersection
1085 or colinear, parallel or otherwise
1086 non-intersecting lines.
1088 @raise TypeError: A B{C{start1}}, B{C{end1}}, B{C{start2}}
1089 or B{C{end2}} point not L{LatLon}.
1091 @raise ValueError: Invalid B{C{height}} or C{null} line.
1093 @example:
1095 >>> p = LatLon(51.8853, 0.2545)
1096 >>> s = LatLon(49.0034, 2.5735)
1097 >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E'
1098 '''
1099 s1 = _T00.others(start1=start1)
1100 s2 = _T00.others(start2=start2)
1101 try:
1102 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap,
1103 LatLon=LatLon, **LatLon_kwds)
1104 except (TypeError, ValueError) as x:
1105 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1108def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1109 height=None, wrap=False, # was=True
1110 LatLon=LatLon, **LatLon_kwds):
1111 '''Compute the intersection points of two circles each defined
1112 by a center point and a radius.
1114 @arg center1: Center of the first circle (L{LatLon}).
1115 @arg rad1: Radius of the first circle (C{meter} or C{radians},
1116 see B{C{radius}}).
1117 @arg center2: Center of the second circle (L{LatLon}).
1118 @arg rad2: Radius of the second circle (C{meter} or C{radians},
1119 see B{C{radius}}).
1120 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1121 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1122 @kwarg eps: Required overlap (C{meter} or C{radians}, see
1123 B{C{radius}}).
1124 @kwarg height: Optional height for the intersection points (C{meter},
1125 conventionally) or C{None} for the I{"radical height"}
1126 at the I{radical line} between both centers.
1127 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{center2}}
1128 (C{bool}).
1129 @kwarg LatLon: Optional class to return the intersection
1130 points (L{LatLon}) or C{None}.
1131 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1132 arguments, ignored if C{B{LatLon} is None}.
1134 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1135 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1136 lon, height)}. For abutting circles, both intersection
1137 points are the same instance, aka the I{radical center}.
1139 @raise IntersectionError: Concentric, antipodal, invalid or
1140 non-intersecting circles.
1142 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1144 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1145 B{C{eps}} or B{C{height}}.
1147 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1149 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1150 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1151 '''
1152 c1 = _T00.others(center1=center1)
1153 c2 = _T00.others(center2=center2)
1154 try:
1155 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1156 height=height, wrap=wrap,
1157 LatLon=LatLon, **LatLon_kwds)
1158 except (TypeError, ValueError) as x:
1159 raise _xError(x, center1=center1, rad1=rad1,
1160 center2=center2, rad2=rad2, wrap=wrap)
1163def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1164 height=None, too_d=None, wrap=False, # was=True
1165 LatLon=LatLon, **LatLon_kwds):
1166 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1167 # above, separated to allow callers to embellish any exceptions
1169 def _dest3(bearing, h):
1170 a, b = _destination2(a1, b1, r1, bearing)
1171 return _LL3Tuple(degrees90(a), degrees180(b), h,
1172 intersections2, LatLon, LatLon_kwds)
1174 a1, b1 = c1.philam
1175 a2, b2 = c2.philam
1176 if wrap:
1177 a2, b2 = _Wrap.philam(a2, b2)
1179 r1, r2, f = _rads3(rad1, rad2, radius)
1180 if f: # swapped radii, swap centers
1181 a1, a2 = a2, a1 # PYCHOK swap!
1182 b1, b2 = b2, b1 # PYCHOK swap!
1184 db, b2 = unrollPI(b1, b2, wrap=wrap)
1185 d = vincentys_(a2, a1, db) # radians
1186 if d < max(r1 - r2, EPS):
1187 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1189 r = eps if radius is None else (m2radians(
1190 eps, radius=radius) if eps else _0_0)
1191 if r < _0_0:
1192 raise _ValueError(eps=r)
1194 x = fsumf_(r1, r2, -d) # overlap
1195 if x > max(r, EPS):
1196 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1197 x = sd * sr1
1198 if isnear0(x):
1199 raise _ValueError(_invalid_)
1200 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1202 elif x < r: # PYCHOK no cover
1203 t = (d * radius) if too_d is None else too_d
1204 raise IntersectionError(_too_(_Fmt.distant(t)))
1206 if height is None: # "radical height"
1207 f = _radical2(d, r1, r2).ratio
1208 h = Height(favg(c1.height, c2.height, f=f))
1209 else:
1210 h = Height(height)
1212 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1213 if x < EPS4: # externally ...
1214 r = _dest3(b, h)
1215 elif x > _PI_EPS4: # internally ...
1216 r = _dest3(b + PI, h)
1217 else:
1218 return _dest3(b + x, h), _dest3(b - x, h)
1219 return r, r # ... abutting circles
1222@deprecated_function
1223def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1224 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1225 '''
1226 return ispolar(points, wrap=wrap)
1229def _LL3Tuple(lat, lon, height, where, LatLon, LatLon_kwds):
1230 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1231 '''
1232 n = where.__name__
1233 if LatLon is None:
1234 r = LatLon3Tuple(lat, lon, height, name=n)
1235 else:
1236 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1237 r = LatLon(lat, lon, **kwds)
1238 return r
1241def meanOf(points, height=None, wrap=False, LatLon=LatLon, **LatLon_kwds):
1242 '''Compute the I{geographic} mean of several points.
1244 @arg points: Points to be averaged (L{LatLon}[]).
1245 @kwarg height: Optional height at mean point, overriding the mean
1246 height (C{meter}).
1247 @kwarg wrap: If C{True}, wrap or I{normalize} the B{C{points}}
1248 (C{bool}).
1249 @kwarg LatLon: Optional class to return the mean point (L{LatLon})
1250 or C{None}.
1251 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1252 arguments, ignored if C{B{LatLon} is None}.
1254 @return: The geographic mean and height (B{C{LatLon}}) or a
1255 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}}
1256 is C{None}.
1258 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1260 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1261 '''
1262 def _N_vs(ps, w):
1263 Ps = _T00.PointsIter(ps, wrap=w)
1264 for p in Ps.iterate(closed=False):
1265 yield p._N_vector
1267 m = _MODS.nvectorBase
1268 # geographic, vectorial mean
1269 n = m.sumOf(_N_vs(points, wrap), h=height, Vector=m.NvectorBase)
1270 lat, lon, h = n.latlonheight
1271 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1274@deprecated_function
1275def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1276 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1278 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1279 point (L{LatLon}) on the polygon and the C{distance}
1280 between the C{closest} and the given B{C{point}}. The
1281 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1282 lon)} if B{C{LatLon}} is C{None} ...
1283 '''
1284 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1285 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1286 ll = LatLon2Tuple(ll.lat, ll.lon)
1287 return ll, d
1290def nearestOn3(point, points, closed=False, radius=R_M, wrap=False, adjust=True,
1291 limit=9, **LatLon_and_kwds):
1292 '''Locate the point on a path or polygon closest to a reference point.
1294 Distances are I{approximated} using function L{pygeodesy.equirectangular_},
1295 subject to the supplied B{C{options}}.
1297 @arg point: The reference point (L{LatLon}).
1298 @arg points: The path or polygon points (L{LatLon}[]).
1299 @kwarg closed: Optionally, close the polygon (C{bool}).
1300 @kwarg radius: Mean earth radius (C{meter}).
1301 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1302 B{C{points}} (C{bool}).
1303 @kwarg adjust: See function L{pygeodesy.equirectangular_} (C{bool}).
1304 @kwarg limit: See function L{pygeodesy.equirectangular_} (C{degrees}),
1305 default C{9 degrees} is about C{1,000 Kmeter} (for mean
1306 spherical earth radius L{R_KM}).
1307 @kwarg LatLon: Optional class to return the closest point (L{LatLon})
1308 or C{None}.
1309 @kwarg options: Optional keyword arguments for function
1310 L{pygeodesy.equirectangular_}.
1312 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1313 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1314 lon, height)} if B{C{LatLon}} is C{None}. The C{distance}
1315 is the L{pygeodesy.equirectangular_} distance between the
1316 C{closest} and the given B{C{point}} converted to C{meter},
1317 same units as B{C{radius}}. The C{angle} from the given
1318 B{C{point}} to the C{closest} is in compass C{degrees360},
1319 like function L{pygeodesy.compassAngle}. The C{height} is
1320 the (interpolated) height at the C{closest} point.
1322 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1323 see function L{pygeodesy.equirectangular_}.
1325 @raise PointsError: Insufficient number of B{C{points}}.
1327 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1329 @raise ValueError: Invalid B{C{radius}}.
1331 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}.
1332 '''
1333 t = _nearestOn5(point, points, closed=closed, wrap=wrap,
1334 adjust=adjust, limit=limit)
1335 d = degrees2m(t.distance, radius=radius)
1336 h = t.height
1337 n = nearestOn3.__name__
1339 kwds = _xkwds(LatLon_and_kwds, height=h, name=n)
1340 LL = _xkwds_pop(kwds, LatLon=LatLon)
1341 r = LatLon3Tuple(t.lat, t.lon, h, name=n) if LL is None else \
1342 LL(t.lat, t.lon, **kwds)
1343 return NearestOn3Tuple(r, d, t.angle, name=n)
1346def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1347 '''Compute the perimeter of a (spherical) polygon or composite
1348 (with great circle arcs joining the points).
1350 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1351 or L{BooleanGH}).
1352 @kwarg closed: Optionally, close the polygon (C{bool}).
1353 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1354 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll the
1355 B{C{points}} (C{bool}).
1357 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1358 or C{radians} if B{C{radius}} is C{None}).
1360 @raise PointsError: Insufficient number of B{C{points}}.
1362 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1364 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1365 C{B{points}} a composite.
1367 @note: Distances are based on function L{pygeodesy.vincentys_}.
1369 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf}
1370 and L{ellipsoidalKarney.perimeterOf}.
1371 '''
1372 def _rads(ps, c, w): # angular edge lengths in radians
1373 Ps = _T00.PointsIter(ps, loop=1, wrap=w)
1374 a1, b1 = Ps[0].philam
1375 for p in Ps.iterate(closed=c):
1376 a2, b2 = p.philam
1377 db, b2 = unrollPI(b1, b2, wrap=w and not (c and Ps.looped))
1378 yield vincentys_(a2, a1, db)
1379 a1, b1 = a2, b2
1381 if _MODS.booleans.isBoolean(points):
1382 if not closed:
1383 raise _ValueError(closed=closed, points=_composite_)
1384 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1385 else:
1386 r = fsum(_rads(points, closed, wrap), floats=True)
1387 return _r2m(r, radius)
1390def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1391 excess=excessAbc_,
1392 wrap=False):
1393 '''Compute the angles, sides, and area of a (spherical) triangle.
1395 @arg latA: First corner latitude (C{degrees}).
1396 @arg lonA: First corner longitude (C{degrees}).
1397 @arg latB: Second corner latitude (C{degrees}).
1398 @arg lonB: Second corner longitude (C{degrees}).
1399 @arg latC: Third corner latitude (C{degrees}).
1400 @arg lonC: Third corner longitude (C{degrees}).
1401 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1402 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1403 or C{None}.
1404 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1405 L{excessGirard_} or L{excessLHuilier_}).
1406 @kwarg wrap: If C{True}, wrap and L{pygeodesy.unroll180}
1407 longitudes (C{bool}).
1409 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1410 spherical angles C{A}, C{B} and C{C}, angular sides
1411 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1412 in I{square} C{meter} or same units as B{C{radius}}
1413 I{squared} or if C{B{radius}=0} or C{None}, a
1414 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in
1415 C{radians} with the I{spherical excess} C{E} as the
1416 C{unit area} in C{radians}.
1417 '''
1418 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA),
1419 Phi_(latB=latB), Lam_(lonB=lonB),
1420 Phi_(latC=latC), Lam_(lonC=lonC),
1421 excess=excess, wrap=wrap)
1422 return _t7Tuple(t, radius)
1425def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_,
1426 wrap=False):
1427 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1428 excess} of a (spherical) triangle.
1430 @arg phiA: First corner latitude (C{radians}).
1431 @arg lamA: First corner longitude (C{radians}).
1432 @arg phiB: Second corner latitude (C{radians}).
1433 @arg lamB: Second corner longitude (C{radians}).
1434 @arg phiC: Third corner latitude (C{radians}).
1435 @arg lamC: Third corner longitude (C{radians}).
1436 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1437 L{excessGirard_} or L{excessLHuilier_}).
1438 @kwarg wrap: If C{True}, L{pygeodesy.unrollPI} the
1439 longitudinal deltas (C{bool}).
1441 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1442 spherical angles C{A}, C{B} and C{C}, angular sides
1443 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1444 I{spherical excess} C{E}, all in C{radians}.
1445 '''
1446 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1447 d, _ = unrollPI(lamB, lamC, wrap=w)
1448 a = vincentys_(phiC, phiB, d)
1449 return a, (phiB, lamB, phiC, lamC, phiA, lamA)
1451 def _A_r(a, sa, ca, sb, cb, sc, cc):
1452 s = sb * sc
1453 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1454 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2's
1456 # notation: side C{a} is oposite to corner C{A}, etc.
1457 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1458 b, r = _a_r(wrap, *r)
1459 c, _ = _a_r(wrap, *r)
1461 A, r = _A_r(a, *sincos2_(a, b, c))
1462 B, r = _A_r(b, *r)
1463 C, _ = _A_r(c, *r)
1465 D = fsumf_(PI2, -a, -b, -c) # deficit aka defect
1466 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else (
1467 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else
1468 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b))))
1470 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1473def _t7Tuple(t, radius):
1474 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1475 '''
1476 if radius: # not in (None, _0_0)
1477 r = radius if isscalar(radius) else \
1478 _ellipsoidal_datum(radius).ellipsoid.Rmean
1479 A, B, C = map1(degrees, t.A, t.B, t.C)
1480 t = Triangle7Tuple(A, (r * t.a),
1481 B, (r * t.b),
1482 C, (r * t.c), t.E * r**2)
1483 return t
1486__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1487 areaOf, # functions
1488 intersecant2, intersection, intersections2, ispolar,
1489 isPoleEnclosedBy, # DEPRECATED, use ispolar
1490 meanOf,
1491 nearestOn2, nearestOn3,
1492 perimeterOf,
1493 sumOf, # XXX == vector3d.sumOf
1494 triangle7, triangle8_)
1496# **) MIT License
1497#
1498# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1499#
1500# Permission is hereby granted, free of charge, to any person obtaining a
1501# copy of this software and associated documentation files (the "Software"),
1502# to deal in the Software without restriction, including without limitation
1503# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1504# and/or sell copies of the Software, and to permit persons to whom the
1505# Software is furnished to do so, subject to the following conditions:
1506#
1507# The above copyright notice and this permission notice shall be included
1508# in all copies or substantial portions of the Software.
1509#
1510# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1511# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1512# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1513# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1514# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1515# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1516# OTHER DEALINGS IN THE SOFTWARE.