python3-geographiclib 1.45-2 / usr / lib / python3 / dist-packages / geographiclib / polygonarea.py

This file is indexed.

/usr/lib/python3/dist-packages/geographiclib/polygonarea.py is in python3-geographiclib 1.45-2.

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
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
"""polygonarea.py: transcription of GeographicLib::PolygonArea class."""
# 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
#
# The algorithms are derived in
#
#    Charles F. F. Karney,
#    Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
#    https://dx.doi.org/10.1007/s00190-012-0578-z
#    Addenda: http://geographiclib.sf.net/geod-addenda.html
#
# Copyright (c) Charles Karney (2011-2014) <charles@karney.com> and licensed
# under the MIT/X11 License.  For more information, see
# http://geographiclib.sourceforge.net/
######################################################################

import math
from geographiclib.geomath import Math
from geographiclib.accumulator import Accumulator

class PolygonArea(object):
  """Area of a geodesic polygon"""

  def transit(lon1, lon2):
    """Count crossings of prime meridian for AddPoint."""
    # Return 1 or -1 if crossing prime meridian in east or west direction.
    # Otherwise return zero.
    # Compute lon12 the same way as Geodesic::Inverse.
    lon1 = Math.AngNormalize(lon1)
    lon2 = Math.AngNormalize(lon2)
    lon12 = Math.AngDiff(lon1, lon2)
    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 transitdirect(lon1, lon2):
    """Count crossings of prime meridian for AddEdge."""
    # We want to compute exactly
    #   int(floor(lon2 / 360)) - int(floor(lon1 / 360))
    # Since we only need the parity of the result we can use std::remquo but
    # this is buggy with g++ 4.8.3 and requires C++11.  So instead we do
    lon1 = math.fmod(lon1, 720.0); lon2 = math.fmod(lon2, 720.0)
    return ( (0 if ((lon2 >= 0 and lon2 < 360) or lon2 < -360) else 1) -
             (0 if ((lon1 >= 0 and lon1 < 360) or lon1 < -360) else 1) )
  transitdirect = staticmethod(transitdirect)

  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.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.Clear()

  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 AddPoint(self, lat, lon):
    """Add a vertex to the polygon."""
    if self._num == 0:
      self._lat0 = self._lat1 = lat
      self._lon0 = self._lon1 = lon
    else:
      _, s12, _, _, _, _, _, 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

  def AddEdge(self, azi, s):
    """Add an edge to the polygon."""
    if self._num != 0:
      _, lat, lon, _, _, _, _, _, S12 = self._earth.GenDirect(
        self._lat1, self._lon1, azi, False, s, self._mask)
      self._perimetersum.Add(s)
      if not self._polyline:
        self._areasum.Add(S12)
        self._crossings += PolygonArea.transitdirect(self._lon1, lon)
      self._lat1 = lat
      self._lon1 = lon
      self._num += 1

  # return number, perimeter, area
  def Compute(self, reverse, sign):
    """Return the number, perimeter, and area for the polygon."""
    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

    _, s12, _, _, _, _, _, 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.Sum() < 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 TestPoint(self, lat, lon, reverse, sign):
    """Return the results for a tentative additional vertex."""
    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]):
      _, 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
        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

  # return number, perimeter, area (for backward compatibility)
  def TestCompute(self, lat, lon, reverse, sign):
    return self.TestPoint(lat, lon, reverse, sign)

  # return num, perimeter, area
  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.transitdirect(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)

    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 CurrentPoint(self):
    """Return the current point as a lat, lon tuple."""
    return self._lat1, self._lon1

  def Area(earth, points, polyline):
    """Return the number, perimeter, and area for a set of vertices."""
    poly = PolygonArea(earth, polyline)
    for p in points:
      poly.AddPoint(p['lat'], p['lon'])
    return poly.Compute(False, True)
  Area = staticmethod(Area)

APT Browse - Built by Thomas Orozco.