Coverage for pygeodesy/sphericalTrigonometry.py : 97%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
Trigonometric classes geodetic (lat-/longitude) L{LatLon} and geocentric (ECEF) L{Cartesian} and functions L{areaOf}, L{intersection}, L{intersections2}, L{isPoleEnclosedBy}, L{meanOf}, L{nearestOn3} and L{perimeterOf}, I{all spherical}.
Pure Python implementation of geodetic (lat-/longitude) methods using spherical trigonometry, transcoded from JavaScript originals by I{(C) Chris Veness 2011-2016} published under the same MIT Licence**, see U{Latitude/Longitude<https://www.Movable-Type.co.UK/scripts/latlong.html>}. ''' # make sure int/int division yields float quotient, see .basics
map1, signOf _ValueError, IntersectionError, _xError, \ _xkwds, _xkwds_get excessGirard, excessLHuilier, opposing_, _radical2, \ vincentys_ _coincident_, _colinear_, _concentric_, _convex_, \ _end_, _infinite_, _invalid_, _LatLon_, _near_, \ _not_, _null_, _points_, _SPACE_, _too_, _1_, _2_, \ _0_0, _0_5, _1_0, _2_0, _90_0 # from pygeodesy.named import notImplemented # from .points NearestOn3Tuple, Triangle7Tuple, \ Triangle8Tuple notImplemented, Fmt as _Fmt # XXX shadowed LatLonSphericalBase, _rads3, _trilaterate5 # from pygeodesy.streprs import Fmt as _Fmt # from .points XXX shadowed Radius_, Scalar m2radians, radiansPI2, sincos2_, tan_2, unrollPI, \ wrap180, wrapPI
raise _AssertionError(EPS4=EPS4, PI=PI, PI_EPS4=_PI_EPS4)
'''(INTERNAL) Destination lat- and longitude in C{radians}.
@arg a: Latitude (C{radians}). @arg b: Longitude (C{radians}). @arg r: Angular distance (C{radians}). @arg t: Bearing (compass C{radians}).
@return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}). ''' # see <https://www.EdWilliams.org/avform.htm#LL>
# note, in EdWilliams.org/avform.htm W is + and E is -
'''(INTERNAL) Angular distance in C{radians} to C{meter}. '''
'''Extended to convert geocentric, L{Cartesian} points to spherical, geodetic L{LatLon}. '''
'''Convert this cartesian point to a C{spherical} geodetic point.
@kwarg LatLon_and_kwds: Optional L{LatLon} and L{LatLon} keyword arguments. Use C{B{LatLon}=...} to override this L{LatLon} class or specify C{B{LatLon}=None}.
@return: The geodetic point (L{LatLon}) or if B{C{LatLon}} is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with C{C} and C{M} if available.
@raise TypeError: Invalid B{C{LatLon_and_kwds}} argument. '''
'''New point on spherical model earth model.
@example:
>>> p = LatLon(52.205, 0.119) # height=0 '''
'''Compute the (angular) distance (signed) from the start to the closest point on the great circle path defined by a start and an end point.
That is, if a perpendicular is drawn from this point to the great circle path, the along-track distance is the distance from the start point to the point where the perpendicular crosses the path.
@arg start: Start point of great circle path (L{LatLon}). @arg end: End point of great circle path (L{LatLon}). @kwarg radius: Mean earth radius (C{meter}) or C{None}. @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: Distance along the great circle path (C{meter}, same units as B{C{radius}}) or C{radians} if C{B{radius} is None}, positive if after the B{C{start}} toward the B{C{end}} point of the path or negative if before the B{C{start}} point.
@raise TypeError: Invalid B{C{start}} or B{C{end}} point.
@raise ValueError: Invalid B{C{radius}}.
@example:
>>> p = LatLon(53.2611, -0.7972)
>>> s = LatLon(53.3206, -1.7297) >>> e = LatLon(53.1887, 0.1334) >>> d = p.alongTrackDistanceTo(s, e) # 62331.58 ''' _r2m(copysign0(acos1(cos(r) / cx), cos(b)), radius)
def bearingTo(self, other, wrap=False, raiser=False): # PYCHOK no cover '''DEPRECATED, use method L{initialBearingTo}. ''' return self.initialBearingTo(other, wrap=wrap, raiser=raiser)
'''Return the pair of meridians at which a great circle defined by this and an other point crosses the given latitude.
@arg other: The other point defining great circle (L{LatLon}). @arg lat: Latitude at the crossing (C{degrees}). @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or C{None} if the great circle doesn't reach B{C{lat}}. '''
sa2, ca2, sdb, cdb = sincos2_(a, a1, a2, db)
if h < EPS or abs(z) > h: # PYCHOK no cover return None # great circle doesn't reach latitude
'''Compute the (angular) distance (signed) from this point to the great circle defined by a start and an end point.
@arg start: Start point of great circle path (L{LatLon}). @arg end: End point of great circle path (L{LatLon}). @kwarg radius: Mean earth radius (C{meter}) or C{None}. @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: Distance to great circle (negative if to the left or positive if to the right of the path) (C{meter}, same units as B{C{radius}} or C{radians} if B{C{radius}} is C{None}).
@raise TypeError: The B{C{start}} or B{C{end}} point is not L{LatLon}.
@raise ValueError: Invalid B{C{radius}}.
@example:
>>> p = LatLon(53.2611, -0.7972)
>>> s = LatLon(53.3206, -1.7297) >>> e = LatLon(53.1887, 0.1334) >>> d = p.crossTrackDistanceTo(s, e) # -307.5 '''
'''Locate the destination from this point after having travelled the given distance on the given initial bearing.
@arg distance: Distance travelled (C{meter}, same units as B{C{radius}}). @arg bearing: Bearing from this point (compass C{degrees360}). @kwarg radius: Mean earth radius (C{meter}). @kwarg height: Optional height at destination (C{meter}, same units a B{C{radius}}).
@return: Destination point (L{LatLon}).
@raise ValueError: Invalid B{C{distance}}, B{C{bearing}}, B{C{radius}} or B{C{height}}.
@example:
>>> p1 = LatLon(51.4778, -0.0015) >>> p2 = p1.destination(7794, 300.7) >>> p2.toStr() # '51.5135°N, 000.0983°W'
@JSname: I{destinationPoint}. '''
'''Compute the (angular) distance from this to an other point.
@arg other: The other point (L{LatLon}). @kwarg radius: Mean earth radius (C{meter}) or C{None}. @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: Distance between this and the B{C{other}} point (C{meter}, same units as B{C{radius}} or C{radians} if B{C{radius}} is C{None}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@raise ValueError: Invalid B{C{radius}}.
@example:
>>> p1 = LatLon(52.205, 0.119) >>> p2 = LatLon(48.857, 2.351); >>> d = p1.distanceTo(p2) # 404300 '''
# @Property_RO # def Ecef(self): # '''Get the ECEF I{class} (L{EcefVeness}), I{lazily}. # ''' # return _MODS.ecef.EcefKarney
'''Compute the vector normal to great circle obtained by heading on the given initial bearing from this point.
Direction of vector is such that initial bearing vector b = c × n, where n is an n-vector representing this point.
@arg bearing: Bearing from this point (compass C{degrees360}).
@return: Vector representing great circle (L{Vector3d}).
@raise ValueError: Invalid B{C{bearing}}.
@example:
>>> p = LatLon(53.3206, -1.7297) >>> g = p.greatCircle(96.0) >>> g.toStr() # (-0.794, 0.129, 0.594) '''
-cb * ct - sb * sa * st, ca * st) # XXX .unit()?
'''Compute the initial bearing (forward azimuth) from this to an other point.
@arg other: The other point (spherical L{LatLon}). @kwarg wrap: Wrap and unroll longitudes (C{bool}). @kwarg raiser: Optionally, raise L{CrossError} (C{bool}), use C{B{raiser}=True} for behavior like C{sphericalNvector.LatLon.initialBearingTo}.
@return: Initial bearing (compass C{degrees360}).
@raise CrossError: If this and the B{C{other}} point coincide, provided both B{C{raiser}} is C{True} and L{pygeodesy.crosserrors} is C{True}.
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@example:
>>> p1 = LatLon(52.205, 0.119) >>> p2 = LatLon(48.857, 2.351) >>> b = p1.initialBearingTo(p2) # 156.2
@JSname: I{bearingTo}. '''
# XXX behavior like sphericalNvector.LatLon.initialBearingTo raise CrossError(_points_, self, txt=_coincident_)
'''Locate the point at given fraction between (or along) this and an other point.
@arg other: The other point (L{LatLon}). @arg fraction: Fraction between both points (C{scalar}, 0.0 at this and 1.0 at the other point). @kwarg height: Optional height, overriding the intermediate height (C{meter}). @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: Intermediate point (L{LatLon}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@raise ValueError: Invalid B{C{fraction}} or B{C{height}}.
@example:
>>> p1 = LatLon(52.205, 0.119) >>> p2 = LatLon(48.857, 2.351) >>> p = p1.intermediateTo(p2, 0.25) # 51.3721°N, 000.7073°E
@JSname: I{intermediatePointTo}. ''' if isnear0(f): # PYCHOK no cover r = self elif isnear1(f) and not wrap: # PYCHOK no cover r = self.others(other) else:
sb1, cb1, sb2, cb2 = sincos2_(a1, a2, b1, b2)
else: # PYCHOK no cover a = favg(a1, a2, f=f) # coincident b = favg(b1, b2, f=f)
'''Compute the intersection point of two paths, each defined by two points or a start point and bearing from North.
@arg end1: End point of this path (L{LatLon}) or the initial bearing at this point (compass C{degrees360}). @arg other: Start point of the other path (L{LatLon}). @arg end2: End point of the other path (L{LatLon}) or the initial bearing at the other start point (compass C{degrees360}). @kwarg height: Optional height for intersection point, overriding the mean height (C{meter}). @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: The intersection point (L{LatLon}). An alternate intersection point might be the L{antipode} to the returned result.
@raise IntersectionError: Ambiguous or infinite intersection or colinear, parallel or otherwise non-intersecting paths.
@raise TypeError: If B{C{other}} is not L{LatLon} or B{C{end1}} or B{C{end2}} not C{scalar} nor L{LatLon}.
@raise ValueError: Invalid B{C{height}} or C{null} path.
@example:
>>> p = LatLon(51.8853, 0.2545) >>> s = LatLon(49.0034, 2.5735) >>> i = p.intersection(108.547, s, 32.435) # '50.9078°N, 004.5084°E' ''' LatLon=self.classof) except (TypeError, ValueError) as x: raise _xError(x, start1=self, end1=end1, other=other, end2=end2)
height=None, wrap=True): '''Compute the intersection points of two circles, each defined by a center point and radius.
@arg rad1: Radius of the this circle (C{meter} or C{radians}, see B{C{radius}}). @arg other: Center point of the other circle (L{LatLon}). @arg rad2: Radius of the other circle (C{meter} or C{radians}, see B{C{radius}}). @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}}, B{C{rad2}} and B{C{eps}} are given in C{radians}). @kwarg eps: Required overlap (C{meter} or C{radians}, see B{C{radius}}). @kwarg height: Optional height for the intersection points (C{meter}, conventionally) or C{None} for the I{"radical height"} at the I{radical line} between both centers. @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: 2-Tuple of the intersection points, each a L{LatLon} instance. For abutting circles, both intersection points are the same instance, aka I{radical center}.
@raise IntersectionError: Concentric, antipodal, invalid or non-intersecting circles.
@raise TypeError: If B{C{other}} is not L{LatLon}.
@raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}}, B{C{eps}} or B{C{height}}. ''' height=height, wrap=wrap, LatLon=self.classof) except (TypeError, ValueError) as x: raise _xError(x, center=self, rad1=rad1, other=other, rad2=rad2)
'''Check whether a (convex) polygon encloses this point.
@arg points: The polygon points (L{LatLon}[]).
@return: C{True} if the polygon encloses this point, C{False} otherwise.
@raise PointsError: Insufficient number of B{C{points}}.
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@raise ValueError: Invalid B{C{points}}, non-convex polygon.
@see: Functions L{pygeodesy.isconvex}, L{pygeodesy.isenclosedBy} and L{pygeodesy.ispolar} especially if the B{C{points}} may enclose a pole or wrap around the earth longitudinally.
@example:
>>> b = LatLon(45,1), LatLon(45,2), LatLon(46,2), LatLon(46,1) >>> p = LatLon(45,1, 1.1) >>> inside = p.isEnclosedBy(b) # True '''
# check whether this point on same side of all # polygon edges (to the left or right depending # on anti-/clockwise polygon direction) # get great-circle vector for each edge
return False # outside
# check for convex polygon: angle between # gc vectors, signed by direction of n0 # (otherwise the test above is not reliable) else: raise _ValueError(t, p, txt=_not_(_convex_))
def isEnclosedBy(self, points): # PYCHOK no cover '''DEPRECATED, use method C{isenclosedBy}.''' return self.isenclosedBy(points)
'''Find the midpoint between this and an other point.
@arg other: The other point (L{LatLon}). @kwarg height: Optional height for midpoint, overriding the mean height (C{meter}). @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: Midpoint (L{LatLon}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@raise ValueError: Invalid B{C{height}}.
@example:
>>> p1 = LatLon(52.205, 0.119) >>> p2 = LatLon(48.857, 2.351) >>> m = p1.midpointTo(p2) # '50.5363°N, 001.2746°E' '''
# see <https://MathForum.org/library/drmath/view/51822.html>
'''Locate the point between two points closest to this point.
Distances are approximated by function L{pygeodesy.equirectangular_}, subject to the supplied B{C{options}}.
@arg point1: Start point (L{LatLon}). @arg point2: End point (L{LatLon}). @kwarg radius: Mean earth radius (C{meter}). @kwarg options: Optional keyword arguments for function L{pygeodesy.equirectangular_}.
@return: Closest point on the great circle path (L{LatLon}).
@raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}}, see function L{pygeodesy.equirectangular_}.
@raise NotImplementedError: Keyword argument C{B{within}=False} is not (yet) supported.
@raise TypeError: Invalid B{C{point1}} or B{C{point2}}.
@raise ValueError: Invalid B{C{radius}} or B{C{options}}.
@see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5} and method L{sphericalTrigonometry.LatLon.nearestOn3}. ''' notImplemented(self, within=within)
# # UNTESTED - handle C{B{within}=False} and C{B{within}=True} # wrap = options.get('wrap', False) # a = self.alongTrackDistanceTo(point1, point2, radius=radius, wrap=wrap) # if abs(a) < EPS or (within and a < EPS): # return point1 # d = point1.distanceTo(point2, radius=radius, wrap=wrap) # if isnear0(d): # return point1 # or point2 # elif abs(d - a) < EPS or (a + EPS) > d: # return point2 # f = a / d # if within: # if f > EPS1: # return point2 # elif f < EPS: # return point1 # return point1.intermediateTo(point2, f, wrap=wrap)
except KeyError: pass # without kwarg B{C{within}}, use backward compatible .nearestOn3 **options)[0]
def nearestOn2(self, points, closed=False, radius=R_M, **options): # PYCHOK no cover '''DEPRECATED, use method L{sphericalTrigonometry.LatLon.nearestOn3}.
@return: ... 2-Tuple C{(closest, distance)} of the closest point (L{LatLon}) on the polygon and the distance to that point from this point in C{meter}, same units of B{C{radius}}. ''' r = self.nearestOn3(points, closed=closed, radius=radius, **options) return r.closest, r.distance
'''Locate the point on a polygon closest to this point.
Distances are approximated by function L{pygeodesy.equirectangular_}, subject to the supplied B{C{options}}.
@arg points: The polygon points (L{LatLon}[]). @kwarg closed: Optionally, close the polygon (C{bool}). @kwarg radius: Mean earth radius (C{meter}). @kwarg options: Optional keyword arguments for function L{pygeodesy.equirectangular_}.
@return: A L{NearestOn3Tuple}C{(closest, distance, angle)} of the C{closest} point (L{LatLon}), the L{pygeodesy.equirectangular_} C{distance} between this and the C{closest} point converted to C{meter}, same units as B{C{radius}}. The C{angle} from this to the C{closest} point is in compass C{degrees360}, like function L{pygeodesy.compassAngle}.
@raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}}, see function L{pygeodesy.equirectangular_}.
@raise PointsError: Insufficient number of B{C{points}}.
@raise TypeError: Some B{C{points}} are not C{LatLon}.
@raise ValueError: Invalid B{C{radius}} or B{C{options}}.
@see: Functions L{pygeodesy.compassAngle}, L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}. ''' raise _ValueError(**options)
degrees2m(d, radius=radius), c)
'''Convert this point to C{Karney}-based cartesian (ECEF) coordinates.
@kwarg Cartesian_datum_kwds: Optional L{Cartesian}, B{C{datum}} and other keyword arguments, ignored if C{B{Cartesian} is None}. Use C{B{Cartesian}=...} to override this L{Cartesian} class or specify C{B{Cartesian}=None}.
@return: The cartesian point (L{Cartesian}) or if B{C{Cartesian}} is C{None}, an L{Ecef9Tuple}C{(x, y, z, lat, lon, height, C, M, datum)} with C{C} and C{M} if available.
@raise TypeError: Invalid B{C{Cartesian_datum_kwds}} argument. '''
'''(INTERNAL) Helper for .along-/crossTrackDistanceTo. '''
'''Compute the angles, sides and area of a spherical triangle.
@arg otherB: Second triangle point (C{LatLon}). @arg otherC: Third triangle point (C{LatLon}). @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) or C{None}. @kwarg wrap: Wrap/unroll angular distances (C{bool}).
@return: L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} or if B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)}.
@see: Function L{triangle7} and U{Spherical trigonometry <https://WikiPedia.org/wiki/Spherical_trigonometry>}. '''
area=True, eps=EPS1, radius=R_M, wrap=False): '''Trilaterate three points by area overlap or perimeter intersection of three corresponding circles.
@arg distance1: Distance to this point (C{meter}, same units as B{C{radius}}). @arg point2: Second center point (C{LatLon}). @arg distance2: Distance to point2 (C{meter}, same units as B{C{radius}}). @arg point3: Third center point (C{LatLon}). @arg distance3: Distance to point3 (C{meter}, same units as B{C{radius}}). @kwarg area: If C{True} compute the area overlap, otherwise the perimeter intersection of the circles (C{bool}). @kwarg eps: The required I{minimal overlap} for C{B{area}=True} or the I{intersection margin} for C{B{area}=False} (C{meter}, same units as B{C{radius}}). @kwarg radius: Mean earth radius (C{meter}, conventionally). @kwarg wrap: Wrap/unroll angular distances (C{bool}).
@return: A L{Trilaterate5Tuple}C{(min, minPoint, max, maxPoint, n)} with C{min} and C{max} in C{meter}, same units as B{C{eps}}, the corresponding trilaterated points C{minPoint} and C{maxPoint} as I{spherical} C{LatLon} and C{n}, the number of trilatered points found for the given B{C{eps}}.
If only a single trilaterated point is found, C{min I{is} max}, C{minPoint I{is} maxPoint} and C{n = 1}.
For C{B{area}=True}, C{min} and C{max} are the smallest respectively largest I{radial} overlap found.
For C{B{area}=False}, C{min} and C{max} represent the nearest respectively farthest intersection margin.
If C{B{area}=True} and all 3 circles are concentric, C{n = 0} and C{minPoint} and C{maxPoint} are both the B{C{point#}} with the smallest B{C{distance#}} C{min} and C{max} the largest B{C{distance#}}.
@raise IntersectionError: Trilateration failed for the given B{C{eps}}, insufficient overlap for C{B{area}=True} or no intersection or all (near-)concentric for C{B{area}=False}.
@raise TypeError: Invalid B{C{point2}} or B{C{point3}}.
@raise ValueError: Coincident B{C{points}} or invalid B{C{distance1}}, B{C{distance2}}, B{C{distance3}} or B{C{radius}}. ''' self.others(point2=point2), distance2, self.others(point3=point3), distance3, area=area, radius=radius, eps=eps, wrap=wrap)
'''Calculate the area of a (spherical) polygon (with the points joined by great circle arcs).
@arg points: The polygon points (L{LatLon}[]). @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) or C{None}. @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: Polygon area (C{meter} I{quared}, same units as B{C{radius}} or C{radians} if B{C{radius}} is C{None}).
@raise PointsError: Insufficient number of B{C{points}}.
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@raise ValueError: Invalid B{C{radius}} or semi-circular polygon edge.
@note: The area is based on I{Karney}'s U{'Area of a spherical polygon'<https://MathOverflow.net/questions/97711/ the-area-of-spherical-polygons>}, 3rd Answer.
@see: Funxtions L{pygeodesy.areaOf}, L{sphericalNvector.areaOf}, L{ellipsoidalKarney.areaOf} and L{pygeodesy.excessKarney}.
@example:
>>> b = LatLon(45, 1), LatLon(45, 2), LatLon(46, 2), LatLon(46, 1) >>> areaOf(b) # 8666058750.718977
>>> c = LatLon(0, 0), LatLon(1, 0), LatLon(0, 1) >>> areaOf(c) # 6.18e9 '''
# ispolar: Summation of course deltas around pole is 0° rather than normally ±360° # <https://blog.Element84.com/determining-if-a-spherical-polygon-contains-a-pole.html> # XXX duplicate of function C{points.ispolar} to avoid copying all iterated points
wrap180(z2 - z)) # (z2 - z + 540) % 360 - 180
# see <https://www.EdWilliams.org/intersect.htm> (5) ff # and similar logic in .ellipsoidalBaseDI._intersect3
else: # must be a point
raise _ValueError(_SPACE_(_path_ + _i_, _null_)) # note, in EdWilliams.org/avform.htm W is + and E is -
sa21, _, sa12, _ = sincos2_(b21, b12, a1 - a2, a1 + a2)
sa21 * cb12 * cb21 + sa12 * sb12 * sb21, cos(a1) * cos(a2) * sin(db)) # ll=start
# compute dot product ds . (-b + b1, a - a1)
LatLon=None, **LatLon_kwds): # (INTERNAL) Intersect two (spherical) path, see L{intersection} # above, separated to allow callers to embellish any exceptions
a, b = favg(a1, a2), favg(b1, b2)
# see <https://www.EdWilliams.org/avform.htm#Intersection>
raise IntersectionError(_parallel_) # handle domain error for equivalent longitudes, # see also functions asin_safe and acos_safe at # <https://www.EdWilliams.org/avform.htm#Math> acos1((sa1 - sa2 * cr12) / x2) else: wrapPI(t21 - t23)) # angle 1-2-3) raise IntersectionError(_infinite_) # XXX if sx3 < 0: # XXX raise ValueError(_ambiguous_)
# like .ellipsoidalBaseDI,_intersect3, if this intersection # is "before" the first point, use the antipodal intersection
else: # end point(s) or bearing(s)
raise IntersectionError(_colinear_) # choose intersection similar to sphericalNvector _intdot(d2, a2, b2, a, b, wrap)) > 0:
intersection, LatLon, LatLon_kwds)
LatLon=LatLon, **LatLon_kwds): '''Compute the intersection point of two paths, each defined by two points or a start point and bearing from North.
@arg start1: Start point of the first path (L{LatLon}). @arg end1: End point of the first path (L{LatLon}) or the initial bearing at the first start point (compass C{degrees360}). @arg start2: Start point of the second path (L{LatLon}). @arg end2: End point of the second path (L{LatLon}) or the initial bearing at the second start point (compass C{degrees360}). @kwarg height: Optional height for the intersection point, overriding the mean height (C{meter}). @kwarg wrap: Wrap and unroll longitudes (C{bool}). @kwarg LatLon: Optional class to return the intersection point (L{LatLon}) or C{None}. @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments, ignored if C{B{LatLon} is None}.
@return: The intersection point as a (B{C{LatLon}}) or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon, height)}. An alternate intersection point might be the L{antipode} to the returned result.
@raise IntersectionError: Ambiguous or infinite intersection or colinear, parallel or otherwise non-intersecting paths.
@raise TypeError: A B{C{start}} or B{C{end}} point not L{LatLon}.
@raise ValueError: Invalid B{C{height}} or C{null} path.
@example:
>>> p = LatLon(51.8853, 0.2545) >>> s = LatLon(49.0034, 2.5735) >>> i = intersection(p, 108.547, s, 32.435) # '50.9078°N, 004.5084°E' ''' s1 = _T00.others(start1=start1) s2 = _T00.others(start2=start2) try: return _intersect(s1, end1, s2, end2, height=height, wrap=wrap, LatLon=LatLon, **LatLon_kwds) except (TypeError, ValueError) as x: raise _xError(x, start1=start1, end1=end1, start2=start2, end2=end2)
height=None, wrap=True, LatLon=LatLon, **LatLon_kwds): '''Compute the intersection points of two circles each defined by a center point and a radius.
@arg center1: Center of the first circle (L{LatLon}). @arg rad1: Radius of the first circle (C{meter} or C{radians}, see B{C{radius}}). @arg center2: Center of the second circle (L{LatLon}). @arg rad2: Radius of the second circle (C{meter} or C{radians}, see B{C{radius}}). @kwarg radius: Mean earth radius (C{meter} or C{None} if B{C{rad1}}, B{C{rad2}} and B{C{eps}} are given in C{radians}). @kwarg eps: Required overlap (C{meter} or C{radians}, see B{C{radius}}). @kwarg height: Optional height for the intersection points (C{meter}, conventionally) or C{None} for the I{"radical height"} at the I{radical line} between both centers. @kwarg wrap: Wrap and unroll longitudes (C{bool}). @kwarg LatLon: Optional class to return the intersection points (L{LatLon}) or C{None}. @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments, ignored if C{B{LatLon} is None}.
@return: 2-Tuple of the intersection points, each a B{C{LatLon}} instance or if C{B{LatLon} is None} a L{LatLon3Tuple}C{(lat, lon, height)}. For abutting circles, both intersection points are the same instance, aka I{radical center}.
@raise IntersectionError: Concentric, antipodal, invalid or non-intersecting circles.
@raise TypeError: If B{C{center1}} or B{C{center2}} not L{LatLon}.
@raise ValueError: Invalid B{C{rad1}}, B{C{rad2}}, B{C{radius}}, B{C{eps}} or B{C{height}}.
@note: Courtesy of U{Samuel Čavoj<https://GitHub.com/mrJean1/PyGeodesy/issues/41>}.
@see: This U{Answer<https://StackOverflow.com/questions/53324667/ find-intersection-coordinates-of-two-circles-on-earth/53331953>}. ''' height=height, wrap=wrap, LatLon=LatLon, **LatLon_kwds) except (TypeError, ValueError) as x: raise _xError(x, center1=center1, rad1=rad1, center2=center2, rad2=rad2)
height=None, too_d=None, wrap=True, LatLon=LatLon, **LatLon_kwds): # (INTERNAL) Intersect two spherical circles, see L{intersections2} # above, separated to allow callers to embellish any exceptions
intersections2, LatLon, LatLon_kwds)
raise IntersectionError(_near_(_concentric_))
eps, radius=radius) if eps else _0_0) raise _ValueError(eps=r)
raise _ValueError(_invalid_)
elif x < r: # PYCHOK no cover t = (d * radius) if too_d is None else too_d raise IntersectionError(_too_(_Fmt.distant(t)))
else:
r = _dest3(b, h) else:
def isPoleEnclosedBy(points, wrap=False): # PYCHOK no cover '''DEPRECATED, use function L{pygeodesy.ispolar}. ''' return ispolar(points, wrap=wrap)
'''(INTERNAL) Helper for L{intersection}, L{intersections2} and L{meanOf}. ''' r = LatLon3Tuple(lat, lon, height, name=n) else:
'''Compute the geographic mean of several points.
@arg points: Points to be averaged (L{LatLon}[]). @kwarg height: Optional height at mean point, overriding the mean height (C{meter}). @kwarg LatLon: Optional class to return the mean point (L{LatLon}) or C{None}. @kwarg LatLon_kwds: Optional, additional B{C{LatLon}} keyword arguments, ignored if C{B{LatLon} is None}.
@return: The geographic mean and height (B{C{LatLon}}) or a L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}} is C{None}.
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@raise ValueError: No B{C{points}} or invalid B{C{height}}. ''' # geographic mean
else: # PYCHOK no cover h = Height(height)
def nearestOn2(point, points, **closed_radius_LatLon_options): # PYCHOK no cover '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
@return: ... 2-tuple C{(closest, distance)} of the C{closest} point (L{LatLon}) on the polygon and the C{distance} between the C{closest} and the given B{C{point}}. The C{closest} is a B{C{LatLon}} or a L{LatLon2Tuple}C{(lat, lon)} if B{C{LatLon}} is C{None} ... ''' ll, d, _ = nearestOn3(point, points, **closed_radius_LatLon_options) # PYCHOK 3-tuple if _xkwds_get(closed_radius_LatLon_options, LatLon=LatLon) is None: ll = LatLon2Tuple(ll.lat, ll.lon) return ll, d
LatLon=LatLon, **options): '''Locate the point on a path or polygon closest to a reference point.
Distances are I{approximated} using function L{pygeodesy.equirectangular_}, subject to the supplied B{C{options}}.
@arg point: The reference point (L{LatLon}). @arg points: The path or polygon points (L{LatLon}[]). @kwarg closed: Optionally, close the polygon (C{bool}). @kwarg radius: Mean earth radius (C{meter}). @kwarg LatLon: Optional class to return the closest point (L{LatLon}) or C{None}. @kwarg options: Optional keyword arguments for function L{pygeodesy.equirectangular_}.
@return: A L{NearestOn3Tuple}C{(closest, distance, angle)} with the C{closest} point as B{C{LatLon}} or L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}} is C{None}. The C{distance} is the L{pygeodesy.equirectangular_} distance between the C{closest} and the given B{C{point}} converted to C{meter}, same units as B{C{radius}}. The C{angle} from the given B{C{point}} to the C{closest} is in compass C{degrees360}, like function L{pygeodesy.compassAngle}. The C{height} is the (interpolated) height at the C{closest} point.
@raise LimitError: Lat- and/or longitudinal delta exceeds the B{C{limit}}, see function L{pygeodesy.equirectangular_}.
@raise PointsError: Insufficient number of B{C{points}}.
@raise TypeError: Some B{C{points}} are not C{LatLon}.
@raise ValueError: Invalid B{C{radius}}.
@see: Functions L{pygeodesy.equirectangular_} and L{pygeodesy.nearestOn5}. ''' LatLon=None, **options) LatLon(lat, lon, height=h, name=n)
'''Compute the perimeter of a (spherical) polygon (with great circle arcs joining the points).
@arg points: The polygon points (L{LatLon}[]). @kwarg closed: Optionally, close the polygon (C{bool}). @kwarg radius: Mean earth radius (C{meter}) or C{None}. @kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: Polygon perimeter (C{meter}, same units as B{C{radius}} or C{radians} if B{C{radius}} is C{None}).
@raise PointsError: Insufficient number of B{C{points}}.
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@raise ValueError: Invalid B{C{radius}}.
@note: This perimeter is based on the L{pygeodesy.haversine} formula.
@see: Functions L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf} and L{ellipsoidalKarney.perimeterOf}. '''
excess=excessAbc, wrap=False): '''Compute the angles, sides, and area of a (spherical) triangle.
@arg latA: First corner latitude (C{degrees}). @arg lonA: First corner longitude (C{degrees}). @arg latB: Second corner latitude (C{degrees}). @arg lonB: Second corner longitude (C{degrees}). @arg latC: Third corner latitude (C{degrees}). @arg lonC: Third corner longitude (C{degrees}). @kwarg radius: Mean earth radius, ellipsoid or datum (C{meter}, L{Ellipsoid}, L{Ellipsoid2}, L{Datum} or L{a_f2Tuple}) or C{None}. @kwarg excess: I{Spherical excess} callable (L{excessAbc}, L{excessGirard} or L{excessLHuilier}). @kwarg wrap: Wrap and L{pygeodesy.unroll180} longitudes (C{bool}).
@return: A L{Triangle7Tuple}C{(A, a, B, b, C, c, area)} with spherical angles C{A}, C{B} and C{C}, angular sides C{a}, C{b} and C{c} all in C{degrees} and C{area} in I{square} C{meter} or units of B{C{radius}} I{squared} or if B{C{radius}} is C{None}, a L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} all in C{radians} with I{spherical excess} C{E} as the C{unit area} in C{radians}. ''' t = triangle8_(Phi_(latA=latA), Lam_(lonA=lonA), Phi_(latB=latB), Lam_(lonB=lonB), Phi_(latC=latC), Lam_(lonC=lonC), excess=excess, wrap=wrap) return _t7Tuple(t, radius)
wrap=False): '''Compute the angles, sides, I{spherical deficit} and I{spherical excess} of a (spherical) triangle.
@arg phiA: First corner latitude (C{radians}). @arg lamA: First corner longitude (C{radians}). @arg phiB: Second corner latitude (C{radians}). @arg lamB: Second corner longitude (C{radians}). @arg phiC: Third corner latitude (C{radians}). @arg lamC: Third corner longitude (C{radians}). @kwarg excess: I{Spherical excess} callable (L{excessAbc}, L{excessGirard} or L{excessLHuilier}). @kwarg wrap: Wrap and L{pygeodesy.unrollPI} longitudes (C{bool}).
@return: A L{Triangle8Tuple}C{(A, a, B, b, C, c, D, E)} with spherical angles C{A}, C{B} and C{C}, angular sides C{a}, C{b} and C{c}, I{spherical deficit} C{D} and I{spherical excess} C{E}, all in C{radians}. '''
# notation: side C{a} is oposite to corner C{A}, etc.
excessLHuilier(a, b, c) if excess in (excessLHuilier, False) else excessAbc(*max((A, b, c), (B, c, a), (C, a, b))))
'''(INTERNAL) Convert a L{Triangle8Tuple} to L{Triangle7Tuple}. ''' _ellipsoidal_datum(radius).ellipsoid.Rmean B, (r * t.b), C, (r * t.c), t.E * r**2)
areaOf, # functions intersection, intersections2, ispolar, isPoleEnclosedBy, # DEPRECATED, use ispolar meanOf, nearestOn2, nearestOn3, perimeterOf, sumOf, # == vector3d.sumOf triangle7, triangle8_)
# **) MIT License # # Copyright (C) 2016-2022 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE. |