Coverage for pygeodesy/sphericalTrigonometry.py : 91%

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 -*-
geocentric (ECEF) L{Cartesian} and functionsL{areaOf}, L{intersection}, L{isPoleEnclosedBy}, L{meanOf}, L{nearestOn2} and L{perimeterOf}, I{all spherical}.
Pure Python implementation of geodetic (lat-/longitude) methods using spherical trigonometry, transcribed 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>}.
@newfield example: Example, Examples '''
isscalar, map1 degrees2m, iterNumpy2, radiansPI2, \ sincos2, tan_2, unrollPI, wrapPI
radians, sin
# all public contants, classes and functions 'Cartesian', 'LatLon', # classes 'areaOf', # functions 'intersection', 'ispolar', 'isPoleEnclosedBy', # DEPRECATED, use ispolar 'meanOf', 'nearestOn2', 'nearestOn3', 'perimeterOf', 'sumOf') # == vector3d.sumOf
'''Extended to convert geocentric, L{Cartesian} points to spherical, geodetic L{LatLon}. '''
'''Convert this cartesian point to an C{Nvector}-based geodetic point.
@keyword kwds: Optional, additional B{C{LatLon}} keyword arguments, ignored if C{B{LatLon}=None}. For example, use C{LatLon=...} to override the L{LatLon} (sub-)class or specify
@return: The B{C{LatLon}} point (L{LatLon}) or if C{B{LatLon}=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}}, B{C{datum}} or B{C{kwds}}. ''' kwds = _2kwds(kwds, LatLon=LatLon, datum=self.datum) return CartesianSphericalBase.toLatLon(self, **kwds)
'''New point on spherical model earth model.
@example:
>>> p = LatLon(52.205, 0.119) # height=0 '''
'''(INTERNAL) Helper for along-/crossTrackDistanceTo. '''
raise ValueError('%s invalid: %r' % ('radius', radius))
'''Compute the (signed) distance 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.
@param start: Start point of great circle path (L{LatLon}). @param end: End point of great circle path (L{LatLon}). @keyword radius: Optional, mean earth radius (C{meter}). @keyword wrap: Wrap and unroll longitudes (C{bool}).
@return: Distance along the great circle path (C{meter}, same units as B{C{radius}}), 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: The B{C{start}} or B{C{end}} point is not L{LatLon}.
@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 ''' else: return 0.0
'''DEPRECATED, use method C{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.
@param other: The other point defining great circle (L{LatLon}). @param lat: Latitude at the crossing (C{degrees}). @keyword 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)
return None # great circle doesn't reach latitude
'''Compute the (signed) distance from this point to the great circle defined by a start and an end point.
@param start: Start point of great circle path (L{LatLon}). @param end: End point of great circle path (L{LatLon}). @keyword radius: Optional, mean earth radius (C{meter}). @keyword 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).
@raise TypeError: The B{C{start}} or B{C{end}} point is not L{LatLon}.
@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.
@param distance: Distance travelled (C{meter}, same units as B{C{radius}}). @param bearing: Bearing from this point (compass C{degrees360}). @keyword radius: Optional, mean earth radius (C{meter}). @keyword height: Optional height at destination (C{meter}, same units a B{C{radius}}).
@return: Destination point (L{LatLon}).
@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 distance from this to an other point.
@param other: The other point (L{LatLon}). @keyword radius: Optional, mean earth radius (C{meter}). @keyword 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}}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@example:
>>> p1 = LatLon(52.205, 0.119) >>> p2 = LatLon(48.857, 2.351); >>> d = p1.distanceTo(p2) # 404300 '''
'''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.
@param bearing: Bearing from this point (compass C{degrees360}).
@return: Vector representing great circle (L{Vector3d}).
@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.
@param other: The other point (spherical L{LatLon}). @keyword wrap: Wrap and unroll longitudes (C{bool}). @keyword raiser: Optionally, raise L{CrossError} (C{bool}), use B{C{raiser}}=C{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 B{C{raiser}} is C{True} and L{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('%s %s: %r' % ('coincident', 'points', other))
'''Locate the point at given fraction between this and an other point.
@param other: The other point (L{LatLon}). @param fraction: Fraction between both points (float, 0.0 = this point, 1.0 = the other point). @keyword height: Optional height, overriding the fractional height (C{meter}). @keyword wrap: Wrap and unroll longitudes (C{bool}).
@return: Intermediate point (L{LatLon}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@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}. '''
sb1, cb1, sb2, cb2 = sincos2(a1, a2, b1, b2)
else: # points too close a = favg(a1, a2, f=fraction) b = favg(b1, b2, f=fraction)
else: h = height
height=None, wrap=False): '''Locate the intersection point of two paths both defined by two points or a start point and bearing from North.
@param end1: End point of the first path (L{LatLon}) or the initial bearing at this point (compass C{degrees360}). @param start2: Start point of the second path (L{LatLon}). @param end2: End point of the second path (L{LatLon}) or the initial bearing at the second point (compass C{degrees360}). @keyword height: Optional height for intersection point, overriding the mean height (C{meter}). @keyword 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 TypeError: The B{C{start2}}, B{C{end1}} or B{C{end2}} is not L{LatLon}.
@raise ValueError: Intersection is ambiguous or infinite or the paths are parallel, coincident or null.
@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' ''' height=height, wrap=wrap, LatLon=self.classof)
'''Check whether a (convex) polygon encloses this point.
@param points: The polygon points (L{LatLon}[]).
@return: C{True} if the polygon encloses this point, C{False} otherwise.
@raise ValueError: Insufficient number of B{C{points}} or non-convex polygon.
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@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 '''
return False # outside
raise ValueError('non-convex: %r...' % (points[:2],))
else: # get great-circle vector for each edge
# check whether this point on same side of all # polygon edges (to the left or right depending # on anti-/clockwise polygon direction) return False # outside
# check for convex polygon (otherwise # the test above is not reliable) # angle between gc vectors, signed by direction of n0 raise ValueError('non-convex: %r...' % (points[:2],))
'''DEPRECATED, use method C{isenclosedBy}. ''' return self.isenclosedBy(points)
'''Find the midpoint between this and an other point.
@param other: The other point (L{LatLon}). @keyword height: Optional height for midpoint, overriding the mean height (C{meter}). @keyword wrap: Wrap and unroll longitudes (C{bool}).
@return: Midpoint (L{LatLon}).
@raise TypeError: The B{C{other}} point is not L{LatLon}.
@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 and this point.
Distances are approximated by function L{equirectangular_}, subject to the supplied B{C{options}}.
@param point1: Start point (L{LatLon}). @param point2: End point (L{LatLon}). @keyword radius: Optional, mean earth radius (C{meter}). @keyword options: Optional keyword arguments for function L{equirectangular_}.
@return: Closest point on the arc (L{LatLon}).
@raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}}, see function L{equirectangular_}.
@raise TypeError: If B{C{point1}} or B{C{point2}} is not L{LatLon}.
@see: Functions L{equirectangular_} and L{nearestOn5} and method L{sphericalTrigonometry.LatLon.nearestOn3}. ''' **options)[0]
'''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 ... ''' **options)[:2]
'''Locate the point on a polygon closest to this point.
Distances are approximated by function L{equirectangular_}, subject to the supplied B{C{options}}.
@param points: The polygon points (L{LatLon}[]). @keyword closed: Optionally, close the polygon (C{bool}). @keyword radius: Optional, mean earth radius (C{meter}). @keyword options: Optional keyword arguments for function L{equirectangular_}.
@return: A L{NearestOn3Tuple}C{(closest, distance, angle)} where C{distance} is the L{equirectangular_} distance between this and the C{closest} point in 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{compassAngle}.
@raise LimitError: Lat- and/or longitudinal delta exceeds B{C{limit}}, see function L{equirectangular_}.
@raise TypeError: Some B{C{points}} are not C{LatLon}.
@raise ValueError: Insufficient number of B{C{points}}.
@see: Functions L{compassAngle}, L{equirectangular_} and L{nearestOn5}. ''' degrees2m(d, radius=radius), c)
'''Convert this point to C{Karney}-based cartesian (ECEF) coordinates.
@keyword kwds: Optional, additional B{C{Cartesian}} keyword arguments, ignored if C{B{Cartesian}=None}. For example, use C{Cartesian=...} to override the L{Cartesian} (sub-)class or specify
@return: The B{C{Cartesian}} point (L{Cartesian}) or if C{B{Cartesian}=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}}, B{C{datum}} or B{C{kwds}}. ''' kwds = _2kwds(kwds, Cartesian=Cartesian, datum=self.datum) return LatLonSphericalBase.toCartesian(self, **kwds)
'''(INTERNAL) Computes destination lat- and longitude.
@param a: Latitude (C{radians}). @param b: Longitude (C{radians}). @param r: Angular distance (C{radians}). @param t: Bearing (compass C{radians}).
@return: 2-Tuple (lat, lon) of (radians, radians). ''' # see <https://www.EdWilliams.org/avform.htm#LL>
# note, in EdWilliams.org/avform.htm W is + and E is -
'''(INTERNAL) Computes destination lat- and longitude.
@param a: Latitude (C{radians}). @param b: Longitude (C{radians}). @param r: Angular distance (C{radians}). @param t: Bearing (compass C{radians}).
@return: 2-Tuple (lat, lon) of (C{degrees90}, C{degrees180}). '''
'''Calculate the area of a (spherical) polygon (with great circle arcs joining the points).
@param points: The polygon points (L{LatLon}[]). @keyword radius: Optional, mean earth radius (C{meter}). @keyword wrap: Wrap and unroll longitudes (C{bool}).
@return: Polygon area (C{meter}, same units as B{C{radius}}, squared).
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@raise ValueError: Insufficient number of B{C{points}}.
@note: The area is based on Karney's U{'Area of a spherical polygon' <https://OSGeo-org.1560.x6.nabble.com/ Area-of-a-spherical-polygon-td3841625.html>}.
@see: L{pygeodesy.areaOf}, L{sphericalNvector.areaOf} and L{ellipsoidalKarney.areaOf}.
@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 '''
# Area method due to Karney: for each edge of the polygon, # # tan(Δλ/2) · (tan(φ1/2) + tan(φ2/2)) # tan(E/2) = ------------------------------------ # 1 + tan(φ1/2) · tan(φ2/2) # # where E is the spherical excess of the trapezium obtained by # extending the edge to the equator-circle vector for each edge
# see <https://www.EdWilliams.org/intersect.htm> (5) ff
else: # must be a point
raise ValueError('intersection %s%s null: %r' % ('path', n, (start, end)))
# 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)
# difference between the bearing to (a, b) and the given # bearing is negative if both are in opposite directions
# compute dot product d . (-b + b1, a - a1)
height=None, wrap=False, LatLon=LatLon): '''Compute the intersection point of two paths both defined by two points or a start point and bearing from North.
@param start1: Start point of the first path (L{LatLon}). @param end1: End point ofthe first path (L{LatLon}) or the initial bearing at the first start point (compass C{degrees360}). @param start2: Start point of the second path (L{LatLon}). @param end2: End point of the second path (L{LatLon}) or the initial bearing at the second start point (compass C{degrees360}). @keyword height: Optional height for the intersection point, overriding the mean height (C{meter}). @keyword wrap: Wrap and unroll longitudes (C{bool}). @keyword LatLon: Optional (sub-)class to return the intersection point (L{LatLon}) or C{None}.
@return: The intersection point (B{C{LatLon}}) or a L{LatLon3Tuple}C{(lat, lon, height)} if B{C{LatLon}} is C{None}. An alternate intersection point might be the L{antipode} to the returned result.
@raise TypeError: A B{C{start}} or B{C{end}} point not L{LatLon}.
@raise ValueError: Intersection is ambiguous or infinite or the paths are parallel, coincident or null.
@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' '''
a, b = map1(degrees, favg(a1, a2), favg(b1, b2))
# see <https://www.EdWilliams.org/avform.htm#Intersection>
raise ValueError('intersection %s: %r vs %r' % ('parallel', (start1, end1), (start2, end2)))
# handle domain error for equivalent longitudes, # see also functions asin_safe and acos_safe at # <https://www.EdWilliams.org/avform.htm#Math> (sa1 - sa2 * cr12) / x2) else:
t21 - t23) # angle 1-2-3 raise ValueError('intersection %s: %r vs %r' % ('infinite', (start1, end1), (start2, end2))) # if sx3 < 0: # raise ValueError('intersection %s: %r vs %r' % ('ambiguous', # (start1, end1), (start2, end2)))
# choose antipode for opposing bearings _xb(a2, b2, end2, a, b, wrap) < 0:
else: # end point(s) or bearing(s) raise ValueError('intersection %s: %r vs %r' % ('colinear', (start1, end1), (start2, end2))) # choose intersection similar to sphericalNvector
LatLon(a, b, height=h)
'''DEPRECATED, use function L{ispolar}. '''
'''Compute the geographic mean of several points.
@param points: Points to be averaged (L{LatLon}[]). @keyword height: Optional height at mean point, overriding the mean height (C{meter}). @keyword LatLon: Optional (sub-)class to return the mean point (L{LatLon}) or C{None}.
@return: Point at 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}}. ''' # geographic mean n, points = _Trll.points2(points, closed=False)
m = sumOf(points[i].toVector3d() for i in range(n)) a, b = m.to2ll()
if height is None: h = fmean(points[i].height for i in range(n)) else: h = height r = LatLon3Tuple(a, b, h) if LatLon is None else \ LatLon(a, b, height=h) return _xnamed(r, meanOf.__name__)
LatLon=LatLon, **options): '''DEPRECATED, use function L{sphericalTrigonometry.nearestOn3}.
@return: ... C{closest} as B{C{LatLon}} or a 2-tuple C{(lat, lon)} without the height if B{C{LatLon}} is C{None} ... '''
LatLon=LatLon, **options): '''Locate the point on a polygon closest to an other, reference point.
Distances are approximated by function L{equirectangular_}, subject to the supplied B{C{options}}.
@param point: The other, reference point (L{LatLon}). @param points: The polygon points (L{LatLon}[]). @keyword closed: Optionally, close the polygon (C{bool}). @keyword radius: Optional, mean earth radius (C{meter}). @keyword LatLon: Optional (sub-)class to return the closest point (L{LatLon}) or C{None}. @keyword options: Optional keyword arguments for function L{equirectangular_}.
@return: A L{NearestOn3Tuple}C{(closest, distance, angle)}. The C{distance} is the L{equirectangular_} distance between the C{closest} and reference B{C{point}} in C{meter}, same units as B{C{radius}}. The C{angle} from the reference B{C{point}} to the C{closest} is in compass C{degrees360}, like function L{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{equirectangular_}.
@raise TypeError: Some I{points} are not C{LatLon}.
@raise ValueError: Insufficient number of B{C{points}}.
@see: Functions L{equirectangular_} and L{nearestOn5}. ''' a, b, d, c, h = _nearestOn5(point, points, closed=closed, LatLon=None, **options) r = LatLon3Tuple(a, b, h) if LatLon is None else \ LatLon(a, b, height=h) r = NearestOn3Tuple(r, degrees2m(d, radius=radius), c) return _xnamed(r, nearestOn3.__name__)
'''Compute the perimeter of a (spherical) polygon (with great circle arcs joining the points).
@param points: The polygon points (L{LatLon}[]). @keyword closed: Optionally, close the polygon (C{bool}). @keyword radius: Optional, mean earth radius (C{meter}). @keyword wrap: Wrap and unroll longitudes (C{bool}).
@return: Polygon perimeter (C{meter}, same units as B{C{radius}}).
@raise TypeError: Some B{C{points}} are not L{LatLon}.
@raise ValueError: Insufficient number of B{C{points}}.
@note: This perimeter is based on the L{haversine} formula.
@see: L{pygeodesy.perimeterOf}, L{sphericalNvector.perimeterOf} and L{ellipsoidalKarney.perimeterOf}. '''
# **) MIT License # # Copyright (C) 2016-2020 -- 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. |