/usr/share/pyshared/geographiclib/polygonarea.py is in python-geographiclib 1.21-1ubuntu1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 | # polygonarea.py
#
# This is a rather literal translation of the GeographicLib::PolygonArea
# class to python. See the documentation for the C++ class for more
# information at
#
# http://geographiclib.sourceforge.net/html/annotated.html
#
# Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed
# under the MIT/X11 License. For more information, see
# http://geographiclib.sourceforge.net/
#
# $Id: 4a2cf9d0cb7eeb7c9c3a2248f44e84a04e4e0b35 $
######################################################################
import math
from geographiclib.geomath import Math
from geographiclib.accumulator import Accumulator
class PolygonArea(object):
"""Area of a geodesic polygon"""
def transit(lon1, lon2):
# Return 1 or -1 if crossing prime meridian in east or west direction.
# Otherwise return zero.
from geographiclib.geodesic import Geodesic
lon1 = Geodesic.AngNormalize(lon1)
lon2 = Geodesic.AngNormalize(lon2)
# treat lon12 = -180 as an eastward geodesic, so convert to 180.
lon12 = -Geodesic.AngNormalize(lon1 - lon2) # In (-180, 180]
cross = (1 if lon1 < 0 and lon2 >= 0 and lon12 > 0
else (-1 if lon2 < 0 and lon1 >= 0 and lon12 < 0 else 0))
return cross
transit = staticmethod(transit)
def __init__(self, earth, polyline = False):
from geographiclib.geodesic import Geodesic
self._earth = earth
self._area0 = 4 * math.pi * earth._c2
self._polyline = polyline
self._mask = Geodesic.DISTANCE | (0 if self._polyline else Geodesic.AREA)
if not self._polyline: self._areasum = Accumulator()
self._perimetersum = Accumulator()
self.Clear()
def Clear(self):
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 = 0
def AddPoint(self, lat, lon):
if self._num == 0:
self._lat0 = self._lat1 = lat
self._lon0 = self._lon1 = lon
else:
t, s12, t, t, t, t, t, S12 = self._earth.GenInverse(
self._lat1, self._lon1, lat, lon, self._mask)
self._perimetersum.Add(s12)
if not self._polyline:
self._areasum.Add(S12)
self._crossings += PolygonArea.transit(self._lon1, lon)
self._lat1 = lat
self._lon1 = lon
self._num += 1
# return number, perimeter, area
def Compute(self, reverse, sign):
if self._polyline: area = Math.nan
if self._num < 2:
perimeter = 0
if not self._polyline: area = 0
return self._num, perimeter, area
if self._polyline:
perimeter = self._perimetersum.Sum()
return self._num, perimeter, area
t, s12, t, t, t, t, t, S12 = self._earth.GenInverse(
self._lat1, self._lon1, self._lat0, self._lon0, self._mask)
perimeter = self._perimetersum.Sum(s12)
tempsum = Accumulator(self._areasum)
tempsum.Add(S12)
crossings = self._crossings + PolygonArea.transit(self._lon1, self._lon0)
if crossings & 1:
tempsum.Add( (1 if tempsum < 0 else -1) * self._area0/2 )
# area is with the clockwise sense. If !reverse convert to
# counter-clockwise convention.
if not reverse: tempsum.Negate()
# If sign put area in (-area0/2, area0/2], else put area in [0, area0)
if sign:
if tempsum.Sum() > self._area0/2:
tempsum.Add( -self._area0 )
elif tempsum.Sum() <= -self._area0/2:
tempsum.Add( +self._area0 )
else:
if tempsum.Sum() >= self._area0:
tempsum.Add( -self._area0 )
elif tempsum.Sum() < 0:
tempsum.Add( +self._area0 )
area = 0 + tempsum.Sum()
return self._num, perimeter, area
# return number, perimeter, area
def TestCompute(self, lat, lon, reverse, sign):
if self._polyline: area = Math.nan
if self._num == 0:
perimeter = 0
if not self._polyline: area = 0
return 1, perimeter, area
perimeter = self._perimetersum.Sum()
tempsum = 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]):
t, s12, t, t, t, t, t, 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
crossings += PolygonArea.transit(self._lon1 if i == 0 else lon,
self._lon0 if i != 0 else lon)
if self._polyline:
return num, perimeter, area
if crossings & 1:
tempsum += (1 if tempsum < 0 else -1) * self._area0/2
# area is with the clockwise sense. If !reverse convert to
# counter-clockwise convention.
if not reverse: tempsum *= -1
# If sign put area in (-area0/2, area0/2], else put area in [0, area0)
if sign:
if tempsum > self._area0/2:
tempsum -= self._area0
elif tempsum <= -self._area0/2:
tempsum += self._area0
else:
if tempsum >= self._area0:
tempsum -= self._area0
elif tempsum < 0:
tempsum += self._area0
area = 0 + tempsum
return num, perimeter, area
def Area(earth, points, polyline):
poly = PolygonArea(earth, polyline)
for p in points:
poly.AddPoint(p['lat'], p['lon'])
return poly.Compute(False, True)
Area = staticmethod(Area)
|