Coverage for pygeodesy/sphericalTrigonometry.py: 96%
374 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-03-31 10:52 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-03-31 10:52 -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, fsum_
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_, _LatLon_, _near_, _not_, _null_, \
35 _points_, _SPACE_, _too_
36from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS, _ALL_OTHER
37# from pygeodesy.named import notImplemented # from .points
38from pygeodesy.namedTuples import LatLon2Tuple, LatLon3Tuple, \
39 NearestOn3Tuple, Triangle7Tuple, \
40 Triangle8Tuple
41from pygeodesy.points import ispolar, nearestOn5 as _nearestOn5, \
42 notImplemented, Fmt as _Fmt # XXX shadowed
43from pygeodesy.props import deprecated_function, deprecated_method
44from pygeodesy.sphericalBase import _angular, CartesianSphericalBase, _intersecant2, \
45 LatLonSphericalBase, _rads3, _trilaterate5
46# from pygeodesy.streprs import Fmt as _Fmt # from .points XXX shadowed
47from pygeodesy.units import Bearing_, Height, Lam_, Phi_, Radius, \
48 Radius_, Scalar
49from pygeodesy.utily import acos1, asin1, degrees90, degrees180, degrees2m, \
50 m2radians, radiansPI2, sincos2_, tan_2, unrollPI, \
51 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.03.30'
59_parallel_ = 'parallel'
60_path_ = 'path'
62_PI_EPS4 = PI - EPS4
63if _PI_EPS4 >= PI:
64 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4)
67class Cartesian(CartesianSphericalBase):
68 '''Extended to convert geocentric, L{Cartesian} points to
69 spherical, geodetic L{LatLon}.
70 '''
72 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon
73 '''Convert this cartesian point to a C{spherical} geodetic point.
75 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
76 arguments. Use C{B{LatLon}=...} to override
77 this L{LatLon} class or specify C{B{LatLon}=None}.
79 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is C{None},
80 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
81 with C{C} and C{M} if available.
83 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
84 '''
85 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
86 return CartesianSphericalBase.toLatLon(self, **kwds)
89class LatLon(LatLonSphericalBase):
90 '''New point on spherical model earth model.
92 @example:
94 >>> p = LatLon(52.205, 0.119) # height=0
95 '''
97 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
98 '''Compute the (angular) distance (signed) from the start to
99 the closest point on the great circle path defined by a
100 start and an end point.
102 That is, if a perpendicular is drawn from this point to the
103 great circle path, the along-track distance is the distance
104 from the start point to the point where the perpendicular
105 crosses the path.
107 @arg start: Start point of great circle path (L{LatLon}).
108 @arg end: End point of great circle path (L{LatLon}).
109 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
110 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
112 @return: Distance along the great circle path (C{meter},
113 same units as B{C{radius}}) or C{radians} if
114 C{B{radius} is None}, positive if after the B{C{start}}
115 toward the B{C{end}} point of the path or negative
116 if before the B{C{start}} point.
118 @raise TypeError: Invalid B{C{start}} or B{C{end}} point.
120 @raise ValueError: Invalid B{C{radius}}.
122 @example:
124 >>> p = LatLon(53.2611, -0.7972)
126 >>> s = LatLon(53.3206, -1.7297)
127 >>> e = LatLon(53.1887, 0.1334)
128 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58
129 '''
130 r, x, b = self._trackDistanceTo3(start, end, radius, wrap)
131 cx = cos(x)
132 return _0_0 if isnear0(cx) else \
133 _r2m(copysign0(acos1(cos(r) / cx), cos(b)), radius)
135 @deprecated_method
136 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover
137 '''DEPRECATED, use method L{initialBearingTo}.
138 '''
139 return self.initialBearingTo(other, wrap=wrap, raiser=raiser)
141 def crossingParallels(self, other, lat, wrap=False):
142 '''Return the pair of meridians at which a great circle defined
143 by this and an other point crosses the given latitude.
145 @arg other: The other point defining great circle (L{LatLon}).
146 @arg lat: Latitude at the crossing (C{degrees}).
147 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
149 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
150 C{None} if the great circle doesn't reach B{C{lat}}.
151 '''
152 self.others(other)
154 a1, b1 = self.philam
155 a2, b2 = other.philam
157 a = radians(lat)
158 db, b2 = unrollPI(b1, b2, wrap=wrap)
160 sa, ca, sa1, ca1, \
161 sa2, ca2, sdb, cdb = sincos2_(a, a1, a2, db)
162 sa1 *= ca2 * ca
164 x = sa1 * sdb
165 y = sa1 * cdb - ca1 * sa2 * ca
166 z = ca1 * sdb * ca2 * sa
168 h = hypot(x, y)
169 if h < EPS or fabs(z) > h: # PYCHOK no cover
170 return None # great circle doesn't reach latitude
172 m = atan2(-y, x) + b1 # longitude at max latitude
173 d = acos1(z / h) # delta longitude to intersections
174 return degrees180(m - d), degrees180(m + d)
176 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
177 '''Compute the (angular) distance (signed) from this point to
178 the great circle defined by a start and an end point.
180 @arg start: Start point of great circle path (L{LatLon}).
181 @arg end: End point of great circle path (L{LatLon}).
182 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
183 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
185 @return: Distance to great circle (negative if to the left
186 or positive if to the right of the path) (C{meter},
187 same units as B{C{radius}} or C{radians} if
188 B{C{radius}} is C{None}).
190 @raise TypeError: The B{C{start}} or B{C{end}} point is not L{LatLon}.
192 @raise ValueError: Invalid B{C{radius}}.
194 @example:
196 >>> p = LatLon(53.2611, -0.7972)
198 >>> s = LatLon(53.3206, -1.7297)
199 >>> e = LatLon(53.1887, 0.1334)
200 >>> d = p.crossTrackDistanceTo(s, e) # -307.5
201 '''
202 _, x, _ = self._trackDistanceTo3(start, end, radius, wrap)
203 return _r2m(x, radius)
205 def destination(self, distance, bearing, radius=R_M, height=None):
206 '''Locate the destination from this point after having
207 travelled the given distance on the given initial bearing.
209 @arg distance: Distance travelled (C{meter}, same units as
210 B{C{radius}}).
211 @arg bearing: Bearing from this point (compass C{degrees360}).
212 @kwarg radius: Mean earth radius (C{meter}).
213 @kwarg height: Optional height at destination (C{meter}, same
214 units a B{C{radius}}).
216 @return: Destination point (L{LatLon}).
218 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
219 B{C{radius}} or B{C{height}}.
221 @example:
223 >>> p1 = LatLon(51.4778, -0.0015)
224 >>> p2 = p1.destination(7794, 300.7)
225 >>> p2.toStr() # '51.5135°N, 000.0983°W'
226 '''
227 a, b = self.philam
228 r, t = _angular(distance, radius), Bearing_(bearing)
230 a, b = _destination2(a, b, r, t)
231 h = self.height if height is None else Height(height)
232 return self.classof(degrees90(a), degrees180(b), height=h)
234 def distanceTo(self, other, radius=R_M, wrap=False):
235 '''Compute the (angular) distance from this to an other point.
237 @arg other: The other point (L{LatLon}).
238 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
239 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
241 @return: Distance between this and the B{C{other}} point
242 (C{meter}, same units as B{C{radius}} or
243 C{radians} if B{C{radius}} is C{None}).
245 @raise TypeError: The B{C{other}} point is not L{LatLon}.
247 @raise ValueError: Invalid B{C{radius}}.
249 @example:
251 >>> p1 = LatLon(52.205, 0.119)
252 >>> p2 = LatLon(48.857, 2.351);
253 >>> d = p1.distanceTo(p2) # 404300
254 '''
255 self.others(other)
257 a1, b1 = self.philam
258 a2, b2 = other.philam
260 db, _ = unrollPI(b1, b2, wrap=wrap)
261 return _r2m(vincentys_(a2, a1, db), radius)
263# @Property_RO
264# def Ecef(self):
265# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
266# '''
267# return _MODS.ecef.EcefKarney
269 def greatCircle(self, bearing):
270 '''Compute the vector normal to great circle obtained by heading
271 on the given initial bearing from this point.
273 Direction of vector is such that initial bearing vector
274 b = c × n, where n is an n-vector representing this point.
276 @arg bearing: Bearing from this point (compass C{degrees360}).
278 @return: Vector representing great circle (L{Vector3d}).
280 @raise ValueError: Invalid B{C{bearing}}.
282 @example:
284 >>> p = LatLon(53.3206, -1.7297)
285 >>> g = p.greatCircle(96.0)
286 >>> g.toStr() # (-0.794, 0.129, 0.594)
287 '''
288 a, b = self.philam
290 t = Bearing_(bearing)
291 sa, ca, sb, cb, st, ct = sincos2_(a, b, t)
293 return Vector3d(sb * ct - cb * sa * st,
294 -cb * ct - sb * sa * st,
295 ca * st) # XXX .unit()?
297 def initialBearingTo(self, other, wrap=False, raiser=False):
298 '''Compute the initial bearing (forward azimuth) from this
299 to an other point.
301 @arg other: The other point (spherical L{LatLon}).
302 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
303 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}),
304 use C{B{raiser}=True} for behavior like
305 C{sphericalNvector.LatLon.initialBearingTo}.
307 @return: Initial bearing (compass C{degrees360}).
309 @raise CrossError: If this and the B{C{other}} point coincide,
310 provided both B{C{raiser}} is C{True} and
311 L{pygeodesy.crosserrors} is C{True}.
313 @raise TypeError: The B{C{other}} point is not L{LatLon}.
315 @example:
317 >>> p1 = LatLon(52.205, 0.119)
318 >>> p2 = LatLon(48.857, 2.351)
319 >>> b = p1.initialBearingTo(p2) # 156.2
320 '''
321 self.others(other)
323 a1, b1 = self.philam
324 a2, b2 = other.philam
326 # XXX behavior like sphericalNvector.LatLon.initialBearingTo
327 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(b2 - b1)) < EPS:
328 raise CrossError(_points_, self, txt=_coincident_)
330 return degrees(bearing_(a1, b1, a2, b2, final=False, wrap=wrap))
332 def intermediateTo(self, other, fraction, height=None, wrap=False):
333 '''Locate the point at given fraction between (or along) this
334 and an other point.
336 @arg other: The other point (L{LatLon}).
337 @arg fraction: Fraction between both points (C{scalar},
338 0.0 at this and 1.0 at the other point).
339 @kwarg height: Optional height, overriding the intermediate
340 height (C{meter}).
341 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
343 @return: Intermediate point (L{LatLon}).
345 @raise TypeError: The B{C{other}} point is not L{LatLon}.
347 @raise ValueError: Invalid B{C{fraction}} or B{C{height}}.
349 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
351 @example:
353 >>> p1 = LatLon(52.205, 0.119)
354 >>> p2 = LatLon(48.857, 2.351)
355 >>> p = p1.intermediateTo(p2, 0.25) # 51.3721°N, 000.7073°E
356 '''
357 f = Scalar(fraction=fraction) # not high=_1_0
358 if isnear0(f): # PYCHOK no cover
359 r = self
360 elif isnear1(f) and not wrap: # PYCHOK no cover
361 r = self.others(other)
362 else:
363 a1, b1 = self.philam
364 a2, b2 = self.others(other).philam
366 db, b2 = unrollPI(b1, b2, wrap=wrap)
367 r = vincentys_(a2, a1, db)
368 sr = sin(r)
369 if isnon0(sr):
370 sa1, ca1, sa2, ca2, \
371 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2)
373 t = f * r
374 a = sin(r - t) # / sr superflous
375 b = sin( t) # / sr superflous
377 x = a * ca1 * cb1 + b * ca2 * cb2
378 y = a * ca1 * sb1 + b * ca2 * sb2
379 z = a * sa1 + b * sa2
381 a = atan2(z, hypot(x, y))
382 b = atan2(y, x)
384 else: # PYCHOK no cover
385 a = favg(a1, a2, f=f) # coincident
386 b = favg(b1, b2, f=f)
388 h = self._havg(other, f=f) if height is None else Height(height)
389 r = self.classof(degrees90(a), degrees180(b), height=h)
390 return r
392 def intersection(self, end1, other, end2, height=None, wrap=False):
393 '''Compute the intersection point of two paths, each defined
394 by two points or a start point and bearing from North.
396 @arg end1: End point of this path (L{LatLon}) or the
397 initial bearing at this point (compass
398 C{degrees360}).
399 @arg other: Start point of the other path (L{LatLon}).
400 @arg end2: End point of the other path (L{LatLon}) or
401 the initial bearing at the other start point
402 (compass C{degrees360}).
403 @kwarg height: Optional height for intersection point,
404 overriding the mean height (C{meter}).
405 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
407 @return: The intersection point (L{LatLon}). An alternate
408 intersection point might be the L{antipode} to
409 the returned result.
411 @raise IntersectionError: Ambiguous or infinite intersection
412 or colinear, parallel or otherwise
413 non-intersecting paths.
415 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}}
416 or B{C{end2}} not C{scalar} nor L{LatLon}.
418 @raise ValueError: Invalid B{C{height}} or C{null} path.
420 @example:
422 >>> p = LatLon(51.8853, 0.2545)
423 >>> s = LatLon(49.0034, 2.5735)
424 >>> i = p.intersection(108.547, s, 32.435) # '50.9078°N, 004.5084°E'
425 '''
426 try:
427 s2 = self.others(other)
428 return _intersect(self, end1, s2, end2, height=height, wrap=wrap,
429 LatLon=self.classof)
430 except (TypeError, ValueError) as x:
431 raise _xError(x, start1=self, end1=end1, other=other, end2=end2)
433 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0,
434 height=None, wrap=True):
435 '''Compute the intersection points of two circles, each defined
436 by a center point and radius.
438 @arg rad1: Radius of the this circle (C{meter} or C{radians},
439 see B{C{radius}}).
440 @arg other: Center point of the other circle (L{LatLon}).
441 @arg rad2: Radius of the other circle (C{meter} or C{radians},
442 see B{C{radius}}).
443 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
444 B{C{rad2}} and B{C{eps}} are given in C{radians}).
445 @kwarg eps: Required overlap (C{meter} or C{radians}, see
446 B{C{radius}}).
447 @kwarg height: Optional height for the intersection points (C{meter},
448 conventionally) or C{None} for the I{"radical height"}
449 at the I{radical line} between both centers.
450 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
452 @return: 2-Tuple of the intersection points, each a L{LatLon}
453 instance. For abutting circles, both intersection
454 points are the same instance, aka I{radical center}.
456 @raise IntersectionError: Concentric, antipodal, invalid or
457 non-intersecting circles.
459 @raise TypeError: If B{C{other}} is not L{LatLon}.
461 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
462 B{C{eps}} or B{C{height}}.
463 '''
464 try:
465 c2 = self.others(other)
466 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps,
467 height=height, wrap=wrap,
468 LatLon=self.classof)
469 except (TypeError, ValueError) as x:
470 raise _xError(x, center=self, rad1=rad1, other=other, rad2=rad2)
472 def isenclosedBy(self, points):
473 '''Check whether a (convex) polygon encloses this point.
475 @arg points: The polygon points (L{LatLon}[]).
477 @return: C{True} if the polygon encloses this point,
478 C{False} otherwise.
480 @raise PointsError: Insufficient number of B{C{points}}.
482 @raise TypeError: Some B{C{points}} are not L{LatLon}.
484 @raise ValueError: Invalid B{C{points}}, non-convex polygon.
486 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
487 and L{pygeodesy.ispolar} especially if the B{C{points}} may
488 enclose a pole or wrap around the earth longitudinally.
490 @example:
492 >>> b = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1)
493 >>> p = LatLon(45,1, 1.1)
494 >>> inside = p.isEnclosedBy(b) # True
495 '''
496 Ps = self.PointsIter(points, loop=2, dedup=True)
497 n0 = self._N_vector
499 v2 = Ps[0]._N_vector
500 v1 = Ps[1]._N_vector
501 gc1 = v2.cross(v1)
502 # check whether this point on same side of all
503 # polygon edges (to the left or right depending
504 # on anti-/clockwise polygon direction)
505 t0 = gc1.angleTo(n0) > PI_2
506 s0 = None
507 # get great-circle vector for each edge
508 for i, p in Ps.enumerate(closed=True):
509 v2 = p._N_vector
510 gc = v1.cross(v2)
511 v1 = v2
513 t = gc.angleTo(n0) > PI_2
514 if t != t0: # different sides of edge i
515 return False # outside
517 # check for convex polygon: angle between
518 # gc vectors, signed by direction of n0
519 # (otherwise the test above is not reliable)
520 s = signOf(gc1.angleTo(gc, vSign=n0))
521 if s != s0:
522 if s0 is None:
523 s0 = s
524 else:
525 t = _Fmt.SQUARE(points=i)
526 raise _ValueError(t, p, txt=_not_(_convex_))
527 gc1 = gc
529 return True # inside
531 @deprecated_method
532 def isEnclosedBy(self, points): # PYCHOK no cover
533 '''DEPRECATED, use method C{isenclosedBy}.'''
534 return self.isenclosedBy(points)
536 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
537 '''Find the midpoint between this and an other point.
539 @arg other: The other point (L{LatLon}).
540 @kwarg height: Optional height for midpoint, overriding
541 the mean height (C{meter}).
542 @kwarg fraction: Midpoint location from this point (C{scalar}),
543 may be negative or greater than 1.0.
544 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
546 @return: Midpoint (L{LatLon}).
548 @raise TypeError: The B{C{other}} point is not L{LatLon}.
550 @raise ValueError: Invalid B{C{height}}.
552 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
554 @example:
556 >>> p1 = LatLon(52.205, 0.119)
557 >>> p2 = LatLon(48.857, 2.351)
558 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
559 '''
560 if fraction is _0_5:
561 self.others(other)
563 # see <https://MathForum.org/library/drmath/view/51822.html>
564 a1, b1 = self.philam
565 a2, b2 = other.philam
567 db, b2 = unrollPI(b1, b2, wrap=wrap)
569 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db)
571 x = ca2 * cdb + ca1
572 y = ca2 * sdb
574 a = atan2(sa1 + sa2, hypot(x, y))
575 b = atan2(y, x) + b1
577 h = self._havg(other) if height is None else Height(height)
578 r = self.classof(degrees90(a), degrees180(b), height=h)
579 else:
580 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
581 return r
583 def nearestOn(self, point1, point2, radius=R_M, **options):
584 '''Locate the point between two points closest to this point.
586 Distances are approximated by function L{pygeodesy.equirectangular_},
587 subject to the supplied B{C{options}}.
589 @arg point1: Start point (L{LatLon}).
590 @arg point2: End point (L{LatLon}).
591 @kwarg radius: Mean earth radius (C{meter}).
592 @kwarg options: Optional keyword arguments for function
593 L{pygeodesy.equirectangular_}.
595 @return: Closest point on the great circle path (L{LatLon}).
597 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
598 see function L{pygeodesy.equirectangular_}.
600 @raise NotImplementedError: Keyword argument C{B{within}=False}
601 is not (yet) supported.
603 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
605 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
607 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}
608 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
609 '''
610 # remove kwarg B{C{within}} if present
611 within = _xkwds_pop(options, within=True)
612 if not within:
613 notImplemented(self, within=within)
615# # UNTESTED - handle C{B{within}=False} and C{B{within}=True}
616# wrap = _xkwds_get(options, wrap=False)
617# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap)
618# if fabs(a) < EPS or (within and a < EPS):
619# return point1
620# d = point1.distanceTo(point2, radius=radius, wrap=wrap)
621# if isnear0(d):
622# return point1 # or point2
623# elif fabs(d - a) < EPS or (a + EPS) > d:
624# return point2
625# f = a / d
626# if within:
627# if f > EPS1:
628# return point2
629# elif f < EPS:
630# return point1
631# return point1.intermediateTo(point2, f, wrap=wrap)
633 # without kwarg B{C{within}}, use backward compatible .nearestOn3
634 return self.nearestOn3([point1, point2], closed=False, radius=radius,
635 **options)[0]
637 @deprecated_method
638 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover
639 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
641 @return: ... 2-Tuple C{(closest, distance)} of the closest
642 point (L{LatLon}) on the polygon and the distance
643 to that point from this point in C{meter}, same
644 units of B{C{radius}}.
645 '''
646 r = self.nearestOn3(points, closed=closed, radius=radius, **options)
647 return r.closest, r.distance
649 def nearestOn3(self, points, closed=False, radius=R_M, **options):
650 '''Locate the point on a polygon closest to this point.
652 Distances are approximated by function L{pygeodesy.equirectangular_},
653 subject to the supplied B{C{options}}.
655 @arg points: The polygon points (L{LatLon}[]).
656 @kwarg closed: Optionally, close the polygon (C{bool}).
657 @kwarg radius: Mean earth radius (C{meter}).
658 @kwarg options: Optional keyword arguments for function
659 L{pygeodesy.equirectangular_}.
661 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
662 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_}
663 C{distance} between this and the C{closest} point converted to
664 C{meter}, same units as B{C{radius}}. The C{angle} from this
665 to the C{closest} point is in compass C{degrees360}, like
666 function L{pygeodesy.compassAngle}.
668 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
669 see function L{pygeodesy.equirectangular_}.
671 @raise PointsError: Insufficient number of B{C{points}}.
673 @raise TypeError: Some B{C{points}} are not C{LatLon}.
675 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
677 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_}
678 and L{pygeodesy.nearestOn5}.
679 '''
680 if _LatLon_ in options:
681 raise _ValueError(**options)
683 lat, lon, d, c, h = _nearestOn5(self, points, closed=closed, **options)
684 return NearestOn3Tuple(self.classof(lat, lon, height=h),
685 degrees2m(d, radius=radius), c)
687 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
688 '''Convert this point to C{Karney}-based cartesian (ECEF)
689 coordinates.
691 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}}
692 and other keyword arguments, ignored
693 if C{B{Cartesian} is None}. Use
694 C{B{Cartesian}=...} to override
695 this L{Cartesian} class or specify
696 C{B{Cartesian}=None}.
698 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
699 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
700 C, M, datum)} with C{C} and C{M} if available.
702 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
703 '''
704 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
705 return LatLonSphericalBase.toCartesian(self, **kwds)
707 def _trackDistanceTo3(self, start, end, radius, wrap):
708 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo.
709 '''
710 self.others(start=start)
711 self.others(end=end)
713 r = Radius_(radius)
714 r = start.distanceTo(self, r, wrap=wrap) / r
716 b = radians(start.initialBearingTo(self, wrap=wrap))
717 e = radians(start.initialBearingTo(end, wrap=wrap))
718 x = asin(sin(r) * sin(b - e))
719 return r, x, (e - b)
721 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
722 '''Compute the angles, sides and area of a spherical triangle.
724 @arg otherB: Second triangle point (C{LatLon}).
725 @arg otherC: Third triangle point (C{LatLon}).
726 @kwarg radius: Mean earth radius, ellipsoid or datum
727 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},
728 L{Datum} or L{a_f2Tuple}) or C{None}.
729 @kwarg wrap: Wrap/unroll angular distances (C{bool}).
731 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if
732 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A,
733 a, B, b, C, c, D, E)}.
735 @see: Function L{triangle7} and U{Spherical trigonometry
736 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
737 '''
738 self.others(otherB=otherB)
739 self.others(otherC=otherC)
741 r = self.philam + otherB.philam + otherC.philam
742 t = triangle8_(*r, wrap=wrap)
743 return self._xnamed(_t7Tuple(t, radius))
745 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
746 area=True, eps=EPS1, radius=R_M, wrap=False):
747 '''Trilaterate three points by area overlap or perimeter intersection
748 of three corresponding circles.
750 @arg distance1: Distance to this point (C{meter}, same units
751 as B{C{radius}}).
752 @arg point2: Second center point (C{LatLon}).
753 @arg distance2: Distance to point2 (C{meter}, same units as
754 B{C{radius}}).
755 @arg point3: Third center point (C{LatLon}).
756 @arg distance3: Distance to point3 (C{meter}, same units as
757 B{C{radius}}).
758 @kwarg area: If C{True} compute the area overlap, otherwise the
759 perimeter intersection of the circles (C{bool}).
760 @kwarg eps: The required I{minimal overlap} for C{B{area}=True}
761 or the I{intersection margin} for C{B{area}=False}
762 (C{meter}, same units as B{C{radius}}).
763 @kwarg radius: Mean earth radius (C{meter}, conventionally).
764 @kwarg wrap: Wrap/unroll angular distances (C{bool}).
766 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
767 with C{min} and C{max} in C{meter}, same units as B{C{eps}},
768 the corresponding trilaterated points C{minPoint} and
769 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number
770 of trilatered points found for the given B{C{eps}}.
772 If only a single trilaterated point is found, C{min I{is}
773 max}, C{minPoint I{is} maxPoint} and C{n = 1}.
775 For C{B{area}=True}, C{min} and C{max} are the smallest
776 respectively largest I{radial} overlap found.
778 For C{B{area}=False}, C{min} and C{max} represent the
779 nearest respectively farthest intersection margin.
781 If C{B{area}=True} and all 3 circles are concentric, C{n =
782 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}}
783 with the smallest B{C{distance#}} C{min} and C{max} the
784 largest B{C{distance#}}.
786 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
787 insufficient overlap for C{B{area}=True} or
788 no intersection or all (near-)concentric for
789 C{B{area}=False}.
791 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
793 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}},
794 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.
795 '''
796 return _trilaterate5(self, distance1,
797 self.others(point2=point2), distance2,
798 self.others(point3=point3), distance3,
799 area=area, radius=radius, eps=eps, wrap=wrap)
802_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
805def areaOf(points, radius=R_M, wrap=True):
806 '''Calculate the area of a (spherical) polygon or composite
807 (with the pointsjoined by great circle arcs).
809 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
810 or L{BooleanGH}).
811 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
812 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
813 or C{None}.
814 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
816 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
817 or C{radians} if B{C{radius}} is C{None}).
819 @raise PointsError: Insufficient number of B{C{points}}.
821 @raise TypeError: Some B{C{points}} are not L{LatLon}.
823 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
825 @note: The area is based on I{Karney}'s U{'Area of a spherical
826 polygon'<https://MathOverflow.net/questions/97711/
827 the-area-of-spherical-polygons>}, 3rd Answer.
829 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
830 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
831 L{ellipsoidalKarney.areaOf}.
833 @example:
835 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
836 >>> areaOf(b) # 8666058750.718977
838 >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1)
839 >>> areaOf(c) # 6.18e9
840 '''
841 if _MODS.booleans.isBoolean(points):
842 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
844 Ps = _T00.PointsIter(points, loop=1)
845 p1 = p2 = Ps[0]
846 a1, b1 = p1.philam
847 ta1, z1 = tan_2(a1), None
849 A = Fsum() # mean phi
850 E = Fsum() # see L{pygeodesy.excessKarney_}
851 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
852 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
853 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
854 D = Fsum()
855 for i, p2 in Ps.enumerate(closed=True):
856 a2, b2 = p2.philam
857 db, b2 = unrollPI(b1, b2, wrap=wrap if i else False)
858 ta2 = tan_2(a2)
859 A += a2
860 E += atan2(tan_2(db, points=i) * (ta1 + ta2),
861 _1_0 + ta1 * ta2)
862 ta1, b1 = ta2, b2
864 if not p2.isequalTo(p1, eps=EPS):
865 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
866 if z1 is not None:
867 D += wrap180(z - z1) # (z - z1 + 540) ...
868 D += wrap180(z2 - z) # (z2 - z + 540) % 360 - 180
869 p1, z1 = p2, z2
871 R = fabs(E * _2_0)
872 if fabs(D) < _90_0: # ispolar(points)
873 R = fabs(R - PI2)
874 if radius:
875 a = degrees(A.fover(len(A))) # mean lat
876 R *= _mean_radius(radius, a)**2
877 return float(R)
880def _destination2(a, b, r, t):
881 '''(INTERNAL) Destination lat- and longitude in C{radians}.
883 @arg a: Latitude (C{radians}).
884 @arg b: Longitude (C{radians}).
885 @arg r: Angular distance (C{radians}).
886 @arg t: Bearing (compass C{radians}).
888 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
889 '''
890 # see <https://www.EdWilliams.org/avform.htm#LL>
891 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
892 ca *= sr
894 a = asin1(ct * ca + cr * sa)
895 d = atan2(st * ca, cr - sa * sin(a))
896 # note, in EdWilliams.org/avform.htm W is + and E is -
897 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
900def _int3d2(start, end, wrap, _i_, Vector, hs):
901 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
902 # and similar logic in .ellipsoidalBaseDI._intersect3
903 a1, b1 = start.philam
905 if isscalar(end): # bearing, get pseudo-end point
906 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
907 else: # must be a point
908 start.others(end, name=_end_ + _i_)
909 hs.append(end.height)
910 a2, b2 = end.philam
912 db, b2 = unrollPI(b1, b2, wrap=wrap)
913 if max(fabs(db), fabs(a2 - a1)) < EPS:
914 raise _ValueError(_SPACE_(_path_ + _i_, _null_))
915 # note, in EdWilliams.org/avform.htm W is + and E is -
916 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
917 -(b1 + b2) * _0_5)
918 cb21 *= sin(a1 - a2) # sa21
919 sb21 *= sin(a1 + a2) # sa12
920 x = Vector(sb12 * cb21 - cb12 * sb21,
921 cb12 * cb21 + sb12 * sb21,
922 cos(a1) * cos(a2) * sin(db)) # ll=start
923 return x.unit(), (db, (a2 - a1)) # negated d
926def _intdot(ds, a1, b1, a, b, wrap):
927 # compute dot product ds . (-b + b1, a - a1)
928 db, _ = unrollPI(b1, b, wrap=wrap)
929 return fdot(ds, db, a - a1)
932def intersecant2(center, circle, point, bearing, radius=R_M, exact=False,
933 height=None, wrap=True):
934 '''Compute the intersections of a circle and a line.
936 @arg center: Center of the circle (L{LatLon}).
937 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
938 or a point on the circle (L{LatLon}).
939 @arg point: A point in- or outside the circle (L{LatLon}).
940 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or
941 a second point on the line (L{LatLon}).
942 @kwarg radius: Mean earth radius (C{meter}, conventionally).
943 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
944 destination and distance, if C{False} use the basic
945 rhumb methods (C{bool}) or if C{None} use the I{great
946 circle} methods.
947 @kwarg height: Optional height for the intersection points (C{meter},
948 conventionally) or C{None}.
949 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
951 @return: 2-Tuple of the intersection points (representing a chord),
952 each an instance of this class. For a tangent line, each
953 point C{is} this very instance.
955 @raise IntersectionError: The circle and line do not intersect.
957 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
958 B{C{circle}} or B{C{bearing}} invalid.
960 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
961 B{C{exact}} or B{C{height}}.
962 '''
963 c = _T00.others(center=center)
964 p = _T00.others(point=point)
965 try:
966 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact,
967 height=height, wrap=wrap)
968 except (TypeError, ValueError) as x:
969 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing, exact=exact)
972def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
973 LatLon=None, **LatLon_kwds):
974 # (INTERNAL) Intersect two (spherical) path, see L{intersection}
975 # above, separated to allow callers to embellish any exceptions
977 hs = [start1.height, start2.height]
979 a1, b1 = start1.philam
980 a2, b2 = start2.philam
982 db, b2 = unrollPI(b1, b2, wrap=wrap)
983 r12 = vincentys_(a2, a1, db)
984 if fabs(r12) < EPS: # [nearly] coincident points
985 a, b = favg(a1, a2), favg(b1, b2)
987 # see <https://www.EdWilliams.org/avform.htm#Intersection>
988 elif isscalar(end1) and isscalar(end2): # both bearings
989 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
991 x1, x2 = (sr12 * ca1), (sr12 * ca2)
992 if isnear0(x1) or isnear0(x2):
993 raise IntersectionError(_parallel_)
994 # handle domain error for equivalent longitudes,
995 # see also functions asin_safe and acos_safe at
996 # <https://www.EdWilliams.org/avform.htm#Math>
997 t1, t2 = acos1((sa2 - sa1 * cr12) / x1), \
998 acos1((sa1 - sa2 * cr12) / x2)
999 if sin(db) > 0:
1000 t12, t21 = t1, PI2 - t2
1001 else:
1002 t12, t21 = PI2 - t1, t2
1003 t13, t23 = radiansPI2(end1), radiansPI2(end2)
1004 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
1005 wrapPI(t21 - t23)) # angle 1-2-3)
1006 if isnear0(sx1) and isnear0(sx2):
1007 raise IntersectionError(_infinite_)
1008 sx3 = sx1 * sx2
1009# XXX if sx3 < 0:
1010# XXX raise ValueError(_ambiguous_)
1011 x3 = acos1(cr12 * sx3 - cx2 * cx1)
1012 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
1014 a, b = _destination2(a1, b1, r13, t13)
1015 # like .ellipsoidalBaseDI,_intersect3, if this intersection
1016 # is "before" the first point, use the antipodal intersection
1017 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
1018 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1020 else: # end point(s) or bearing(s)
1021 _N_vector_ = _MODS.nvectorBase._N_vector_
1023 x1, d1 = _int3d2(start1, end1, wrap, _1_, _N_vector_, hs)
1024 x2, d2 = _int3d2(start2, end2, wrap, _2_, _N_vector_, hs)
1025 x = x1.cross(x2)
1026 if x.length < EPS: # [nearly] colinear or parallel paths
1027 raise IntersectionError(_colinear_)
1028 a, b = x.philam
1029 # choose intersection similar to sphericalNvector
1030 if not (_intdot(d1, a1, b1, a, b, wrap) *
1031 _intdot(d2, a2, b2, a, b, wrap)) > 0:
1032 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1034 h = fmean(hs) if height is None else Height(height)
1035 return _LL3Tuple(degrees90(a), degrees180(b), h,
1036 intersection, LatLon, LatLon_kwds)
1039def intersection(start1, end1, start2, end2, height=None, wrap=False,
1040 LatLon=LatLon, **LatLon_kwds):
1041 '''Compute the intersection point of two paths, each defined
1042 by two points or a start point and bearing from North.
1044 @arg start1: Start point of the first path (L{LatLon}).
1045 @arg end1: End point of the first path (L{LatLon}) or
1046 the initial bearing at the first start point
1047 (compass C{degrees360}).
1048 @arg start2: Start point of the second path (L{LatLon}).
1049 @arg end2: End point of the second path (L{LatLon}) or
1050 the initial bearing at the second start point
1051 (compass C{degrees360}).
1052 @kwarg height: Optional height for the intersection point,
1053 overriding the mean height (C{meter}).
1054 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1055 @kwarg LatLon: Optional class to return the intersection
1056 point (L{LatLon}) or C{None}.
1057 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1058 arguments, ignored if C{B{LatLon} is None}.
1060 @return: The intersection point as a (B{C{LatLon}}) or if
1061 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon,
1062 height)}. An alternate intersection point might
1063 be the L{antipode} to the returned result.
1065 @raise IntersectionError: Ambiguous or infinite intersection
1066 or colinear, parallel or otherwise
1067 non-intersecting paths.
1069 @raise TypeError: A B{C{start}} or B{C{end}} point not L{LatLon}.
1071 @raise ValueError: Invalid B{C{height}} or C{null} path.
1073 @example:
1075 >>> p = LatLon(51.8853, 0.2545)
1076 >>> s = LatLon(49.0034, 2.5735)
1077 >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E'
1078 '''
1079 s1 = _T00.others(start1=start1)
1080 s2 = _T00.others(start2=start2)
1081 try:
1082 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap,
1083 LatLon=LatLon, **LatLon_kwds)
1084 except (TypeError, ValueError) as x:
1085 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1088def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1089 height=None, wrap=True,
1090 LatLon=LatLon, **LatLon_kwds):
1091 '''Compute the intersection points of two circles each defined
1092 by a center point and a radius.
1094 @arg center1: Center of the first circle (L{LatLon}).
1095 @arg rad1: Radius of the first circle (C{meter} or C{radians},
1096 see B{C{radius}}).
1097 @arg center2: Center of the second circle (L{LatLon}).
1098 @arg rad2: Radius of the second circle (C{meter} or C{radians},
1099 see B{C{radius}}).
1100 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1101 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1102 @kwarg eps: Required overlap (C{meter} or C{radians}, see
1103 B{C{radius}}).
1104 @kwarg height: Optional height for the intersection points (C{meter},
1105 conventionally) or C{None} for the I{"radical height"}
1106 at the I{radical line} between both centers.
1107 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1108 @kwarg LatLon: Optional class to return the intersection
1109 points (L{LatLon}) or C{None}.
1110 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1111 arguments, ignored if C{B{LatLon} is None}.
1113 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1114 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1115 lon, height)}. For abutting circles, both intersection
1116 points are the same instance, aka I{radical center}.
1118 @raise IntersectionError: Concentric, antipodal, invalid or
1119 non-intersecting circles.
1121 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1123 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1124 B{C{eps}} or B{C{height}}.
1126 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1128 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1129 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1130 '''
1131 c1 = _T00.others(center1=center1)
1132 c2 = _T00.others(center2=center2)
1133 try:
1134 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1135 height=height, wrap=wrap,
1136 LatLon=LatLon, **LatLon_kwds)
1137 except (TypeError, ValueError) as x:
1138 raise _xError(x, center1=center1, rad1=rad1, center2=center2, rad2=rad2)
1141def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1142 height=None, too_d=None, wrap=True,
1143 LatLon=LatLon, **LatLon_kwds):
1144 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1145 # above, separated to allow callers to embellish any exceptions
1147 def _dest3(bearing, h):
1148 a, b = _destination2(a1, b1, r1, bearing)
1149 return _LL3Tuple(degrees90(a), degrees180(b), h,
1150 intersections2, LatLon, LatLon_kwds)
1152 r1, r2, f = _rads3(rad1, rad2, radius)
1153 if f: # swapped
1154 c1, c2 = c2, c1 # PYCHOK swap
1156 a1, b1 = c1.philam
1157 a2, b2 = c2.philam
1159 db, b2 = unrollPI(b1, b2, wrap=wrap)
1160 d = vincentys_(a2, a1, db) # radians
1161 if d < max(r1 - r2, EPS):
1162 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1164 r = eps if radius is None else (m2radians(
1165 eps, radius=radius) if eps else _0_0)
1166 if r < _0_0:
1167 raise _ValueError(eps=r)
1169 x = fsum_(r1, r2, -d) # overlap
1170 if x > max(r, EPS):
1171 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1172 x = sd * sr1
1173 if isnear0(x):
1174 raise _ValueError(_invalid_)
1175 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1177 elif x < r: # PYCHOK no cover
1178 t = (d * radius) if too_d is None else too_d
1179 raise IntersectionError(_too_(_Fmt.distant(t)))
1181 if height is None: # "radical height"
1182 f = _radical2(d, r1, r2).ratio
1183 h = Height(favg(c1.height, c2.height, f=f))
1184 else:
1185 h = Height(height)
1187 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1188 if x < EPS4: # externally ...
1189 r = _dest3(b, h)
1190 elif x > _PI_EPS4: # internally ...
1191 r = _dest3(b + PI, h)
1192 else:
1193 return _dest3(b + x, h), _dest3(b - x, h)
1194 return r, r # ... abutting circles
1197@deprecated_function
1198def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1199 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1200 '''
1201 return ispolar(points, wrap=wrap)
1204def _LL3Tuple(lat, lon, height, func, LatLon, LatLon_kwds):
1205 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1206 '''
1207 n = func.__name__
1208 if LatLon is None:
1209 r = LatLon3Tuple(lat, lon, height, name=n)
1210 else:
1211 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1212 r = LatLon(lat, lon, **kwds)
1213 return r
1216def meanOf(points, height=None, LatLon=LatLon, **LatLon_kwds):
1217 '''Compute the geographic mean of several points.
1219 @arg points: Points to be averaged (L{LatLon}[]).
1220 @kwarg height: Optional height at mean point, overriding the mean
1221 height (C{meter}).
1222 @kwarg LatLon: Optional class to return the mean point (L{LatLon})
1223 or C{None}.
1224 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1225 arguments, ignored if C{B{LatLon} is None}.
1227 @return: The geographic mean and height (B{C{LatLon}}) or a
1228 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}}
1229 is C{None}.
1231 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1233 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1234 '''
1235 # geographic mean
1236 m = sumOf(p._N_vector for p in _T00.PointsIter(points).iterate(closed=False))
1237 lat, lon = m._N_vector.latlon
1239 if height is None:
1240 h = fmean(p.height for p in _T00.PointsIter(points).iterate(closed=False))
1241 else: # PYCHOK no cover
1242 h = Height(height)
1243 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1246@deprecated_function
1247def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1248 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1250 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1251 point (L{LatLon}) on the polygon and the C{distance}
1252 between the C{closest} and the given B{C{point}}. The
1253 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1254 lon)} if B{C{LatLon}} is C{None} ...
1255 '''
1256 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1257 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1258 ll = LatLon2Tuple(ll.lat, ll.lon)
1259 return ll, d
1262def nearestOn3(point, points, closed=False, radius=R_M,
1263 LatLon=LatLon, **options):
1264 '''Locate the point on a path or polygon closest to a reference point.
1266 Distances are I{approximated} using function L{pygeodesy.equirectangular_},
1267 subject to the supplied B{C{options}}.
1269 @arg point: The reference point (L{LatLon}).
1270 @arg points: The path or polygon points (L{LatLon}[]).
1271 @kwarg closed: Optionally, close the polygon (C{bool}).
1272 @kwarg radius: Mean earth radius (C{meter}).
1273 @kwarg LatLon: Optional class to return the closest point (L{LatLon})
1274 or C{None}.
1275 @kwarg options: Optional keyword arguments for function
1276 L{pygeodesy.equirectangular_}.
1278 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1279 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1280 lon, height)} if B{C{LatLon}} is C{None}. The C{distance}
1281 is the L{pygeodesy.equirectangular_} distance between the
1282 C{closest} and the given B{C{point}} converted to C{meter},
1283 same units as B{C{radius}}. The C{angle} from the given
1284 B{C{point}} to the C{closest} is in compass C{degrees360},
1285 like function L{pygeodesy.compassAngle}. The C{height} is
1286 the (interpolated) height at the C{closest} point.
1288 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1289 see function L{pygeodesy.equirectangular_}.
1291 @raise PointsError: Insufficient number of B{C{points}}.
1293 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1295 @raise ValueError: Invalid B{C{radius}}.
1297 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}.
1298 '''
1299 lat, lon, d, c, h = _nearestOn5(point, points, closed=closed,
1300 LatLon=None, **options)
1301 d = degrees2m(d, radius=radius)
1302 n = nearestOn3.__name__
1303 r = LatLon3Tuple(lat, lon, h, name=n) if LatLon is None else \
1304 LatLon(lat, lon, height=h, name=n)
1305 return NearestOn3Tuple(r, d, c, name=n)
1308def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1309 '''Compute the perimeter of a (spherical) polygon or composite
1310 (with great circle arcs joining the points).
1312 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1313 or L{BooleanGH}).
1314 @kwarg closed: Optionally, close the polygon (C{bool}).
1315 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1316 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1318 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1319 or C{radians} if B{C{radius}} is C{None}).
1321 @raise PointsError: Insufficient number of B{C{points}}.
1323 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1325 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1326 C{B{points}} a composite.
1328 @note: Distances are based on function L{pygeodesy.vincentys_}.
1330 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf}
1331 and L{ellipsoidalKarney.perimeterOf}.
1332 '''
1333 def _rads(points, closed, wrap): # angular edge lengths in radians
1334 Ps = _T00.PointsIter(points, loop=1)
1335 a1, b1 = Ps[0].philam
1336 for i, p in Ps.enumerate(closed=closed):
1337 a2, b2 = p.philam
1338 db, b2 = unrollPI(b1, b2, wrap=wrap if i else False)
1339 yield vincentys_(a2, a1, db)
1340 a1, b1 = a2, b2
1342 if _MODS.booleans.isBoolean(points):
1343 if not closed:
1344 raise _ValueError(closed=closed, points=_composite_)
1345 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1346 else:
1347 r = fsum(_rads(points, closed, wrap), floats=True)
1348 return _r2m(r, radius)
1351def _r2m(r, radius):
1352 '''(INTERNAL) Angular distance in C{radians} to C{meter}.
1353 '''
1354 if radius is not None: # not in (None, _0_0)
1355 r *= R_M if radius is R_M else Radius(radius)
1356 return r
1359def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1360 excess=excessAbc,
1361 wrap=False):
1362 '''Compute the angles, sides, and area of a (spherical) triangle.
1364 @arg latA: First corner latitude (C{degrees}).
1365 @arg lonA: First corner longitude (C{degrees}).
1366 @arg latB: Second corner latitude (C{degrees}).
1367 @arg lonB: Second corner longitude (C{degrees}).
1368 @arg latC: Third corner latitude (C{degrees}).
1369 @arg lonC: Third corner longitude (C{degrees}).
1370 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1371 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1372 or C{None}.
1373 @kwarg excess: I{Spherical excess} callable (L{excessAbc},
1374 L{excessGirard} or L{excessLHuilier}).
1375 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
1377 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1378 spherical angles C{A}, C{B} and C{C}, angular sides
1379 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1380 in I{square} C{meter} or units of B{C{radius}}
1381 I{squared} or if B{C{radius}} is C{None}, a
1382 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in
1383 C{radians} with I{spherical excess} C{E} as the
1384 C{unit area} in C{radians}.
1385 '''
1386 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA),
1387 Phi_(latB=latB), Lam_(lonB=lonB),
1388 Phi_(latC=latC), Lam_(lonC=lonC),
1389 excess=excess, wrap=wrap)
1390 return _t7Tuple(t, radius)
1393def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc,
1394 wrap=False):
1395 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1396 excess} of a (spherical) triangle.
1398 @arg phiA: First corner latitude (C{radians}).
1399 @arg lamA: First corner longitude (C{radians}).
1400 @arg phiB: Second corner latitude (C{radians}).
1401 @arg lamB: Second corner longitude (C{radians}).
1402 @arg phiC: Third corner latitude (C{radians}).
1403 @arg lamC: Third corner longitude (C{radians}).
1404 @kwarg excess: I{Spherical excess} callable (L{excessAbc},
1405 L{excessGirard} or L{excessLHuilier}).
1406 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
1408 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1409 spherical angles C{A}, C{B} and C{C}, angular sides
1410 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1411 I{spherical excess} C{E}, all in C{radians}.
1412 '''
1413 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1414 d, _ = unrollPI(lamB, lamC, wrap=w)
1415 a = vincentys_(phiC, phiB, d)
1416 return a, (phiB, lamB, phiC, lamC, phiA, lamA)
1418 def _A_r(a, sa, ca, sb, cb, sc, cc):
1419 s = sb * sc
1420 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1421 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2's
1423 # notation: side C{a} is oposite to corner C{A}, etc.
1424 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1425 b, r = _a_r(wrap, *r)
1426 c, _ = _a_r(wrap, *r)
1428 A, r = _A_r(a, *sincos2_(a, b, c))
1429 B, r = _A_r(b, *r)
1430 C, _ = _A_r(c, *r)
1432 D = fsum_(PI2, -a, -b, -c, floats=True) # deficit aka defect
1433 E = excessGirard(A, B, C) if excess in (excessGirard, True) else (
1434 excessLHuilier(a, b, c) if excess in (excessLHuilier, False) else
1435 excessAbc(*max((A, b, c), (B, c, a), (C, a, b))))
1437 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1440def _t7Tuple(t, radius):
1441 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1442 '''
1443 if radius: # not in (None, _0_0)
1444 r = radius if isscalar(radius) else \
1445 _ellipsoidal_datum(radius).ellipsoid.Rmean
1446 A, B, C = map1(degrees, t.A, t.B, t.C)
1447 t = Triangle7Tuple(A, (r * t.a),
1448 B, (r * t.b),
1449 C, (r * t.c), t.E * r**2)
1450 return t
1453__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1454 areaOf, # functions
1455 intersecant2, intersection, intersections2, ispolar,
1456 isPoleEnclosedBy, # DEPRECATED, use ispolar
1457 meanOf,
1458 nearestOn2, nearestOn3,
1459 perimeterOf,
1460 sumOf, # == vector3d.sumOf
1461 triangle7, triangle8_)
1463# **) MIT License
1464#
1465# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1466#
1467# Permission is hereby granted, free of charge, to any person obtaining a
1468# copy of this software and associated documentation files (the "Software"),
1469# to deal in the Software without restriction, including without limitation
1470# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1471# and/or sell copies of the Software, and to permit persons to whom the
1472# Software is furnished to do so, subject to the following conditions:
1473#
1474# The above copyright notice and this permission notice shall be included
1475# in all copies or substantial portions of the Software.
1476#
1477# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1478# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1479# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1480# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1481# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1482# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1483# OTHER DEALINGS IN THE SOFTWARE.