Coverage for pygeodesy/sphericalTrigonometry.py: 95%
375 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-11 14:35 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-11 14:35 -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_, _invalid_, \
34 _LatLon_, _line_, _near_, _not_, _null_, _points_, \
35 _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.11'
59_parallel_ = 'parallel'
61_PI_EPS4 = PI - EPS4
62if _PI_EPS4 >= PI:
63 raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4)
66class Cartesian(CartesianSphericalBase):
67 '''Extended to convert geocentric, L{Cartesian} points to
68 spherical, geodetic L{LatLon}.
69 '''
71 def toLatLon(self, **LatLon_and_kwds): # PYCHOK LatLon=LatLon
72 '''Convert this cartesian point to a C{spherical} geodetic point.
74 @kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword
75 arguments. Use C{B{LatLon}=...} to override
76 this L{LatLon} class or specify C{B{LatLon}=None}.
78 @return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is C{None},
79 an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)}
80 with C{C} and C{M} if available.
82 @raise TypeError: Invalid B{C{LatLon_and_kwds}} argument.
83 '''
84 kwds = _xkwds(LatLon_and_kwds, LatLon=LatLon, datum=self.datum)
85 return CartesianSphericalBase.toLatLon(self, **kwds)
88class LatLon(LatLonSphericalBase):
89 '''New point on spherical model earth model.
91 @example:
93 >>> p = LatLon(52.205, 0.119) # height=0
94 '''
96 def alongTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
97 '''Compute the (angular) distance (signed) from the start to
98 the closest point on the great circle line defined by a
99 start and an end point.
101 That is, if a perpendicular is drawn from this point to the
102 great circle line, the along-track distance is the distance
103 from the start point to the point where the perpendicular
104 crosses the line.
106 @arg start: Start point of great circle line (L{LatLon}).
107 @arg end: End point of great circle line (L{LatLon}).
108 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
109 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
111 @return: Distance along the great circle line (C{meter},
112 same units as B{C{radius}}) or C{radians} if
113 C{B{radius} is None}, positive if after the B{C{start}}
114 toward the B{C{end}} point of the line or negative
115 if before the B{C{start}} point.
117 @raise TypeError: Invalid B{C{start}} or B{C{end}} point.
119 @raise ValueError: Invalid B{C{radius}}.
121 @example:
123 >>> p = LatLon(53.2611, -0.7972)
125 >>> s = LatLon(53.3206, -1.7297)
126 >>> e = LatLon(53.1887, 0.1334)
127 >>> d = p.alongTrackDistanceTo(s, e) # 62331.58
128 '''
129 r, x, b = self._trackDistanceTo3(start, end, radius, wrap)
130 cx = cos(x)
131 return _0_0 if isnear0(cx) else \
132 _r2m(copysign0(acos1(cos(r) / cx), cos(b)), radius)
134 @deprecated_method
135 def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover
136 '''DEPRECATED, use method L{initialBearingTo}.
137 '''
138 return self.initialBearingTo(other, wrap=wrap, raiser=raiser)
140 def crossingParallels(self, other, lat, wrap=False):
141 '''Return the pair of meridians at which a great circle defined
142 by this and an other point crosses the given latitude.
144 @arg other: The other point defining great circle (L{LatLon}).
145 @arg lat: Latitude at the crossing (C{degrees}).
146 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
148 @return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
149 C{None} if the great circle doesn't reach B{C{lat}}.
150 '''
151 self.others(other)
153 a1, b1 = self.philam
154 a2, b2 = other.philam
156 a = radians(lat)
157 db, b2 = unrollPI(b1, b2, wrap=wrap)
159 sa, ca, sa1, ca1, \
160 sa2, ca2, sdb, cdb = sincos2_(a, a1, a2, db)
161 sa1 *= ca2 * ca
163 x = sa1 * sdb
164 y = sa1 * cdb - ca1 * sa2 * ca
165 z = ca1 * sdb * ca2 * sa
167 h = hypot(x, y)
168 if h < EPS or fabs(z) > h: # PYCHOK no cover
169 return None # great circle doesn't reach latitude
171 m = atan2(-y, x) + b1 # longitude at max latitude
172 d = acos1(z / h) # delta longitude to intersections
173 return degrees180(m - d), degrees180(m + d)
175 def crossTrackDistanceTo(self, start, end, radius=R_M, wrap=False):
176 '''Compute the (angular) distance (signed) from this point to
177 the great circle defined by a start and an end point.
179 @arg start: Start point of great circle line (L{LatLon}).
180 @arg end: End point of great circle line (L{LatLon}).
181 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
182 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
184 @return: Distance to great circle (negative if to the left
185 or positive if to the right of the line) (C{meter},
186 same units as B{C{radius}} or C{radians} if
187 B{C{radius}} is C{None}).
189 @raise TypeError: The B{C{start}} or B{C{end}} point is not L{LatLon}.
191 @raise ValueError: Invalid B{C{radius}}.
193 @example:
195 >>> p = LatLon(53.2611, -0.7972)
197 >>> s = LatLon(53.3206, -1.7297)
198 >>> e = LatLon(53.1887, 0.1334)
199 >>> d = p.crossTrackDistanceTo(s, e) # -307.5
200 '''
201 _, x, _ = self._trackDistanceTo3(start, end, radius, wrap)
202 return _r2m(x, radius)
204 def destination(self, distance, bearing, radius=R_M, height=None):
205 '''Locate the destination from this point after having
206 travelled the given distance on the given initial bearing.
208 @arg distance: Distance travelled (C{meter}, same units as
209 B{C{radius}}).
210 @arg bearing: Bearing from this point (compass C{degrees360}).
211 @kwarg radius: Mean earth radius (C{meter}).
212 @kwarg height: Optional height at destination (C{meter}, same
213 units a B{C{radius}}).
215 @return: Destination point (L{LatLon}).
217 @raise ValueError: Invalid B{C{distance}}, B{C{bearing}},
218 B{C{radius}} or B{C{height}}.
220 @example:
222 >>> p1 = LatLon(51.4778, -0.0015)
223 >>> p2 = p1.destination(7794, 300.7)
224 >>> p2.toStr() # '51.5135°N, 000.0983°W'
225 '''
226 a, b = self.philam
227 r, t = _angular(distance, radius), Bearing_(bearing)
229 a, b = _destination2(a, b, r, t)
230 h = self.height if height is None else Height(height)
231 return self.classof(degrees90(a), degrees180(b), height=h)
233 def distanceTo(self, other, radius=R_M, wrap=False):
234 '''Compute the (angular) distance from this to an other point.
236 @arg other: The other point (L{LatLon}).
237 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
238 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
240 @return: Distance between this and the B{C{other}} point
241 (C{meter}, same units as B{C{radius}} or
242 C{radians} if B{C{radius}} is C{None}).
244 @raise TypeError: The B{C{other}} point is not L{LatLon}.
246 @raise ValueError: Invalid B{C{radius}}.
248 @example:
250 >>> p1 = LatLon(52.205, 0.119)
251 >>> p2 = LatLon(48.857, 2.351);
252 >>> d = p1.distanceTo(p2) # 404300
253 '''
254 self.others(other)
256 a1, b1 = self.philam
257 a2, b2 = other.philam
259 db, _ = unrollPI(b1, b2, wrap=wrap)
260 return _r2m(vincentys_(a2, a1, db), radius)
262# @Property_RO
263# def Ecef(self):
264# '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}.
265# '''
266# return _MODS.ecef.EcefKarney
268 def greatCircle(self, bearing):
269 '''Compute the vector normal to great circle obtained by heading
270 on the given initial bearing from this point.
272 Direction of vector is such that initial bearing vector
273 b = c × n, where n is an n-vector representing this point.
275 @arg bearing: Bearing from this point (compass C{degrees360}).
277 @return: Vector representing great circle (L{Vector3d}).
279 @raise ValueError: Invalid B{C{bearing}}.
281 @example:
283 >>> p = LatLon(53.3206, -1.7297)
284 >>> g = p.greatCircle(96.0)
285 >>> g.toStr() # (-0.794, 0.129, 0.594)
286 '''
287 a, b = self.philam
289 t = Bearing_(bearing)
290 sa, ca, sb, cb, st, ct = sincos2_(a, b, t)
292 return Vector3d(sb * ct - cb * sa * st,
293 -cb * ct - sb * sa * st,
294 ca * st) # XXX .unit()?
296 def initialBearingTo(self, other, wrap=False, raiser=False):
297 '''Compute the initial bearing (forward azimuth) from this
298 to an other point.
300 @arg other: The other point (spherical L{LatLon}).
301 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
302 @kwarg raiser: Optionally, raise L{CrossError} (C{bool}),
303 use C{B{raiser}=True} for behavior like
304 C{sphericalNvector.LatLon.initialBearingTo}.
306 @return: Initial bearing (compass C{degrees360}).
308 @raise CrossError: If this and the B{C{other}} point coincide,
309 provided both B{C{raiser}} is C{True} and
310 L{pygeodesy.crosserrors} is C{True}.
312 @raise TypeError: The B{C{other}} point is not L{LatLon}.
314 @example:
316 >>> p1 = LatLon(52.205, 0.119)
317 >>> p2 = LatLon(48.857, 2.351)
318 >>> b = p1.initialBearingTo(p2) # 156.2
319 '''
320 self.others(other)
322 a1, b1 = self.philam
323 a2, b2 = other.philam
325 # XXX behavior like sphericalNvector.LatLon.initialBearingTo
326 if raiser and crosserrors() and max(fabs(a2 - a1), fabs(b2 - b1)) < EPS:
327 raise CrossError(_points_, self, txt=_coincident_)
329 return degrees(bearing_(a1, b1, a2, b2, final=False, wrap=wrap))
331 def intermediateTo(self, other, fraction, height=None, wrap=False):
332 '''Locate the point at given fraction between (or along) this
333 and an other point.
335 @arg other: The other point (L{LatLon}).
336 @arg fraction: Fraction between both points (C{scalar},
337 0.0 at this and 1.0 at the other point).
338 @kwarg height: Optional height, overriding the intermediate
339 height (C{meter}).
340 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
342 @return: Intermediate point (L{LatLon}).
344 @raise TypeError: The B{C{other}} point is not L{LatLon}.
346 @raise ValueError: Invalid B{C{fraction}} or B{C{height}}.
348 @see: Methods C{midpointTo} and C{rhumbMidpointTo}.
350 @example:
352 >>> p1 = LatLon(52.205, 0.119)
353 >>> p2 = LatLon(48.857, 2.351)
354 >>> p = p1.intermediateTo(p2, 0.25) # 51.3721°N, 000.7073°E
355 '''
356 f = Scalar(fraction=fraction) # not high=_1_0
357 if isnear0(f): # PYCHOK no cover
358 r = self
359 elif isnear1(f) and not wrap: # PYCHOK no cover
360 r = self.others(other)
361 else:
362 a1, b1 = self.philam
363 a2, b2 = self.others(other).philam
365 db, b2 = unrollPI(b1, b2, wrap=wrap)
366 r = vincentys_(a2, a1, db)
367 sr = sin(r)
368 if isnon0(sr):
369 sa1, ca1, sa2, ca2, \
370 sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2)
372 t = f * r
373 a = sin(r - t) # / sr superflous
374 b = sin( t) # / sr superflous
376 x = a * ca1 * cb1 + b * ca2 * cb2
377 y = a * ca1 * sb1 + b * ca2 * sb2
378 z = a * sa1 + b * sa2
380 a = atan2(z, hypot(x, y))
381 b = atan2(y, x)
383 else: # PYCHOK no cover
384 a = favg(a1, a2, f=f) # coincident
385 b = favg(b1, b2, f=f)
387 h = self._havg(other, f=f) if height is None else Height(height)
388 r = self.classof(degrees90(a), degrees180(b), height=h)
389 return r
391 def intersection(self, end1, other, end2, height=None, wrap=False):
392 '''Compute the intersection point of two lines, each defined
393 by two points or a start point and bearing from North.
395 @arg end1: End point of this line (L{LatLon}) or the
396 initial bearing at this point (compass
397 C{degrees360}).
398 @arg other: Start point of the other line (L{LatLon}).
399 @arg end2: End point of the other line (L{LatLon}) or
400 the initial bearing at the other start point
401 (compass C{degrees360}).
402 @kwarg height: Optional height for intersection point,
403 overriding the mean height (C{meter}).
404 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
406 @return: The intersection point (L{LatLon}). An alternate
407 intersection point might be the L{antipode} to
408 the returned result.
410 @raise IntersectionError: Ambiguous or infinite intersection
411 or colinear, parallel or otherwise
412 non-intersecting lines.
414 @raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}}
415 or B{C{end2}} not C{scalar} nor L{LatLon}.
417 @raise ValueError: Invalid B{C{height}} or C{null} line.
419 @example:
421 >>> p = LatLon(51.8853, 0.2545)
422 >>> s = LatLon(49.0034, 2.5735)
423 >>> i = p.intersection(108.547, s, 32.435) # '50.9078°N, 004.5084°E'
424 '''
425 try:
426 s2 = self.others(other)
427 return _intersect(self, end1, s2, end2, height=height, wrap=wrap,
428 LatLon=self.classof)
429 except (TypeError, ValueError) as x:
430 raise _xError(x, start1=self, end1=end1, other=other, end2=end2)
432 def intersections2(self, rad1, other, rad2, radius=R_M, eps=_0_0,
433 height=None, wrap=True):
434 '''Compute the intersection points of two circles, each defined
435 by a center point and radius.
437 @arg rad1: Radius of the this circle (C{meter} or C{radians},
438 see B{C{radius}}).
439 @arg other: Center point of the other circle (L{LatLon}).
440 @arg rad2: Radius of the other circle (C{meter} or C{radians},
441 see B{C{radius}}).
442 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
443 B{C{rad2}} and B{C{eps}} are given in C{radians}).
444 @kwarg eps: Required overlap (C{meter} or C{radians}, see
445 B{C{radius}}).
446 @kwarg height: Optional height for the intersection points (C{meter},
447 conventionally) or C{None} for the I{"radical height"}
448 at the I{radical line} between both centers.
449 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
451 @return: 2-Tuple of the intersection points, each a L{LatLon}
452 instance. For abutting circles, both intersection
453 points are the same instance, aka the I{radical center}.
455 @raise IntersectionError: Concentric, antipodal, invalid or
456 non-intersecting circles.
458 @raise TypeError: If B{C{other}} is not L{LatLon}.
460 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
461 B{C{eps}} or B{C{height}}.
462 '''
463 try:
464 c2 = self.others(other)
465 return _intersects2(self, rad1, c2, rad2, radius=radius, eps=eps,
466 height=height, wrap=wrap,
467 LatLon=self.classof)
468 except (TypeError, ValueError) as x:
469 raise _xError(x, center=self, rad1=rad1, other=other, rad2=rad2)
471 def isenclosedBy(self, points):
472 '''Check whether a (convex) polygon or composite encloses this point.
474 @arg points: The polygon points or composite (L{LatLon}[],
475 L{BooleanFHP} or L{BooleanGH}).
477 @return: C{True} if this point is inside the polygon or composite,
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 I{longitudinally}.
489 '''
490 if _MODS.booleans.isBoolean(points):
491 return points._encloses(self.lat, self.lon)
493 Ps = self.PointsIter(points, loop=2, dedup=True)
494 n0 = self._N_vector
496 v2 = Ps[0]._N_vector
497 v1 = Ps[1]._N_vector
498 gc1 = v2.cross(v1)
499 # check whether this point on same side of all
500 # polygon edges (to the left or right depending
501 # on anti-/clockwise polygon direction)
502 t0 = gc1.angleTo(n0) > PI_2
503 s0 = None
504 # get great-circle vector for each edge
505 for i, p in Ps.enumerate(closed=True):
506 v2 = p._N_vector
507 gc = v1.cross(v2)
508 v1 = v2
510 t = gc.angleTo(n0) > PI_2
511 if t != t0: # different sides of edge i
512 return False # outside
514 # check for convex polygon: angle between
515 # gc vectors, signed by direction of n0
516 # (otherwise the test above is not reliable)
517 s = signOf(gc1.angleTo(gc, vSign=n0))
518 if s != s0:
519 if s0 is None:
520 s0 = s
521 else:
522 t = _Fmt.SQUARE(points=i)
523 raise _ValueError(t, p, txt=_not_(_convex_))
524 gc1 = gc
526 return True # inside
528 @deprecated_method
529 def isEnclosedBy(self, points): # PYCHOK no cover
530 '''DEPRECATED, use method C{isenclosedBy}.'''
531 return self.isenclosedBy(points)
533 def midpointTo(self, other, height=None, fraction=_0_5, wrap=False):
534 '''Find the midpoint between this and an other point.
536 @arg other: The other point (L{LatLon}).
537 @kwarg height: Optional height for midpoint, overriding
538 the mean height (C{meter}).
539 @kwarg fraction: Midpoint location from this point (C{scalar}),
540 may be negative or greater than 1.0.
541 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
543 @return: Midpoint (L{LatLon}).
545 @raise TypeError: The B{C{other}} point is not L{LatLon}.
547 @raise ValueError: Invalid B{C{height}}.
549 @see: Methods C{intermediateTo} and C{rhumbMidpointTo}.
551 @example:
553 >>> p1 = LatLon(52.205, 0.119)
554 >>> p2 = LatLon(48.857, 2.351)
555 >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E'
556 '''
557 if fraction is _0_5:
558 self.others(other)
560 # see <https://MathForum.org/library/drmath/view/51822.html>
561 a1, b1 = self.philam
562 a2, b2 = other.philam
564 db, b2 = unrollPI(b1, b2, wrap=wrap)
566 sa1, ca1, sa2, ca2, sdb, cdb = sincos2_(a1, a2, db)
568 x = ca2 * cdb + ca1
569 y = ca2 * sdb
571 a = atan2(sa1 + sa2, hypot(x, y))
572 b = atan2(y, x) + b1
574 h = self._havg(other) if height is None else Height(height)
575 r = self.classof(degrees90(a), degrees180(b), height=h)
576 else:
577 r = self.intermediateTo(other, fraction, height=height, wrap=wrap)
578 return r
580 def nearestOn(self, point1, point2, radius=R_M, **options):
581 '''Locate the point between two points closest to this point.
583 Distances are approximated by function L{pygeodesy.equirectangular_},
584 subject to the supplied B{C{options}}.
586 @arg point1: Start point (L{LatLon}).
587 @arg point2: End point (L{LatLon}).
588 @kwarg radius: Mean earth radius (C{meter}).
589 @kwarg options: Optional keyword arguments for function
590 L{pygeodesy.equirectangular_}.
592 @return: Closest point on the great circle line (L{LatLon}).
594 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
595 see function L{pygeodesy.equirectangular_}.
597 @raise NotImplementedError: Keyword argument C{B{within}=False}
598 is not (yet) supported.
600 @raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
602 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
604 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}
605 and method L{sphericalTrigonometry.LatLon.nearestOn3}.
606 '''
607 # remove kwarg B{C{within}} if present
608 within = _xkwds_pop(options, within=True)
609 if not within:
610 notImplemented(self, within=within)
612# # UNTESTED - handle C{B{within}=False} and C{B{within}=True}
613# wrap = _xkwds_get(options, wrap=False)
614# a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap)
615# if fabs(a) < EPS or (within and a < EPS):
616# return point1
617# d = point1.distanceTo(point2, radius=radius, wrap=wrap)
618# if isnear0(d):
619# return point1 # or point2
620# elif fabs(d - a) < EPS or (a + EPS) > d:
621# return point2
622# f = a / d
623# if within:
624# if f > EPS1:
625# return point2
626# elif f < EPS:
627# return point1
628# return point1.intermediateTo(point2, f, wrap=wrap)
630 # without kwarg B{C{within}}, use backward compatible .nearestOn3
631 return self.nearestOn3([point1, point2], closed=False, radius=radius,
632 **options)[0]
634 @deprecated_method
635 def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover
636 '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
638 @return: ... 2-Tuple C{(closest, distance)} of the closest
639 point (L{LatLon}) on the polygon and the distance
640 to that point from this point in C{meter}, same
641 units of B{C{radius}}.
642 '''
643 r = self.nearestOn3(points, closed=closed, radius=radius, **options)
644 return r.closest, r.distance
646 def nearestOn3(self, points, closed=False, radius=R_M, **options):
647 '''Locate the point on a polygon closest to this point.
649 Distances are approximated by function L{pygeodesy.equirectangular_},
650 subject to the supplied B{C{options}}.
652 @arg points: The polygon points (L{LatLon}[]).
653 @kwarg closed: Optionally, close the polygon (C{bool}).
654 @kwarg radius: Mean earth radius (C{meter}).
655 @kwarg options: Optional keyword arguments for function
656 L{pygeodesy.equirectangular_}.
658 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the
659 C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_}
660 C{distance} between this and the C{closest} point converted to
661 C{meter}, same units as B{C{radius}}. The C{angle} from this
662 to the C{closest} point is in compass C{degrees360}, like
663 function L{pygeodesy.compassAngle}.
665 @raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}},
666 see function L{pygeodesy.equirectangular_}.
668 @raise PointsError: Insufficient number of B{C{points}}.
670 @raise TypeError: Some B{C{points}} are not C{LatLon}.
672 @raise ValueError: Invalid B{C{radius}} or B{C{options}}.
674 @see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_}
675 and L{pygeodesy.nearestOn5}.
676 '''
677 if _LatLon_ in options:
678 raise _ValueError(**options)
680 lat, lon, d, c, h = _nearestOn5(self, points, closed=closed, **options)
681 return NearestOn3Tuple(self.classof(lat, lon, height=h),
682 degrees2m(d, radius=radius), c)
684 def toCartesian(self, **Cartesian_datum_kwds): # PYCHOK Cartesian=Cartesian, datum=None
685 '''Convert this point to C{Karney}-based cartesian (ECEF)
686 coordinates.
688 @kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}}
689 and other keyword arguments, ignored
690 if C{B{Cartesian} is None}. Use
691 C{B{Cartesian}=...} to override
692 this L{Cartesian} class or specify
693 C{B{Cartesian}=None}.
695 @return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}}
696 is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height,
697 C, M, datum)} with C{C} and C{M} if available.
699 @raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument.
700 '''
701 kwds = _xkwds(Cartesian_datum_kwds, Cartesian=Cartesian, datum=self.datum)
702 return LatLonSphericalBase.toCartesian(self, **kwds)
704 def _trackDistanceTo3(self, start, end, radius, wrap):
705 '''(INTERNAL) Helper for .along-/crossTrackDistanceTo.
706 '''
707 self.others(start=start)
708 self.others(end=end)
710 r = Radius_(radius)
711 r = start.distanceTo(self, r, wrap=wrap) / r
713 b = radians(start.initialBearingTo(self, wrap=wrap))
714 e = radians(start.initialBearingTo(end, wrap=wrap))
715 x = asin(sin(r) * sin(b - e))
716 return r, x, (e - b)
718 def triangle7(self, otherB, otherC, radius=R_M, wrap=False):
719 '''Compute the angles, sides and area of a spherical triangle.
721 @arg otherB: Second triangle point (C{LatLon}).
722 @arg otherC: Third triangle point (C{LatLon}).
723 @kwarg radius: Mean earth radius, ellipsoid or datum
724 (C{meter}, L{Ellipsoid}, L{Ellipsoid2},
725 L{Datum} or L{a_f2Tuple}) or C{None}.
726 @kwarg wrap: Wrap/unroll angular distances (C{bool}).
728 @return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if
729 B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A,
730 a, B, b, C, c, D, E)}.
732 @see: Function L{triangle7} and U{Spherical trigonometry
733 <https://WikiPedia.org/wiki/Spherical_trigonometry>}.
734 '''
735 self.others(otherB=otherB)
736 self.others(otherC=otherC)
738 r = self.philam + otherB.philam + otherC.philam
739 t = triangle8_(*r, wrap=wrap)
740 return self._xnamed(_t7Tuple(t, radius))
742 def trilaterate5(self, distance1, point2, distance2, point3, distance3,
743 area=True, eps=EPS1, radius=R_M, wrap=False):
744 '''Trilaterate three points by area overlap or perimeter intersection
745 of three corresponding circles.
747 @arg distance1: Distance to this point (C{meter}, same units
748 as B{C{radius}}).
749 @arg point2: Second center point (C{LatLon}).
750 @arg distance2: Distance to point2 (C{meter}, same units as
751 B{C{radius}}).
752 @arg point3: Third center point (C{LatLon}).
753 @arg distance3: Distance to point3 (C{meter}, same units as
754 B{C{radius}}).
755 @kwarg area: If C{True} compute the area overlap, otherwise the
756 perimeter intersection of the circles (C{bool}).
757 @kwarg eps: The required I{minimal overlap} for C{B{area}=True}
758 or the I{intersection margin} for C{B{area}=False}
759 (C{meter}, same units as B{C{radius}}).
760 @kwarg radius: Mean earth radius (C{meter}, conventionally).
761 @kwarg wrap: Wrap/unroll angular distances (C{bool}).
763 @return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)}
764 with C{min} and C{max} in C{meter}, same units as B{C{eps}},
765 the corresponding trilaterated points C{minPoint} and
766 C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number
767 of trilatered points found for the given B{C{eps}}.
769 If only a single trilaterated point is found, C{min I{is}
770 max}, C{minPoint I{is} maxPoint} and C{n = 1}.
772 For C{B{area}=True}, C{min} and C{max} are the smallest
773 respectively largest I{radial} overlap found.
775 For C{B{area}=False}, C{min} and C{max} represent the
776 nearest respectively farthest intersection margin.
778 If C{B{area}=True} and all 3 circles are concentric, C{n =
779 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}}
780 with the smallest B{C{distance#}} C{min} and C{max} the
781 largest B{C{distance#}}.
783 @raise IntersectionError: Trilateration failed for the given B{C{eps}},
784 insufficient overlap for C{B{area}=True} or
785 no intersection or all (near-)concentric for
786 C{B{area}=False}.
788 @raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
790 @raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}},
791 B{C{distance2}}, B{C{distance3}} or B{C{radius}}.
792 '''
793 return _trilaterate5(self, distance1,
794 self.others(point2=point2), distance2,
795 self.others(point3=point3), distance3,
796 area=area, radius=radius, eps=eps, wrap=wrap)
799_T00 = LatLon(0, 0, name='T00') # reference instance (L{LatLon})
802def areaOf(points, radius=R_M, wrap=True):
803 '''Calculate the area of a (spherical) polygon or composite
804 (with the pointsjoined by great circle arcs).
806 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
807 or L{BooleanGH}).
808 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
809 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
810 or C{None}.
811 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
813 @return: Polygon area (C{meter} I{quared}, same units as B{C{radius}}
814 or C{radians} if B{C{radius}} is C{None}).
816 @raise PointsError: Insufficient number of B{C{points}}.
818 @raise TypeError: Some B{C{points}} are not L{LatLon}.
820 @raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
822 @note: The area is based on I{Karney}'s U{'Area of a spherical
823 polygon'<https://MathOverflow.net/questions/97711/
824 the-area-of-spherical-polygons>}, 3rd Answer.
826 @see: Functions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf},
827 L{pygeodesy.excessKarney}, L{ellipsoidalExact.areaOf} and
828 L{ellipsoidalKarney.areaOf}.
830 @example:
832 >>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1)
833 >>> areaOf(b) # 8666058750.718977
835 >>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1)
836 >>> areaOf(c) # 6.18e9
837 '''
838 if _MODS.booleans.isBoolean(points):
839 return points._sum2(LatLon, areaOf, radius=radius, wrap=wrap)
841 Ps = _T00.PointsIter(points, loop=1)
842 p1 = p2 = Ps[0]
843 a1, b1 = p1.philam
844 ta1, z1 = tan_2(a1), None
846 A = Fsum() # mean phi
847 E = Fsum() # see L{pygeodesy.excessKarney_}
848 # ispolar: Summation of course deltas around pole is 0° rather than normally ±360°
849 # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html>
850 # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
851 D = Fsum()
852 for i, p2 in Ps.enumerate(closed=True):
853 a2, b2 = p2.philam
854 db, b2 = unrollPI(b1, b2, wrap=wrap if i else False)
855 ta2 = tan_2(a2)
856 A += a2
857 E += atan2(tan_2(db, points=i) * (ta1 + ta2),
858 _1_0 + ta1 * ta2)
859 ta1, b1 = ta2, b2
861 if not p2.isequalTo(p1, eps=EPS):
862 z, z2 = _bearingTo2(p1, p2, wrap=wrap)
863 if z1 is not None:
864 D += wrap180(z - z1) # (z - z1 + 540) ...
865 D += wrap180(z2 - z) # (z2 - z + 540) % 360 - 180
866 p1, z1 = p2, z2
868 R = fabs(E * _2_0)
869 if fabs(D) < _90_0: # ispolar(points)
870 R = fabs(R - PI2)
871 if radius:
872 a = degrees(A.fover(len(A))) # mean lat
873 R *= _mean_radius(radius, a)**2
874 return float(R)
877def _destination2(a, b, r, t):
878 '''(INTERNAL) Destination lat- and longitude in C{radians}.
880 @arg a: Latitude (C{radians}).
881 @arg b: Longitude (C{radians}).
882 @arg r: Angular distance (C{radians}).
883 @arg t: Bearing (compass C{radians}).
885 @return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
886 '''
887 # see <https://www.EdWilliams.org/avform.htm#LL>
888 sa, ca, sr, cr, st, ct = sincos2_(a, r, t)
889 ca *= sr
891 a = asin1(ct * ca + cr * sa)
892 d = atan2(st * ca, cr - sa * sin(a))
893 # note, in EdWilliams.org/avform.htm W is + and E is -
894 return a, (b + d) # (mod(b + d + PI, PI2) - PI)
897def _int3d2(start, end, wrap, _i_, Vector, hs):
898 # see <https://www.EdWilliams.org/intersect.htm> (5) ff
899 # and similar logic in .ellipsoidalBaseDI._intersect3
900 a1, b1 = start.philam
902 if isscalar(end): # bearing, get pseudo-end point
903 a2, b2 = _destination2(a1, b1, PI_4, radians(end))
904 else: # must be a point
905 start.others(end, name=_end_ + _i_)
906 hs.append(end.height)
907 a2, b2 = end.philam
909 db, b2 = unrollPI(b1, b2, wrap=wrap)
910 if max(fabs(db), fabs(a2 - a1)) < EPS:
911 raise _ValueError(_SPACE_(_line_ + _i_, _null_))
912 # note, in EdWilliams.org/avform.htm W is + and E is -
913 sb21, cb21, sb12, cb12 = sincos2_(db * _0_5,
914 -(b1 + b2) * _0_5)
915 cb21 *= sin(a1 - a2) # sa21
916 sb21 *= sin(a1 + a2) # sa12
917 x = Vector(sb12 * cb21 - cb12 * sb21,
918 cb12 * cb21 + sb12 * sb21,
919 cos(a1) * cos(a2) * sin(db)) # ll=start
920 return x.unit(), (db, (a2 - a1)) # negated d
923def _intdot(ds, a1, b1, a, b, wrap):
924 # compute dot product ds . (-b + b1, a - a1)
925 db, _ = unrollPI(b1, b, wrap=wrap)
926 return fdot(ds, db, a - a1)
929def intersecant2(center, circle, point, bearing, radius=R_M, exact=False,
930 height=None, wrap=True):
931 '''Compute the intersections of a circle and a line.
933 @arg center: Center of the circle (L{LatLon}).
934 @arg circle: Radius of the circle (C{meter}, same units as B{C{radius}})
935 or a point on the circle (L{LatLon}).
936 @arg point: A point in- or outside the circle (L{LatLon}).
937 @arg bearing: Bearing at the B{C{point}} (compass C{degrees360}) or
938 a second point on the line (L{LatLon}).
939 @kwarg radius: Mean earth radius (C{meter}, conventionally).
940 @kwarg exact: If C{True} use the I{exact} rhumb methods for azimuth,
941 destination and distance, if C{False} use the basic
942 rhumb methods (C{bool}) or if C{None} use the I{great
943 circle} methods.
944 @kwarg height: Optional height for the intersection points (C{meter},
945 conventionally) or C{None}.
946 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
948 @return: 2-Tuple of the intersection points (representing a chord),
949 each an instance of this class. For a tangent line, each
950 point C{is} this very instance.
952 @raise IntersectionError: The circle and line do not intersect.
954 @raise TypeError: If B{C{center}} or B{C{point}} not L{LatLon} or
955 B{C{circle}} or B{C{bearing}} invalid.
957 @raise ValueError: Invalid B{C{circle}}, B{C{bearing}}, B{C{radius}},
958 B{C{exact}} or B{C{height}}.
959 '''
960 c = _T00.others(center=center)
961 p = _T00.others(point=point)
962 try:
963 return _intersecant2(c, circle, p, bearing, radius=radius, exact=exact,
964 height=height, wrap=wrap)
965 except (TypeError, ValueError) as x:
966 raise _xError(x, center=center, circle=circle, point=point, bearing=bearing, exact=exact)
969def _intersect(start1, end1, start2, end2, height=None, wrap=False, # in.ellipsoidalBaseDI._intersect3
970 LatLon=None, **LatLon_kwds):
971 # (INTERNAL) Intersect two (spherical) lines, see L{intersection}
972 # above, separated to allow callers to embellish any exceptions
974 hs = [start1.height, start2.height]
976 a1, b1 = start1.philam
977 a2, b2 = start2.philam
979 db, b2 = unrollPI(b1, b2, wrap=wrap)
980 r12 = vincentys_(a2, a1, db)
981 if fabs(r12) < EPS: # [nearly] coincident points
982 a, b = favg(a1, a2), favg(b1, b2)
984 # see <https://www.EdWilliams.org/avform.htm#Intersection>
985 elif isscalar(end1) and isscalar(end2): # both bearings
986 sa1, ca1, sa2, ca2, sr12, cr12 = sincos2_(a1, a2, r12)
988 x1, x2 = (sr12 * ca1), (sr12 * ca2)
989 if isnear0(x1) or isnear0(x2):
990 raise IntersectionError(_parallel_)
991 # handle domain error for equivalent longitudes,
992 # see also functions asin_safe and acos_safe at
993 # <https://www.EdWilliams.org/avform.htm#Math>
994 t1, t2 = acos1((sa2 - sa1 * cr12) / x1), \
995 acos1((sa1 - sa2 * cr12) / x2)
996 if sin(db) > 0:
997 t12, t21 = t1, PI2 - t2
998 else:
999 t12, t21 = PI2 - t1, t2
1000 t13, t23 = radiansPI2(end1), radiansPI2(end2)
1001 sx1, cx1, sx2, cx2 = sincos2_(wrapPI(t13 - t12), # angle 2-1-3
1002 wrapPI(t21 - t23)) # angle 1-2-3)
1003 if isnear0(sx1) and isnear0(sx2):
1004 raise IntersectionError(_infinite_)
1005 sx3 = sx1 * sx2
1006# XXX if sx3 < 0:
1007# XXX raise ValueError(_ambiguous_)
1008 x3 = acos1(cr12 * sx3 - cx2 * cx1)
1009 r13 = atan2(sr12 * sx3, cx2 + cx1 * cos(x3))
1011 a, b = _destination2(a1, b1, r13, t13)
1012 # like .ellipsoidalBaseDI,_intersect3, if this intersection
1013 # is "before" the first point, use the antipodal intersection
1014 if opposing_(t13, bearing_(a1, b1, a, b, wrap=wrap)):
1015 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1017 else: # end point(s) or bearing(s)
1018 _N_vector_ = _MODS.nvectorBase._N_vector_
1020 x1, d1 = _int3d2(start1, end1, wrap, _1_, _N_vector_, hs)
1021 x2, d2 = _int3d2(start2, end2, wrap, _2_, _N_vector_, hs)
1022 x = x1.cross(x2)
1023 if x.length < EPS: # [nearly] colinear or parallel lines
1024 raise IntersectionError(_colinear_)
1025 a, b = x.philam
1026 # choose intersection similar to sphericalNvector
1027 if not (_intdot(d1, a1, b1, a, b, wrap) *
1028 _intdot(d2, a2, b2, a, b, wrap)) > 0:
1029 a, b = antipode_(a, b) # PYCHOK PhiLam2Tuple
1031 h = fmean(hs) if height is None else Height(height)
1032 return _LL3Tuple(degrees90(a), degrees180(b), h,
1033 intersection, LatLon, LatLon_kwds)
1036def intersection(start1, end1, start2, end2, height=None, wrap=False,
1037 LatLon=LatLon, **LatLon_kwds):
1038 '''Compute the intersection point of two lines, each defined
1039 by two points or a start point and bearing from North.
1041 @arg start1: Start point of the first line (L{LatLon}).
1042 @arg end1: End point of the first line (L{LatLon}) or
1043 the initial bearing at the first start point
1044 (compass C{degrees360}).
1045 @arg start2: Start point of the second line (L{LatLon}).
1046 @arg end2: End point of the second line (L{LatLon}) or
1047 the initial bearing at the second start point
1048 (compass C{degrees360}).
1049 @kwarg height: Optional height for the intersection point,
1050 overriding the mean height (C{meter}).
1051 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1052 @kwarg LatLon: Optional class to return the intersection
1053 point (L{LatLon}) or C{None}.
1054 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1055 arguments, ignored if C{B{LatLon} is None}.
1057 @return: The intersection point as a (B{C{LatLon}}) or if
1058 C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon,
1059 height)}. An alternate intersection point might
1060 be the L{antipode} to the returned result.
1062 @raise IntersectionError: Ambiguous or infinite intersection
1063 or colinear, parallel or otherwise
1064 non-intersecting lines.
1066 @raise TypeError: A B{C{start}} or B{C{end}} point not L{LatLon}.
1068 @raise ValueError: Invalid B{C{height}} or C{null} line.
1070 @example:
1072 >>> p = LatLon(51.8853, 0.2545)
1073 >>> s = LatLon(49.0034, 2.5735)
1074 >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E'
1075 '''
1076 s1 = _T00.others(start1=start1)
1077 s2 = _T00.others(start2=start2)
1078 try:
1079 return _intersect(s1, end1, s2, end2, height=height, wrap=wrap,
1080 LatLon=LatLon, **LatLon_kwds)
1081 except (TypeError, ValueError) as x:
1082 raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
1085def intersections2(center1, rad1, center2, rad2, radius=R_M, eps=_0_0,
1086 height=None, wrap=True,
1087 LatLon=LatLon, **LatLon_kwds):
1088 '''Compute the intersection points of two circles each defined
1089 by a center point and a radius.
1091 @arg center1: Center of the first circle (L{LatLon}).
1092 @arg rad1: Radius of the first circle (C{meter} or C{radians},
1093 see B{C{radius}}).
1094 @arg center2: Center of the second circle (L{LatLon}).
1095 @arg rad2: Radius of the second circle (C{meter} or C{radians},
1096 see B{C{radius}}).
1097 @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}},
1098 B{C{rad2}} and B{C{eps}} are given in C{radians}).
1099 @kwarg eps: Required overlap (C{meter} or C{radians}, see
1100 B{C{radius}}).
1101 @kwarg height: Optional height for the intersection points (C{meter},
1102 conventionally) or C{None} for the I{"radical height"}
1103 at the I{radical line} between both centers.
1104 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1105 @kwarg LatLon: Optional class to return the intersection
1106 points (L{LatLon}) or C{None}.
1107 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1108 arguments, ignored if C{B{LatLon} is None}.
1110 @return: 2-Tuple of the intersection points, each a B{C{LatLon}}
1111 instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat,
1112 lon, height)}. For abutting circles, both intersection
1113 points are the same instance, aka the I{radical center}.
1115 @raise IntersectionError: Concentric, antipodal, invalid or
1116 non-intersecting circles.
1118 @raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
1120 @raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}},
1121 B{C{eps}} or B{C{height}}.
1123 @note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
1125 @see: This U{Answer<https://StackOverflow.com/questions/53324667/
1126 find-intersection-coordinates-of-two-circles-on-earth/53331953>}.
1127 '''
1128 c1 = _T00.others(center1=center1)
1129 c2 = _T00.others(center2=center2)
1130 try:
1131 return _intersects2(c1, rad1, c2, rad2, radius=radius, eps=eps,
1132 height=height, wrap=wrap,
1133 LatLon=LatLon, **LatLon_kwds)
1134 except (TypeError, ValueError) as x:
1135 raise _xError(x, center1=center1, rad1=rad1, center2=center2, rad2=rad2)
1138def _intersects2(c1, rad1, c2, rad2, radius=R_M, eps=_0_0, # in .ellipsoidalBaseDI._intersects2
1139 height=None, too_d=None, wrap=True,
1140 LatLon=LatLon, **LatLon_kwds):
1141 # (INTERNAL) Intersect two spherical circles, see L{intersections2}
1142 # above, separated to allow callers to embellish any exceptions
1144 def _dest3(bearing, h):
1145 a, b = _destination2(a1, b1, r1, bearing)
1146 return _LL3Tuple(degrees90(a), degrees180(b), h,
1147 intersections2, LatLon, LatLon_kwds)
1149 r1, r2, f = _rads3(rad1, rad2, radius)
1150 if f: # swapped
1151 c1, c2 = c2, c1 # PYCHOK swap
1153 a1, b1 = c1.philam
1154 a2, b2 = c2.philam
1156 db, b2 = unrollPI(b1, b2, wrap=wrap)
1157 d = vincentys_(a2, a1, db) # radians
1158 if d < max(r1 - r2, EPS):
1159 raise IntersectionError(_near_(_concentric_)) # XXX ConcentricError?
1161 r = eps if radius is None else (m2radians(
1162 eps, radius=radius) if eps else _0_0)
1163 if r < _0_0:
1164 raise _ValueError(eps=r)
1166 x = fsum_(r1, r2, -d) # overlap
1167 if x > max(r, EPS):
1168 sd, cd, sr1, cr1, _, cr2 = sincos2_(d, r1, r2)
1169 x = sd * sr1
1170 if isnear0(x):
1171 raise _ValueError(_invalid_)
1172 x = acos1((cr2 - cd * cr1) / x) # 0 <= x <= PI
1174 elif x < r: # PYCHOK no cover
1175 t = (d * radius) if too_d is None else too_d
1176 raise IntersectionError(_too_(_Fmt.distant(t)))
1178 if height is None: # "radical height"
1179 f = _radical2(d, r1, r2).ratio
1180 h = Height(favg(c1.height, c2.height, f=f))
1181 else:
1182 h = Height(height)
1184 b = bearing_(a1, b1, a2, b2, final=False, wrap=wrap)
1185 if x < EPS4: # externally ...
1186 r = _dest3(b, h)
1187 elif x > _PI_EPS4: # internally ...
1188 r = _dest3(b + PI, h)
1189 else:
1190 return _dest3(b + x, h), _dest3(b - x, h)
1191 return r, r # ... abutting circles
1194@deprecated_function
1195def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover
1196 '''DEPRECATED, use function L{pygeodesy.ispolar}.
1197 '''
1198 return ispolar(points, wrap=wrap)
1201def _LL3Tuple(lat, lon, height, func, LatLon, LatLon_kwds):
1202 '''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}.
1203 '''
1204 n = func.__name__
1205 if LatLon is None:
1206 r = LatLon3Tuple(lat, lon, height, name=n)
1207 else:
1208 kwds = _xkwds(LatLon_kwds, height=height, name=n)
1209 r = LatLon(lat, lon, **kwds)
1210 return r
1213def meanOf(points, height=None, LatLon=LatLon, **LatLon_kwds):
1214 '''Compute the geographic mean of several points.
1216 @arg points: Points to be averaged (L{LatLon}[]).
1217 @kwarg height: Optional height at mean point, overriding the mean
1218 height (C{meter}).
1219 @kwarg LatLon: Optional class to return the mean point (L{LatLon})
1220 or C{None}.
1221 @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword
1222 arguments, ignored if C{B{LatLon} is None}.
1224 @return: The geographic mean and height (B{C{LatLon}}) or a
1225 L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}}
1226 is C{None}.
1228 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1230 @raise ValueError: No B{C{points}} or invalid B{C{height}}.
1231 '''
1232 # geographic mean
1233 m = sumOf(p._N_vector for p in _T00.PointsIter(points).iterate(closed=False))
1234 lat, lon = m._N_vector.latlon
1236 if height is None:
1237 h = fmean(p.height for p in _T00.PointsIter(points).iterate(closed=False))
1238 else: # PYCHOK no cover
1239 h = Height(height)
1240 return _LL3Tuple(lat, lon, h, meanOf, LatLon, LatLon_kwds)
1243@deprecated_function
1244def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover
1245 '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
1247 @return: ... 2-tuple C{(closest, distance)} of the C{closest}
1248 point (L{LatLon}) on the polygon and the C{distance}
1249 between the C{closest} and the given B{C{point}}. The
1250 C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat,
1251 lon)} if B{C{LatLon}} is C{None} ...
1252 '''
1253 ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple
1254 if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None:
1255 ll = LatLon2Tuple(ll.lat, ll.lon)
1256 return ll, d
1259def nearestOn3(point, points, closed=False, radius=R_M,
1260 LatLon=LatLon, **options):
1261 '''Locate the point on a path or polygon closest to a reference point.
1263 Distances are I{approximated} using function L{pygeodesy.equirectangular_},
1264 subject to the supplied B{C{options}}.
1266 @arg point: The reference point (L{LatLon}).
1267 @arg points: The path or polygon points (L{LatLon}[]).
1268 @kwarg closed: Optionally, close the polygon (C{bool}).
1269 @kwarg radius: Mean earth radius (C{meter}).
1270 @kwarg LatLon: Optional class to return the closest point (L{LatLon})
1271 or C{None}.
1272 @kwarg options: Optional keyword arguments for function
1273 L{pygeodesy.equirectangular_}.
1275 @return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the
1276 C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat,
1277 lon, height)} if B{C{LatLon}} is C{None}. The C{distance}
1278 is the L{pygeodesy.equirectangular_} distance between the
1279 C{closest} and the given B{C{point}} converted to C{meter},
1280 same units as B{C{radius}}. The C{angle} from the given
1281 B{C{point}} to the C{closest} is in compass C{degrees360},
1282 like function L{pygeodesy.compassAngle}. The C{height} is
1283 the (interpolated) height at the C{closest} point.
1285 @raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}},
1286 see function L{pygeodesy.equirectangular_}.
1288 @raise PointsError: Insufficient number of B{C{points}}.
1290 @raise TypeError: Some B{C{points}} are not C{LatLon}.
1292 @raise ValueError: Invalid B{C{radius}}.
1294 @see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}.
1295 '''
1296 lat, lon, d, c, h = _nearestOn5(point, points, closed=closed,
1297 LatLon=None, **options)
1298 d = degrees2m(d, radius=radius)
1299 n = nearestOn3.__name__
1300 r = LatLon3Tuple(lat, lon, h, name=n) if LatLon is None else \
1301 LatLon(lat, lon, height=h, name=n)
1302 return NearestOn3Tuple(r, d, c, name=n)
1305def perimeterOf(points, closed=False, radius=R_M, wrap=True):
1306 '''Compute the perimeter of a (spherical) polygon or composite
1307 (with great circle arcs joining the points).
1309 @arg points: The polygon points or clips (L{LatLon}[], L{BooleanFHP}
1310 or L{BooleanGH}).
1311 @kwarg closed: Optionally, close the polygon (C{bool}).
1312 @kwarg radius: Mean earth radius (C{meter}) or C{None}.
1313 @kwarg wrap: Wrap and unroll longitudes (C{bool}).
1315 @return: Polygon perimeter (C{meter}, same units as B{C{radius}}
1316 or C{radians} if B{C{radius}} is C{None}).
1318 @raise PointsError: Insufficient number of B{C{points}}.
1320 @raise TypeError: Some B{C{points}} are not L{LatLon}.
1322 @raise ValueError: Invalid B{C{radius}} or C{B{closed}=False} with
1323 C{B{points}} a composite.
1325 @note: Distances are based on function L{pygeodesy.vincentys_}.
1327 @see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf}
1328 and L{ellipsoidalKarney.perimeterOf}.
1329 '''
1330 def _rads(points, closed, wrap): # angular edge lengths in radians
1331 Ps = _T00.PointsIter(points, loop=1)
1332 a1, b1 = Ps[0].philam
1333 for i, p in Ps.enumerate(closed=closed):
1334 a2, b2 = p.philam
1335 db, b2 = unrollPI(b1, b2, wrap=wrap if i else False)
1336 yield vincentys_(a2, a1, db)
1337 a1, b1 = a2, b2
1339 if _MODS.booleans.isBoolean(points):
1340 if not closed:
1341 raise _ValueError(closed=closed, points=_composite_)
1342 r = points._sum2(LatLon, perimeterOf, closed=True, radius=radius, wrap=wrap)
1343 else:
1344 r = fsum(_rads(points, closed, wrap), floats=True)
1345 return _r2m(r, radius)
1348def _r2m(r, radius):
1349 '''(INTERNAL) Angular distance in C{radians} to C{meter}.
1350 '''
1351 if radius is not None: # not in (None, _0_0)
1352 r *= R_M if radius is R_M else Radius(radius)
1353 return r
1356def triangle7(latA, lonA, latB, lonB, latC, lonC, radius=R_M,
1357 excess=excessAbc_,
1358 wrap=False):
1359 '''Compute the angles, sides, and area of a (spherical) triangle.
1361 @arg latA: First corner latitude (C{degrees}).
1362 @arg lonA: First corner longitude (C{degrees}).
1363 @arg latB: Second corner latitude (C{degrees}).
1364 @arg lonB: Second corner longitude (C{degrees}).
1365 @arg latC: Third corner latitude (C{degrees}).
1366 @arg lonC: Third corner longitude (C{degrees}).
1367 @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter},
1368 L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple})
1369 or C{None}.
1370 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1371 L{excessGirard_} or L{excessLHuilier_}).
1372 @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
1374 @return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with
1375 spherical angles C{A}, C{B} and C{C}, angular sides
1376 C{a}, C{b} and C{c} all in C{degrees} and C{area}
1377 in I{square} C{meter} or same units as B{C{radius}}
1378 I{squared} or if C{B{radius}=0} or C{None}, a
1379 L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in
1380 C{radians} with the I{spherical excess} C{E} as the
1381 C{unit area} in C{radians}.
1382 '''
1383 t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA),
1384 Phi_(latB=latB), Lam_(lonB=lonB),
1385 Phi_(latC=latC), Lam_(lonC=lonC),
1386 excess=excess, wrap=wrap)
1387 return _t7Tuple(t, radius)
1390def triangle8_(phiA, lamA, phiB, lamB, phiC, lamC, excess=excessAbc_,
1391 wrap=False):
1392 '''Compute the angles, sides, I{spherical deficit} and I{spherical
1393 excess} of a (spherical) triangle.
1395 @arg phiA: First corner latitude (C{radians}).
1396 @arg lamA: First corner longitude (C{radians}).
1397 @arg phiB: Second corner latitude (C{radians}).
1398 @arg lamB: Second corner longitude (C{radians}).
1399 @arg phiC: Third corner latitude (C{radians}).
1400 @arg lamC: Third corner longitude (C{radians}).
1401 @kwarg excess: I{Spherical excess} callable (L{excessAbc_},
1402 L{excessGirard_} or L{excessLHuilier_}).
1403 @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
1405 @return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with
1406 spherical angles C{A}, C{B} and C{C}, angular sides
1407 C{a}, C{b} and C{c}, I{spherical deficit} C{D} and
1408 I{spherical excess} C{E}, all in C{radians}.
1409 '''
1410 def _a_r(w, phiA, lamA, phiB, lamB, phiC, lamC):
1411 d, _ = unrollPI(lamB, lamC, wrap=w)
1412 a = vincentys_(phiC, phiB, d)
1413 return a, (phiB, lamB, phiC, lamC, phiA, lamA)
1415 def _A_r(a, sa, ca, sb, cb, sc, cc):
1416 s = sb * sc
1417 A = acos1((ca - cb * cc) / s) if isnon0(s) else a
1418 return A, (sb, cb, sc, cc, sa, ca) # rotate sincos2's
1420 # notation: side C{a} is oposite to corner C{A}, etc.
1421 a, r = _a_r(wrap, phiA, lamA, phiB, lamB, phiC, lamC)
1422 b, r = _a_r(wrap, *r)
1423 c, _ = _a_r(wrap, *r)
1425 A, r = _A_r(a, *sincos2_(a, b, c))
1426 B, r = _A_r(b, *r)
1427 C, _ = _A_r(c, *r)
1429 D = fsum_(PI2, -a, -b, -c, floats=True) # deficit aka defect
1430 E = excessGirard_(A, B, C) if excess in (excessGirard_, True) else (
1431 excessLHuilier_(a, b, c) if excess in (excessLHuilier_, False) else
1432 excessAbc_(*max((A, b, c), (B, c, a), (C, a, b))))
1434 return Triangle8Tuple(A, a, B, b, C, c, D, E)
1437def _t7Tuple(t, radius):
1438 '''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}.
1439 '''
1440 if radius: # not in (None, _0_0)
1441 r = radius if isscalar(radius) else \
1442 _ellipsoidal_datum(radius).ellipsoid.Rmean
1443 A, B, C = map1(degrees, t.A, t.B, t.C)
1444 t = Triangle7Tuple(A, (r * t.a),
1445 B, (r * t.b),
1446 C, (r * t.c), t.E * r**2)
1447 return t
1450__all__ += _ALL_OTHER(Cartesian, LatLon, # classes
1451 areaOf, # functions
1452 intersecant2, intersection, intersections2, ispolar,
1453 isPoleEnclosedBy, # DEPRECATED, use ispolar
1454 meanOf,
1455 nearestOn2, nearestOn3,
1456 perimeterOf,
1457 sumOf, # == vector3d.sumOf
1458 triangle7, triangle8_)
1460# **) MIT License
1461#
1462# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1463#
1464# Permission is hereby granted, free of charge, to any person obtaining a
1465# copy of this software and associated documentation files (the "Software"),
1466# to deal in the Software without restriction, including without limitation
1467# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1468# and/or sell copies of the Software, and to permit persons to whom the
1469# Software is furnished to do so, subject to the following conditions:
1470#
1471# The above copyright notice and this permission notice shall be included
1472# in all copies or substantial portions of the Software.
1473#
1474# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1475# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1476# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1477# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1478# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1479# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1480# OTHER DEALINGS IN THE SOFTWARE.