Coverage for pygeodesy/sphericalTrigonometry.py: 95%
376 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-05 15:46 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-05 15:46 -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.04.05'
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 or composite encloses this point.
475 @arg points: The polygon points or composite (L{LatLon}[],
476 L{BooleanFHP} or L{BooleanGH}).
478 @return: C{True} if this point is inside the polygon or composite,
479 C{False} otherwise.
481 @raise PointsError: Insufficient number of B{C{points}}.
483 @raise TypeError: Some B{C{points}} are not L{LatLon}.
485 @raise ValueError: Invalid B{C{points}}, non-convex polygon.
487 @see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy}
488 and L{pygeodesy.ispolar} especially if the B{C{points}} may
489 enclose a pole or wrap around the earth I{longitudinally}.
490 '''
491 if _MODS.booleans.isBoolean(points):
492 return points._encloses(self.lat, self.lon)
494 Ps = self.PointsIter(points, loop=2, dedup=True)
495 n0 = self._N_vector
497 v2 = Ps[0]._N_vector
498 v1 = Ps[1]._N_vector
499 gc1 = v2.cross(v1)
500 # check whether this point on same side of all
501 # polygon edges (to the left or right depending
502 # on anti-/clockwise polygon direction)
503 t0 = gc1.angleTo(n0) > PI_2
504 s0 = None
505 # get great-circle vector for each edge
506 for i, p in Ps.enumerate(closed=True):
507 v2 = p._N_vector
508 gc = v1.cross(v2)
509 v1 = v2
511 t = gc.angleTo(n0) > PI_2
512 if t != t0: # different sides of edge i
513 return False # outside
515 # check for convex polygon: angle between
516 # gc vectors, signed by direction of n0
517 # (otherwise the test above is not reliable)
518 s = signOf(gc1.angleTo(gc, vSign=n0))
519 if s != s0:
520 if s0 is None:
521 s0 = s
522 else:
523 t = _Fmt.SQUARE(points=i)
524 raise _ValueError(t, p, txt=_not_(_convex_))
525 gc1 = gc
527 return True # inside
529 @deprecated_method
530 def isEnclosedBy(self, points): # PYCHOK no cover
531 '''DEPRECATED, use method C{isenclosedBy}.'''
532 return self.isenclosedBy(points)
534 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
535 '''Find the midpoint between this and an other point.
537 @arg other: The other point (L{LatLon}).
538 @kwarg height: Optional height for midpoint, overriding
539 the mean height (C{meter}).
540 @kwarg fraction: Midpoint location from this point (C{scalar}),
541 may be negative or greater than 1.0.
542 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
544 @return: Midpoint (L{LatLon}).
546 @raise TypeError: The B{C{other}} point is not L{LatLon}.
548 @raise ValueError: Invalid B{C{height}}.
550 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
552 @example:
554 >>> p1 = LatLon(52.205, 0.119)
555 >>> p2 = LatLon(48.857, 2.351)
556 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
557 '''
558 if fraction is _0_5:
559 self.others(other)
561 # see <https://MathForum.org/library/drmath/view/51822.html>
562 a1, b1 = self.philam
563 a2, b2 = other.philam
565 db, b2 = unrollPI(b1, b2, wrap=wrap)
567 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db)
569 x = ca2 * cdb + ca1
570 y = ca2 * sdb
572 a = atan2(sa1 + sa2, hypot(x, y))
573 b = atan2(y, x) + b1
575 h = self._havg(other) if height is None else Height(height)
576 r = self.classof(degrees90(a), degrees180(b), height=h)
577 else:
578 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
579 return r
581 def nearestOn(self, point1, point2, radius=R_M, **options):
582 '''Locate the point between two points closest to this point.
584 Distances are approximated by function L{pygeodesy.equirectangular_},
585 subject to the supplied B{C{options}}.
587 @arg point1: Start point (L{LatLon}).
588 @arg point2: End point (L{LatLon}).
589 @kwarg radius: Mean earth radius (C{meter}).
590 @kwarg options: Optional keyword arguments for function
591 L{pygeodesy.equirectangular_}.
593 @return: Closest point on the great circle path (L{LatLon}).
595 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
596 see function L{pygeodesy.equirectangular_}.
598 @raise NotImplementedError: Keyword argument C{B{within}=False}
599 is not (yet) supported.
601 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
603 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
605 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}
606 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
607 '''
608 # remove kwarg B{C{within}} if present
609 within = _xkwds_pop(options, within=True)
610 if not within:
611 notImplemented(self, within=within)
613# # UNTESTED - handle C{B{within}=False} and C{B{within}=True}
614# wrap = _xkwds_get(options, wrap=False)
615# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap)
616# if fabs(a) < EPS or (within and a < EPS):
617# return point1
618# d = point1.distanceTo(point2, radius=radius, wrap=wrap)
619# if isnear0(d):
620# return point1 # or point2
621# elif fabs(d - a) < EPS or (a + EPS) > d:
622# return point2
623# f = a / d
624# if within:
625# if f > EPS1:
626# return point2
627# elif f < EPS:
628# return point1
629# return point1.intermediateTo(point2, f, wrap=wrap)
631 # without kwarg B{C{within}}, use backward compatible .nearestOn3
632 return self.nearestOn3([point1, point2], closed=False, radius=radius,
633 **options)[0]
635 @deprecated_method
636 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover
637 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
639 @return: ... 2-Tuple C{(closest, distance)} of the closest
640 point (L{LatLon}) on the polygon and the distance
641 to that point from this point in C{meter}, same
642 units of B{C{radius}}.
643 '''
644 r = self.nearestOn3(points, closed=closed, radius=radius, **options)
645 return r.closest, r.distance
647 def nearestOn3(self, points, closed=False, radius=R_M, **options):
648 '''Locate the point on a polygon closest to this point.
650 Distances are approximated by function L{pygeodesy.equirectangular_},
651 subject to the supplied B{C{options}}.
653 @arg points: The polygon points (L{LatLon}[]).
654 @kwarg closed: Optionally, close the polygon (C{bool}).
655 @kwarg radius: Mean earth radius (C{meter}).
656 @kwarg options: Optional keyword arguments for function
657 L{pygeodesy.equirectangular_}.
659 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
660 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_}
661 C{distance} between this and the C{closest} point converted to
662 C{meter}, same units as B{C{radius}}. The C{angle} from this
663 to the C{closest} point is in compass C{degrees360}, like
664 function L{pygeodesy.compassAngle}.
666 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
667 see function L{pygeodesy.equirectangular_}.
669 @raise PointsError: Insufficient number of B{C{points}}.
671 @raise TypeError: Some B{C{points}} are not C{LatLon}.
673 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
675 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_}
676 and L{pygeodesy.nearestOn5}.
677 '''
678 if _LatLon_ in options:
679 raise _ValueError(**options)
681 lat, lon, d, c, h = _nearestOn5(self, points, closed=closed, **options)
682 return NearestOn3Tuple(self.classof(lat, lon, height=h),
683 degrees2m(d, radius=radius), c)
685 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
686 '''Convert this point to C{Karney}-based cartesian (ECEF)
687 coordinates.
689 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}}
690 and other keyword arguments, ignored
691 if C{B{Cartesian} is None}. Use
692 C{B{Cartesian}=...} to override
693 this L{Cartesian} class or specify
694 C{B{Cartesian}=None}.
696 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
697 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
698 C, M, datum)} with C{C} and C{M} if available.
700 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
701 '''
702 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
703 return LatLonSphericalBase.toCartesian(self, **kwds)
705 def _trackDistanceTo3(self, start, end, radius, wrap):
706 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo.
707 '''
708 self.others(start=start)
709 self.others(end=end)
711 r = Radius_(radius)
712 r = start.distanceTo(self, r, wrap=wrap) / r
714 b = radians(start.initialBearingTo(self, wrap=wrap))
715 e = radians(start.initialBearingTo(end, wrap=wrap))
716 x = asin(sin(r) * sin(b - e))
717 return r, x, (e - b)
719 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
720 '''Compute the angles, sides and area of a spherical triangle.
722 @arg otherB: Second triangle point (C{LatLon}).
723 @arg otherC: Third triangle point (C{LatLon}).
724 @kwarg radius: Mean earth radius, ellipsoid or datum
725 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},
726 L{Datum} or L{a_f2Tuple}) or C{None}.
727 @kwarg wrap: Wrap/unroll angular distances (C{bool}).
729 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if
730 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A,
731 a, B, b, C, c, D, E)}.
733 @see: Function L{triangle7} and U{Spherical trigonometry
734 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
735 '''
736 self.others(otherB=otherB)
737 self.others(otherC=otherC)
739 r = self.philam + otherB.philam + otherC.philam
740 t = triangle8_(*r, wrap=wrap)
741 return self._xnamed(_t7Tuple(t, radius))
743 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
744 area=True, eps=EPS1, radius=R_M, wrap=False):
745 '''Trilaterate three points by area overlap or perimeter intersection
746 of three corresponding circles.
748 @arg distance1: Distance to this point (C{meter}, same units
749 as B{C{radius}}).
750 @arg point2: Second center point (C{LatLon}).
751 @arg distance2: Distance to point2 (C{meter}, same units as
752 B{C{radius}}).
753 @arg point3: Third center point (C{LatLon}).
754 @arg distance3: Distance to point3 (C{meter}, same units as
755 B{C{radius}}).
756 @kwarg area: If C{True} compute the area overlap, otherwise the
757 perimeter intersection of the circles (C{bool}).
758 @kwarg eps: The required I{minimal overlap} for C{B{area}=True}
759 or the I{intersection margin} for C{B{area}=False}
760 (C{meter}, same units as B{C{radius}}).
761 @kwarg radius: Mean earth radius (C{meter}, conventionally).
762 @kwarg wrap: Wrap/unroll angular distances (C{bool}).
764 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
765 with C{min} and C{max} in C{meter}, same units as B{C{eps}},
766 the corresponding trilaterated points C{minPoint} and
767 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number
768 of trilatered points found for the given B{C{eps}}.
770 If only a single trilaterated point is found, C{min I{is}
771 max}, C{minPoint I{is} maxPoint} and C{n = 1}.
773 For C{B{area}=True}, C{min} and C{max} are the smallest
774 respectively largest I{radial} overlap found.
776 For C{B{area}=False}, C{min} and C{max} represent the
777 nearest respectively farthest intersection margin.
779 If C{B{area}=True} and all 3 circles are concentric, C{n =
780 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}}
781 with the smallest B{C{distance#}} C{min} and C{max} the
782 largest B{C{distance#}}.
784 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
785 insufficient overlap for C{B{area}=True} or
786 no intersection or all (near-)concentric for
787 C{B{area}=False}.
789 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
791 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}},
792 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.
793 '''
794 return _trilaterate5(self, distance1,
795 self.others(point2=point2), distance2,
796 self.others(point3=point3), distance3,
797 area=area, radius=radius, eps=eps, wrap=wrap)
800_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
803def areaOf(points, radius=R_M, wrap=True):
804 '''Calculate the area of a (spherical) polygon or composite
805 (with the pointsjoined by great circle arcs).
807 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
808 or L{BooleanGH}).
809 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
810 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
811 or C{None}.
812 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
814 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
815 or C{radians} if B{C{radius}} is C{None}).
817 @raise PointsError: Insufficient number of B{C{points}}.
819 @raise TypeError: Some B{C{points}} are not L{LatLon}.
821 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
823 @note: The area is based on I{Karney}'s U{'Area of a spherical
824 polygon'<https://MathOverflow.net/questions/97711/
825 the-area-of-spherical-polygons>}, 3rd Answer.
827 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
828 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
829 L{ellipsoidalKarney.areaOf}.
831 @example:
833 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
834 >>> areaOf(b) # 8666058750.718977
836 >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1)
837 >>> areaOf(c) # 6.18e9
838 '''
839 if _MODS.booleans.isBoolean(points):
840 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
842 Ps = _T00.PointsIter(points, loop=1)
843 p1 = p2 = Ps[0]
844 a1, b1 = p1.philam
845 ta1, z1 = tan_2(a1), None
847 A = Fsum() # mean phi
848 E = Fsum() # see L{pygeodesy.excessKarney_}
849 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
850 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
851 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
852 D = Fsum()
853 for i, p2 in Ps.enumerate(closed=True):
854 a2, b2 = p2.philam
855 db, b2 = unrollPI(b1, b2, wrap=wrap if i else False)
856 ta2 = tan_2(a2)
857 A += a2
858 E += atan2(tan_2(db, points=i) * (ta1 + ta2),
859 _1_0 + ta1 * ta2)
860 ta1, b1 = ta2, b2
862 if not p2.isequalTo(p1, eps=EPS):
863 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
864 if z1 is not None:
865 D += wrap180(z - z1) # (z - z1 + 540) ...
866 D += wrap180(z2 - z) # (z2 - z + 540) % 360 - 180
867 p1, z1 = p2, z2
869 R = fabs(E * _2_0)
870 if fabs(D) < _90_0: # ispolar(points)
871 R = fabs(R - PI2)
872 if radius:
873 a = degrees(A.fover(len(A))) # mean lat
874 R *= _mean_radius(radius, a)**2
875 return float(R)
878def _destination2(a, b, r, t):
879 '''(INTERNAL) Destination lat- and longitude in C{radians}.
881 @arg a: Latitude (C{radians}).
882 @arg b: Longitude (C{radians}).
883 @arg r: Angular distance (C{radians}).
884 @arg t: Bearing (compass C{radians}).
886 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
887 '''
888 # see <https://www.EdWilliams.org/avform.htm#LL>
889 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
890 ca *= sr
892 a = asin1(ct * ca + cr * sa)
893 d = atan2(st * ca, cr - sa * sin(a))
894 # note, in EdWilliams.org/avform.htm W is + and E is -
895 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
898def _int3d2(start, end, wrap, _i_, Vector, hs):
899 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
900 # and similar logic in .ellipsoidalBaseDI._intersect3
901 a1, b1 = start.philam
903 if isscalar(end): # bearing, get pseudo-end point
904 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
905 else: # must be a point
906 start.others(end, name=_end_ + _i_)
907 hs.append(end.height)
908 a2, b2 = end.philam
910 db, b2 = unrollPI(b1, b2, wrap=wrap)
911 if max(fabs(db), fabs(a2 - a1)) < EPS:
912 raise _ValueError(_SPACE_(_path_ + _i_, _null_))
913 # note, in EdWilliams.org/avform.htm W is + and E is -
914 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
915 -(b1 + b2) * _0_5)
916 cb21 *= sin(a1 - a2) # sa21
917 sb21 *= sin(a1 + a2) # sa12
918 x = Vector(sb12 * cb21 - cb12 * sb21,
919 cb12 * cb21 + sb12 * sb21,
920 cos(a1) * cos(a2) * sin(db)) # ll=start
921 return x.unit(), (db, (a2 - a1)) # negated d
924def _intdot(ds, a1, b1, a, b, wrap):
925 # compute dot product ds . (-b + b1, a - a1)
926 db, _ = unrollPI(b1, b, wrap=wrap)
927 return fdot(ds, db, a - a1)
930def intersecant2(center, circle, point, bearing, radius=R_M, exact=False,
931 height=None, wrap=True):
932 '''Compute the intersections of a circle and a line.
934 @arg center: Center of the circle (L{LatLon}).
935 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
936 or a point on the circle (L{LatLon}).
937 @arg point: A point in- or outside the circle (L{LatLon}).
938 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or
939 a second point on the line (L{LatLon}).
940 @kwarg radius: Mean earth radius (C{meter}, conventionally).
941 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
942 destination and distance, if C{False} use the basic
943 rhumb methods (C{bool}) or if C{None} use the I{great
944 circle} methods.
945 @kwarg height: Optional height for the intersection points (C{meter},
946 conventionally) or C{None}.
947 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
949 @return: 2-Tuple of the intersection points (representing a chord),
950 each an instance of this class. For a tangent line, each
951 point C{is} this very instance.
953 @raise IntersectionError: The circle and line do not intersect.
955 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
956 B{C{circle}} or B{C{bearing}} invalid.
958 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
959 B{C{exact}} or B{C{height}}.
960 '''
961 c = _T00.others(center=center)
962 p = _T00.others(point=point)
963 try:
964 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact,
965 height=height, wrap=wrap)
966 except (TypeError, ValueError) as x:
967 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing, exact=exact)
970def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
971 LatLon=None, **LatLon_kwds):
972 # (INTERNAL) Intersect two (spherical) path, see L{intersection}
973 # above, separated to allow callers to embellish any exceptions
975 hs = [start1.height, start2.height]
977 a1, b1 = start1.philam
978 a2, b2 = start2.philam
980 db, b2 = unrollPI(b1, b2, wrap=wrap)
981 r12 = vincentys_(a2, a1, db)
982 if fabs(r12) < EPS: # [nearly] coincident points
983 a, b = favg(a1, a2), favg(b1, b2)
985 # see <https://www.EdWilliams.org/avform.htm#Intersection>
986 elif isscalar(end1) and isscalar(end2): # both bearings
987 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
989 x1, x2 = (sr12 * ca1), (sr12 * ca2)
990 if isnear0(x1) or isnear0(x2):
991 raise IntersectionError(_parallel_)
992 # handle domain error for equivalent longitudes,
993 # see also functions asin_safe and acos_safe at
994 # <https://www.EdWilliams.org/avform.htm#Math>
995 t1, t2 = acos1((sa2 - sa1 * cr12) / x1), \
996 acos1((sa1 - sa2 * cr12) / x2)
997 if sin(db) > 0:
998 t12, t21 = t1, PI2 - t2
999 else:
1000 t12, t21 = PI2 - t1, t2
1001 t13, t23 = radiansPI2(end1), radiansPI2(end2)
1002 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
1003 wrapPI(t21 - t23)) # angle 1-2-3)
1004 if isnear0(sx1) and isnear0(sx2):
1005 raise IntersectionError(_infinite_)
1006 sx3 = sx1 * sx2
1007# XXX if sx3 < 0:
1008# XXX raise ValueError(_ambiguous_)
1009 x3 = acos1(cr12 * sx3 - cx2 * cx1)
1010 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
1012 a, b = _destination2(a1, b1, r13, t13)
1013 # like .ellipsoidalBaseDI,_intersect3, if this intersection
1014 # is "before" the first point, use the antipodal intersection
1015 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
1016 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1018 else: # end point(s) or bearing(s)
1019 _N_vector_ = _MODS.nvectorBase._N_vector_
1021 x1, d1 = _int3d2(start1, end1, wrap, _1_, _N_vector_, hs)
1022 x2, d2 = _int3d2(start2, end2, wrap, _2_, _N_vector_, hs)
1023 x = x1.cross(x2)
1024 if x.length < EPS: # [nearly] colinear or parallel paths
1025 raise IntersectionError(_colinear_)
1026 a, b = x.philam
1027 # choose intersection similar to sphericalNvector
1028 if not (_intdot(d1, a1, b1, a, b, wrap) *
1029 _intdot(d2, a2, b2, a, b, wrap)) > 0:
1030 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1032 h = fmean(hs) if height is None else Height(height)
1033 return _LL3Tuple(degrees90(a), degrees180(b), h,
1034 intersection, LatLon, LatLon_kwds)
1037def intersection(start1, end1, start2, end2, height=None, wrap=False,
1038 LatLon=LatLon, **LatLon_kwds):
1039 '''Compute the intersection point of two paths, each defined
1040 by two points or a start point and bearing from North.
1042 @arg start1: Start point of the first path (L{LatLon}).
1043 @arg end1: End point of the first path (L{LatLon}) or
1044 the initial bearing at the first start point
1045 (compass C{degrees360}).
1046 @arg start2: Start point of the second path (L{LatLon}).
1047 @arg end2: End point of the second path (L{LatLon}) or
1048 the initial bearing at the second start point
1049 (compass C{degrees360}).
1050 @kwarg height: Optional height for the intersection point,
1051 overriding the mean height (C{meter}).
1052 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1053 @kwarg LatLon: Optional class to return the intersection
1054 point (L{LatLon}) or C{None}.
1055 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1056 arguments, ignored if C{B{LatLon} is None}.
1058 @return: The intersection point as a (B{C{LatLon}}) or if
1059 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon,
1060 height)}. An alternate intersection point might
1061 be the L{antipode} to the returned result.
1063 @raise IntersectionError: Ambiguous or infinite intersection
1064 or colinear, parallel or otherwise
1065 non-intersecting paths.
1067 @raise TypeError: A B{C{start}} or B{C{end}} point not L{LatLon}.
1069 @raise ValueError: Invalid B{C{height}} or C{null} path.
1071 @example:
1073 >>> p = LatLon(51.8853, 0.2545)
1074 >>> s = LatLon(49.0034, 2.5735)
1075 >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E'
1076 '''
1077 s1 = _T00.others(start1=start1)
1078 s2 = _T00.others(start2=start2)
1079 try:
1080 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap,
1081 LatLon=LatLon, **LatLon_kwds)
1082 except (TypeError, ValueError) as x:
1083 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1086def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1087 height=None, wrap=True,
1088 LatLon=LatLon, **LatLon_kwds):
1089 '''Compute the intersection points of two circles each defined
1090 by a center point and a radius.
1092 @arg center1: Center of the first circle (L{LatLon}).
1093 @arg rad1: Radius of the first circle (C{meter} or C{radians},
1094 see B{C{radius}}).
1095 @arg center2: Center of the second circle (L{LatLon}).
1096 @arg rad2: Radius of the second circle (C{meter} or C{radians},
1097 see B{C{radius}}).
1098 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1099 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1100 @kwarg eps: Required overlap (C{meter} or C{radians}, see
1101 B{C{radius}}).
1102 @kwarg height: Optional height for the intersection points (C{meter},
1103 conventionally) or C{None} for the I{"radical height"}
1104 at the I{radical line} between both centers.
1105 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1106 @kwarg LatLon: Optional class to return the intersection
1107 points (L{LatLon}) or C{None}.
1108 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1109 arguments, ignored if C{B{LatLon} is None}.
1111 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1112 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1113 lon, height)}. For abutting circles, both intersection
1114 points are the same instance, aka I{radical center}.
1116 @raise IntersectionError: Concentric, antipodal, invalid or
1117 non-intersecting circles.
1119 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1121 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1122 B{C{eps}} or B{C{height}}.
1124 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1126 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1127 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1128 '''
1129 c1 = _T00.others(center1=center1)
1130 c2 = _T00.others(center2=center2)
1131 try:
1132 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1133 height=height, wrap=wrap,
1134 LatLon=LatLon, **LatLon_kwds)
1135 except (TypeError, ValueError) as x:
1136 raise _xError(x, center1=center1, rad1=rad1, center2=center2, rad2=rad2)
1139def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1140 height=None, too_d=None, wrap=True,
1141 LatLon=LatLon, **LatLon_kwds):
1142 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1143 # above, separated to allow callers to embellish any exceptions
1145 def _dest3(bearing, h):
1146 a, b = _destination2(a1, b1, r1, bearing)
1147 return _LL3Tuple(degrees90(a), degrees180(b), h,
1148 intersections2, LatLon, LatLon_kwds)
1150 r1, r2, f = _rads3(rad1, rad2, radius)
1151 if f: # swapped
1152 c1, c2 = c2, c1 # PYCHOK swap
1154 a1, b1 = c1.philam
1155 a2, b2 = c2.philam
1157 db, b2 = unrollPI(b1, b2, wrap=wrap)
1158 d = vincentys_(a2, a1, db) # radians
1159 if d < max(r1 - r2, EPS):
1160 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1162 r = eps if radius is None else (m2radians(
1163 eps, radius=radius) if eps else _0_0)
1164 if r < _0_0:
1165 raise _ValueError(eps=r)
1167 x = fsum_(r1, r2, -d) # overlap
1168 if x > max(r, EPS):
1169 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1170 x = sd * sr1
1171 if isnear0(x):
1172 raise _ValueError(_invalid_)
1173 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1175 elif x < r: # PYCHOK no cover
1176 t = (d * radius) if too_d is None else too_d
1177 raise IntersectionError(_too_(_Fmt.distant(t)))
1179 if height is None: # "radical height"
1180 f = _radical2(d, r1, r2).ratio
1181 h = Height(favg(c1.height, c2.height, f=f))
1182 else:
1183 h = Height(height)
1185 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1186 if x < EPS4: # externally ...
1187 r = _dest3(b, h)
1188 elif x > _PI_EPS4: # internally ...
1189 r = _dest3(b + PI, h)
1190 else:
1191 return _dest3(b + x, h), _dest3(b - x, h)
1192 return r, r # ... abutting circles
1195@deprecated_function
1196def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1197 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1198 '''
1199 return ispolar(points, wrap=wrap)
1202def _LL3Tuple(lat, lon, height, func, LatLon, LatLon_kwds):
1203 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1204 '''
1205 n = func.__name__
1206 if LatLon is None:
1207 r = LatLon3Tuple(lat, lon, height, name=n)
1208 else:
1209 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1210 r = LatLon(lat, lon, **kwds)
1211 return r
1214def meanOf(points, height=None, LatLon=LatLon, **LatLon_kwds):
1215 '''Compute the geographic mean of several points.
1217 @arg points: Points to be averaged (L{LatLon}[]).
1218 @kwarg height: Optional height at mean point, overriding the mean
1219 height (C{meter}).
1220 @kwarg LatLon: Optional class to return the mean point (L{LatLon})
1221 or C{None}.
1222 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1223 arguments, ignored if C{B{LatLon} is None}.
1225 @return: The geographic mean and height (B{C{LatLon}}) or a
1226 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}}
1227 is C{None}.
1229 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1231 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1232 '''
1233 # geographic mean
1234 m = sumOf(p._N_vector for p in _T00.PointsIter(points).iterate(closed=False))
1235 lat, lon = m._N_vector.latlon
1237 if height is None:
1238 h = fmean(p.height for p in _T00.PointsIter(points).iterate(closed=False))
1239 else: # PYCHOK no cover
1240 h = Height(height)
1241 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1244@deprecated_function
1245def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1246 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1248 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1249 point (L{LatLon}) on the polygon and the C{distance}
1250 between the C{closest} and the given B{C{point}}. The
1251 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1252 lon)} if B{C{LatLon}} is C{None} ...
1253 '''
1254 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1255 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1256 ll = LatLon2Tuple(ll.lat, ll.lon)
1257 return ll, d
1260def nearestOn3(point, points, closed=False, radius=R_M,
1261 LatLon=LatLon, **options):
1262 '''Locate the point on a path or polygon closest to a reference point.
1264 Distances are I{approximated} using function L{pygeodesy.equirectangular_},
1265 subject to the supplied B{C{options}}.
1267 @arg point: The reference point (L{LatLon}).
1268 @arg points: The path or polygon points (L{LatLon}[]).
1269 @kwarg closed: Optionally, close the polygon (C{bool}).
1270 @kwarg radius: Mean earth radius (C{meter}).
1271 @kwarg LatLon: Optional class to return the closest point (L{LatLon})
1272 or C{None}.
1273 @kwarg options: Optional keyword arguments for function
1274 L{pygeodesy.equirectangular_}.
1276 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1277 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1278 lon, height)} if B{C{LatLon}} is C{None}. The C{distance}
1279 is the L{pygeodesy.equirectangular_} distance between the
1280 C{closest} and the given B{C{point}} converted to C{meter},
1281 same units as B{C{radius}}. The C{angle} from the given
1282 B{C{point}} to the C{closest} is in compass C{degrees360},
1283 like function L{pygeodesy.compassAngle}. The C{height} is
1284 the (interpolated) height at the C{closest} point.
1286 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1287 see function L{pygeodesy.equirectangular_}.
1289 @raise PointsError: Insufficient number of B{C{points}}.
1291 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1293 @raise ValueError: Invalid B{C{radius}}.
1295 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}.
1296 '''
1297 lat, lon, d, c, h = _nearestOn5(point, points, closed=closed,
1298 LatLon=None, **options)
1299 d = degrees2m(d, radius=radius)
1300 n = nearestOn3.__name__
1301 r = LatLon3Tuple(lat, lon, h, name=n) if LatLon is None else \
1302 LatLon(lat, lon, height=h, name=n)
1303 return NearestOn3Tuple(r, d, c, name=n)
1306def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1307 '''Compute the perimeter of a (spherical) polygon or composite
1308 (with great circle arcs joining the points).
1310 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1311 or L{BooleanGH}).
1312 @kwarg closed: Optionally, close the polygon (C{bool}).
1313 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1314 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1316 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1317 or C{radians} if B{C{radius}} is C{None}).
1319 @raise PointsError: Insufficient number of B{C{points}}.
1321 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1323 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1324 C{B{points}} a composite.
1326 @note: Distances are based on function L{pygeodesy.vincentys_}.
1328 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf}
1329 and L{ellipsoidalKarney.perimeterOf}.
1330 '''
1331 def _rads(points, closed, wrap): # angular edge lengths in radians
1332 Ps = _T00.PointsIter(points, loop=1)
1333 a1, b1 = Ps[0].philam
1334 for i, p in Ps.enumerate(closed=closed):
1335 a2, b2 = p.philam
1336 db, b2 = unrollPI(b1, b2, wrap=wrap if i else False)
1337 yield vincentys_(a2, a1, db)
1338 a1, b1 = a2, b2
1340 if _MODS.booleans.isBoolean(points):
1341 if not closed:
1342 raise _ValueError(closed=closed, points=_composite_)
1343 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1344 else:
1345 r = fsum(_rads(points, closed, wrap), floats=True)
1346 return _r2m(r, radius)
1349def _r2m(r, radius):
1350 '''(INTERNAL) Angular distance in C{radians} to C{meter}.
1351 '''
1352 if radius is not None: # not in (None, _0_0)
1353 r *= R_M if radius is R_M else Radius(radius)
1354 return r
1357def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1358 excess=excessAbc_,
1359 wrap=False):
1360 '''Compute the angles, sides, and area of a (spherical) triangle.
1362 @arg latA: First corner latitude (C{degrees}).
1363 @arg lonA: First corner longitude (C{degrees}).
1364 @arg latB: Second corner latitude (C{degrees}).
1365 @arg lonB: Second corner longitude (C{degrees}).
1366 @arg latC: Third corner latitude (C{degrees}).
1367 @arg lonC: Third corner longitude (C{degrees}).
1368 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1369 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1370 or C{None}.
1371 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1372 L{excessGirard_} or L{excessLHuilier_}).
1373 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
1375 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1376 spherical angles C{A}, C{B} and C{C}, angular sides
1377 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1378 in I{square} C{meter} or same units as B{C{radius}}
1379 I{squared} or if C{B{radius}=0} or C{None}, a
1380 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in
1381 C{radians} with the I{spherical excess} C{E} as the
1382 C{unit area} in C{radians}.
1383 '''
1384 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA),
1385 Phi_(latB=latB), Lam_(lonB=lonB),
1386 Phi_(latC=latC), Lam_(lonC=lonC),
1387 excess=excess, wrap=wrap)
1388 return _t7Tuple(t, radius)
1391def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_,
1392 wrap=False):
1393 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1394 excess} of a (spherical) triangle.
1396 @arg phiA: First corner latitude (C{radians}).
1397 @arg lamA: First corner longitude (C{radians}).
1398 @arg phiB: Second corner latitude (C{radians}).
1399 @arg lamB: Second corner longitude (C{radians}).
1400 @arg phiC: Third corner latitude (C{radians}).
1401 @arg lamC: Third corner longitude (C{radians}).
1402 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1403 L{excessGirard_} or L{excessLHuilier_}).
1404 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
1406 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1407 spherical angles C{A}, C{B} and C{C}, angular sides
1408 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1409 I{spherical excess} C{E}, all in C{radians}.
1410 '''
1411 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1412 d, _ = unrollPI(lamB, lamC, wrap=w)
1413 a = vincentys_(phiC, phiB, d)
1414 return a, (phiB, lamB, phiC, lamC, phiA, lamA)
1416 def _A_r(a, sa, ca, sb, cb, sc, cc):
1417 s = sb * sc
1418 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1419 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2's
1421 # notation: side C{a} is oposite to corner C{A}, etc.
1422 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1423 b, r = _a_r(wrap, *r)
1424 c, _ = _a_r(wrap, *r)
1426 A, r = _A_r(a, *sincos2_(a, b, c))
1427 B, r = _A_r(b, *r)
1428 C, _ = _A_r(c, *r)
1430 D = fsum_(PI2, -a, -b, -c, floats=True) # deficit aka defect
1431 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else (
1432 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else
1433 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b))))
1435 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1438def _t7Tuple(t, radius):
1439 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1440 '''
1441 if radius: # not in (None, _0_0)
1442 r = radius if isscalar(radius) else \
1443 _ellipsoidal_datum(radius).ellipsoid.Rmean
1444 A, B, C = map1(degrees, t.A, t.B, t.C)
1445 t = Triangle7Tuple(A, (r * t.a),
1446 B, (r * t.b),
1447 C, (r * t.c), t.E * r**2)
1448 return t
1451__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1452 areaOf, # functions
1453 intersecant2, intersection, intersections2, ispolar,
1454 isPoleEnclosedBy, # DEPRECATED, use ispolar
1455 meanOf,
1456 nearestOn2, nearestOn3,
1457 perimeterOf,
1458 sumOf, # == vector3d.sumOf
1459 triangle7, triangle8_)
1461# **) MIT License
1462#
1463# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1464#
1465# Permission is hereby granted, free of charge, to any person obtaining a
1466# copy of this software and associated documentation files (the "Software"),
1467# to deal in the Software without restriction, including without limitation
1468# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1469# and/or sell copies of the Software, and to permit persons to whom the
1470# Software is furnished to do so, subject to the following conditions:
1471#
1472# The above copyright notice and this permission notice shall be included
1473# in all copies or substantial portions of the Software.
1474#
1475# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1476# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1477# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1478# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1479# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1480# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1481# OTHER DEALINGS IN THE SOFTWARE.