Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def TestEdge(self, azi, s, reverse, sign):
"""Return the results for a tentative additional edge."""
if self._num == 0: # we don't have a starting point!
return 0, Math.nan, Math.nan
num = self._num + 1
perimeter = self._perimetersum.Sum() + s
if self._polyline:
return num, perimeter, Math.nan
tempsum = self._areasum.Sum()
crossings = self._crossings
_, lat, lon, _, _, _, _, _, S12 = self._earth.GenDirect(
self._lat1, self._lon1, azi, False, s, self._mask)
tempsum += S12
crossings += PolygonArea.transit(self._lon1, lon)
_, s12, _, _, _, _, _, S12 = self._earth.GenInverse(
lat, lon, self._lat0, self._lon0, self._mask)
perimeter += s12
tempsum += S12
crossings += PolygonArea.transit(lon, self._lon0)
def __init__(self, geod, lat1, lon1, azi1,
caps = GeodesicCapability.STANDARD |
GeodesicCapability.DISTANCE_IN,
salp1 = Math.nan, calp1 = Math.nan):
"""Construct a GeodesicLine object
:param geod: a :class:`~geographiclib.geodesic.Geodesic` object
:param lat1: latitude of the first point in degrees
:param lon1: longitude of the first point in degrees
:param azi1: azimuth at the first point in degrees
:param caps: the :ref:`capabilities `
This creates an object allowing points along a geodesic starting at
(*lat1*, *lon1*), with azimuth *azi1* to be found. The default
value of *caps* is STANDARD | DISTANCE_IN. The optional parameters
*salp1* and *calp1* should not be supplied; they are part of the
private interface.
"""
from geographiclib.geodesic import Geodesic
self.earth = earth
"""The geodesic object (readonly)"""
self.polyline = polyline
"""Is this a polyline? (readonly)"""
self.area0 = 4 * math.pi * earth._c2
"""The total area of the ellipsoid in meter^2 (readonly)"""
self._mask = (Geodesic.LATITUDE | Geodesic.LONGITUDE |
Geodesic.DISTANCE |
(Geodesic.EMPTY if self.polyline else
Geodesic.AREA | Geodesic.LONG_UNROLL))
if not self.polyline: self._areasum = Accumulator()
self._perimetersum = Accumulator()
self.num = 0
"""The current number of points in the polygon (readonly)"""
self.lat1 = Math.nan
"""The current latitude in degrees (readonly)"""
self.lon1 = Math.nan
"""The current longitude in degrees (readonly)"""
self.Clear()
def _Lengths(self, eps, sig12,
ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, outmask,
# Scratch areas of the right size
C1a, C2a):
"""Private: return a bunch of lengths"""
# Return s12b, m12b, m0, M12, M21, where
# m12b = (reduced length)/_b; s12b = distance/_b,
# m0 = coefficient of secular term in expression for reduced length.
outmask &= Geodesic.OUT_MASK
# outmask & DISTANCE: set s12b
# outmask & REDUCEDLENGTH: set m12b & m0
# outmask & GEODESICSCALE: set M12 & M21
s12b = m12b = m0 = M12 = M21 = Math.nan
if outmask & (Geodesic.DISTANCE | Geodesic.REDUCEDLENGTH |
Geodesic.GEODESICSCALE):
A1 = Geodesic._A1m1f(eps)
Geodesic._C1f(eps, C1a)
if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
A2 = Geodesic._A2m1f(eps)
Geodesic._C2f(eps, C2a)
m0x = A1 - A2
A2 = 1 + A2
A1 = 1 + A1
if outmask & Geodesic.DISTANCE:
B1 = (Geodesic._SinCosSeries(True, ssig2, csig2, C1a) -
Geodesic._SinCosSeries(True, ssig1, csig1, C1a))
# Missing a factor of _b
s12b = A1 * (sig12 + B1)
if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
def _InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
lam12, slam12, clam12,
# Scratch areas of the right size
C1a, C2a):
"""Private: Find a starting value for Newton's method."""
# Return a starting point for Newton's method in salp1 and calp1 (function
# value is -1). If Newton's method doesn't need to be used, return also
# salp2 and calp2 and function value is sig12.
sig12 = -1; salp2 = calp2 = dnm = Math.nan # Return values
# bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0]
sbet12 = sbet2 * cbet1 - cbet2 * sbet1
cbet12 = cbet2 * cbet1 + sbet2 * sbet1
# Volatile declaration needed to fix inverse cases
# 88.202499451857 0 -88.202499451857 179.981022032992859592
# 89.262080389218 0 -89.262080389218 179.992207982775375662
# 89.333123580033 0 -89.333123580032997687 179.99295812360148422
# which otherwise fail with g++ 4.4.4 x86 -O3
sbet12a = sbet2 * cbet1
sbet12a += cbet2 * sbet1
shortline = cbet12 >= 0 and sbet12 < 0.5 and cbet2 * lam12 < 0.5
if shortline:
sbetm2 = Math.sq(sbet1 + sbet2)
# sin((bet1+bet2)/2)^2
# = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
def TestPoint(self, lat, lon, reverse = False, sign = True):
"""Compute the properties for a tentative additional vertex
:param lat: the latitude of the point in degrees
:param lon: the longitude of the point in degrees
:param reverse: if true then clockwise (instead of
counter-clockwise) traversal counts as a positive area
:param sign: if true then return a signed result for the area if the
polygon is traversed in the "wrong" direction instead of returning
the area for the rest of the earth
:return: a tuple of number, perimeter (meters), area (meters^2)
"""
if self.polyline: area = Math.nan
if self.num == 0:
perimeter = 0.0
if not self.polyline: area = 0.0
return 1, perimeter, area
perimeter = self._perimetersum.Sum()
tempsum = 0.0 if self.polyline else self._areasum.Sum()
crossings = self._crossings; num = self.num + 1
for i in ([0] if self.polyline else [0, 1]):
_, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
self.lat1 if i == 0 else lat, self.lon1 if i == 0 else lon,
self._lat0 if i != 0 else lat, self._lon0 if i != 0 else lon,
self._mask)
perimeter += s12
if not self.polyline:
tempsum += S12
def _GenInverse(self, lat1, lon1, lat2, lon2, outmask):
"""Private: General version of the inverse problem"""
a12 = s12 = m12 = M12 = M21 = S12 = Math.nan # return vals
outmask &= Geodesic.OUT_MASK
# Compute longitude difference (AngDiff does this carefully). Result is
# in [-180, 180] but -180 is only for west-going geodesics. 180 is for
# east-going and meridional geodesics.
lon12, lon12s = Math.AngDiff(lon1, lon2)
# Make longitude difference positive.
lonsign = 1 if lon12 >= 0 else -1
# If very close to being on the same half-meridian, then make it so.
lon12 = lonsign * Math.AngRound(lon12)
lon12s = Math.AngRound((180 - lon12) - lonsign * lon12s)
lam12 = math.radians(lon12)
if lon12 > 90:
slam12, clam12 = Math.sincosd(lon12s); clam12 = -clam12
else:
slam12, clam12 = Math.sincosd(lon12)
def _Lengths(self, eps, sig12,
ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, outmask,
# Scratch areas of the right size
C1a, C2a):
"""Private: return a bunch of lengths"""
# Return s12b, m12b, m0, M12, M21, where
# m12b = (reduced length)/_b; s12b = distance/_b,
# m0 = coefficient of secular term in expression for reduced length.
outmask &= Geodesic.OUT_MASK
# outmask & DISTANCE: set s12b
# outmask & REDUCEDLENGTH: set m12b & m0
# outmask & GEODESICSCALE: set M12 & M21
s12b = m12b = m0 = M12 = M21 = Math.nan
if outmask & (Geodesic.DISTANCE | Geodesic.REDUCEDLENGTH |
Geodesic.GEODESICSCALE):
A1 = Geodesic._A1m1f(eps)
Geodesic._C1f(eps, C1a)
if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
A2 = Geodesic._A2m1f(eps)
Geodesic._C2f(eps, C2a)
m0x = A1 - A2
A2 = 1 + A2
A1 = 1 + A1
if outmask & Geodesic.DISTANCE:
B1 = (Geodesic._SinCosSeries(True, ssig2, csig2, C1a) -
Geodesic._SinCosSeries(True, ssig1, csig1, C1a))
# Missing a factor of _b
s12b = A1 * (sig12 + B1)
if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
def Clear(self):
"""Reset to empty polygon."""
self._num = 0
self._crossings = 0
if not self._polyline: self._areasum.Set(0)
self._perimetersum.Set(0)
self._lat0 = self._lon0 = self._lat1 = self._lon1 = Math.nan
def _GenInverse(self, lat1, lon1, lat2, lon2, outmask):
"""Private: General version of the inverse problem"""
a12 = s12 = m12 = M12 = M21 = S12 = Math.nan # return vals
outmask &= Geodesic.OUT_MASK
# Compute longitude difference (AngDiff does this carefully). Result is
# in [-180, 180] but -180 is only for west-going geodesics. 180 is for
# east-going and meridional geodesics.
lon12, lon12s = Math.AngDiff(lon1, lon2)
# Make longitude difference positive.
lonsign = 1 if lon12 >= 0 else -1
# If very close to being on the same half-meridian, then make it so.
lon12 = lonsign * Math.AngRound(lon12)
lon12s = Math.AngRound((180 - lon12) - lonsign * lon12s)
lam12 = math.radians(lon12)
if lon12 > 90:
slam12, clam12 = Math.sincosd(lon12s); clam12 = -clam12
else:
slam12, clam12 = Math.sincosd(lon12)