python-geographiclib 1.49-2 / usr / lib / python2.7 / dist-packages / geographiclib / polygonarea.py

This file is indexed.

/usr/lib/python2.7/dist-packages/geographiclib/polygonarea.py is in python-geographiclib 1.49-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
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
"""Define the :class:`~geographiclib.polygonarea.PolygonArea` class

The constructor initializes a empty polygon.  The available methods are

  * :meth:`~geographiclib.polygonarea.PolygonArea.Clear` reset the
    polygon
  * :meth:`~geographiclib.polygonarea.PolygonArea.AddPoint` add a vertex
    to the polygon
  * :meth:`~geographiclib.polygonarea.PolygonArea.AddEdge` add an edge
    to the polygon
  * :meth:`~geographiclib.polygonarea.PolygonArea.Compute` compute the
    properties of the polygon
  * :meth:`~geographiclib.polygonarea.PolygonArea.TestPoint` compute the
    properties of the polygon with a tentative additional vertex
  * :meth:`~geographiclib.polygonarea.PolygonArea.TestEdge` compute the
    properties of the polygon with a tentative additional edge

The public attributes for this class are

  * :attr:`~geographiclib.polygonarea.PolygonArea.earth`
    :attr:`~geographiclib.polygonarea.PolygonArea.polyline`
    :attr:`~geographiclib.polygonarea.PolygonArea.area0`
    :attr:`~geographiclib.polygonarea.PolygonArea.num`
    :attr:`~geographiclib.polygonarea.PolygonArea.lat1`
    :attr:`~geographiclib.polygonarea.PolygonArea.lon1`

"""
# 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
#
#    https://geographiclib.sourceforge.io/html/annotated.html
#
# The algorithms are derived in
#
#    Charles F. F. Karney,
#    Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
#    https://doi.org/10.1007/s00190-012-0578-z
#    Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
#
# Copyright (c) Charles Karney (2011-2017) <charles@karney.com> and licensed
# under the MIT/X11 License.  For more information, see
# https://geographiclib.sourceforge.io/
######################################################################

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):
    """Construct a PolygonArea object

    :param earth: a :class:`~geographiclib.geodesic.Geodesic` object
    :param polyline: if true, treat object as a polyline instead of a polygon

    Initially the polygon has no vertices.
    """

    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 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 the next vertex to the polygon

    :param lat: the latitude of the point in degrees
    :param lon: the longitude of the point in degrees

    This adds an edge from the current vertex to the new vertex.
    """

    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 the next edge to the polygon

    :param azi: the azimuth at the current the point in degrees
    :param s: the length of the edge in meters

    This specifies the new vertex in terms of the edge from the current
    vertex.

    """

    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 = False, sign = True):
    """Compute the properties of the polygon

    :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 the object is a polygon (and not a polygon), the perimeter
    includes the length of a final edge connecting the current point to
    the initial point.  If the object is a polyline, then area is nan.

    More points can be added to the polygon after this call.

    """
    if self.polyline: area = Math.nan
    if self.num < 2:
      perimeter = 0.0
      if not self.polyline: area = 0.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.0 + tempsum.Sum()
    return self.num, perimeter, area

  # return number, perimeter, area
  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
        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.0 + tempsum
    return num, perimeter, area

  # return num, perimeter, area
  def TestEdge(self, azi, s, reverse = False, sign = True):
    """Compute the properties for a tentative additional edge

    :param azi: the azimuth at the current the point in degrees
    :param s: the length of the edge in meters
    :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.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.0 + tempsum
    return num, perimeter, area

APT Browse - Built by Thomas Orozco.