polygonarea.py
Go to the documentation of this file.
1 """Define the :class:`~geographiclib.polygonarea.PolygonArea` class
2 
3 The constructor initializes a empty polygon. The available methods are
4 
5  * :meth:`~geographiclib.polygonarea.PolygonArea.Clear` reset the
6  polygon
7  * :meth:`~geographiclib.polygonarea.PolygonArea.AddPoint` add a vertex
8  to the polygon
9  * :meth:`~geographiclib.polygonarea.PolygonArea.AddEdge` add an edge
10  to the polygon
11  * :meth:`~geographiclib.polygonarea.PolygonArea.Compute` compute the
12  properties of the polygon
13  * :meth:`~geographiclib.polygonarea.PolygonArea.TestPoint` compute the
14  properties of the polygon with a tentative additional vertex
15  * :meth:`~geographiclib.polygonarea.PolygonArea.TestEdge` compute the
16  properties of the polygon with a tentative additional edge
17 
18 The public attributes for this class are
19 
20  * :attr:`~geographiclib.polygonarea.PolygonArea.earth`
21  :attr:`~geographiclib.polygonarea.PolygonArea.polyline`
22  :attr:`~geographiclib.polygonarea.PolygonArea.area0`
23  :attr:`~geographiclib.polygonarea.PolygonArea.num`
24  :attr:`~geographiclib.polygonarea.PolygonArea.lat1`
25  :attr:`~geographiclib.polygonarea.PolygonArea.lon1`
26 
27 """
28 # polygonarea.py
29 #
30 # This is a rather literal translation of the GeographicLib::PolygonArea class
31 # to python. See the documentation for the C++ class for more information at
32 #
33 # https://geographiclib.sourceforge.io/html/annotated.html
34 #
35 # The algorithms are derived in
36 #
37 # Charles F. F. Karney,
38 # Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
39 # https://doi.org/10.1007/s00190-012-0578-z
40 # Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
41 #
42 # Copyright (c) Charles Karney (2011-2017) <charles@karney.com> and licensed
43 # under the MIT/X11 License. For more information, see
44 # https://geographiclib.sourceforge.io/
45 
46 
47 import math
48 from geographiclib.geomath import Math
49 from geographiclib.accumulator import Accumulator
50 
52  """Area of a geodesic polygon"""
53 
54  def _transit(lon1, lon2):
55  """Count crossings of prime meridian for AddPoint."""
56  # Return 1 or -1 if crossing prime meridian in east or west direction.
57  # Otherwise return zero.
58  # Compute lon12 the same way as Geodesic::Inverse.
59  lon1 = Math.AngNormalize(lon1)
60  lon2 = Math.AngNormalize(lon2)
61  lon12, _ = Math.AngDiff(lon1, lon2)
62  cross = (1 if lon1 <= 0 and lon2 > 0 and lon12 > 0
63  else (-1 if lon2 <= 0 and lon1 > 0 and lon12 < 0 else 0))
64  return cross
65  _transit = staticmethod(_transit)
66 
67  def _transitdirect(lon1, lon2):
68  """Count crossings of prime meridian for AddEdge."""
69  # We want to compute exactly
70  # int(floor(lon2 / 360)) - int(floor(lon1 / 360))
71  # Since we only need the parity of the result we can use std::remquo but
72  # this is buggy with g++ 4.8.3 and requires C++11. So instead we do
73  lon1 = math.fmod(lon1, 720.0); lon2 = math.fmod(lon2, 720.0)
74  return ( (0 if ((lon2 >= 0 and lon2 < 360) or lon2 < -360) else 1) -
75  (0 if ((lon1 >= 0 and lon1 < 360) or lon1 < -360) else 1) )
76  _transitdirect = staticmethod(_transitdirect)
77 
78  def __init__(self, earth, polyline = False):
79  """Construct a PolygonArea object
80 
81  :param earth: a :class:`~geographiclib.geodesic.Geodesic` object
82  :param polyline: if true, treat object as a polyline instead of a polygon
83 
84  Initially the polygon has no vertices.
85  """
86 
87  from geographiclib.geodesic import Geodesic
88  self.earth = earth
89  """The geodesic object (readonly)"""
90  self.polyline = polyline
91  """Is this a polyline? (readonly)"""
92  self.area0 = 4 * math.pi * earth._c2
93  """The total area of the ellipsoid in meter^2 (readonly)"""
94  self._mask = (Geodesic.LATITUDE | Geodesic.LONGITUDE |
95  Geodesic.DISTANCE |
96  (Geodesic.EMPTY if self.polyline else
97  Geodesic.AREA | Geodesic.LONG_UNROLL))
98  if not self.polyline: self._areasum = Accumulator()
100  self.num = 0
101  """The current number of points in the polygon (readonly)"""
102  self.lat1 = Math.nan
103  """The current latitude in degrees (readonly)"""
104  self.lon1 = Math.nan
105  """The current longitude in degrees (readonly)"""
106  self.Clear()
107 
108  def Clear(self):
109  """Reset to empty polygon."""
110  self.num = 0
111  self._crossings = 0
112  if not self.polyline: self._areasum.Set(0)
113  self._perimetersum.Set(0)
114  self._lat0 = self._lon0 = self.lat1 = self.lon1 = Math.nan
115 
116  def AddPoint(self, lat, lon):
117  """Add the next vertex to the polygon
118 
119  :param lat: the latitude of the point in degrees
120  :param lon: the longitude of the point in degrees
121 
122  This adds an edge from the current vertex to the new vertex.
123  """
124 
125  if self.num == 0:
126  self._lat0 = self.lat1 = lat
127  self._lon0 = self.lon1 = lon
128  else:
129  _, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
130  self.lat1, self.lon1, lat, lon, self._mask)
131  self._perimetersum.Add(s12)
132  if not self.polyline:
133  self._areasum.Add(S12)
134  self._crossings += PolygonArea._transit(self.lon1, lon)
135  self.lat1 = lat
136  self.lon1 = lon
137  self.num += 1
138 
139  def AddEdge(self, azi, s):
140  """Add the next edge to the polygon
141 
142  :param azi: the azimuth at the current the point in degrees
143  :param s: the length of the edge in meters
144 
145  This specifies the new vertex in terms of the edge from the current
146  vertex.
147 
148  """
149 
150  if self.num != 0:
151  _, lat, lon, _, _, _, _, _, S12 = self.earth._GenDirect(
152  self.lat1, self.lon1, azi, False, s, self._mask)
153  self._perimetersum.Add(s)
154  if not self.polyline:
155  self._areasum.Add(S12)
156  self._crossings += PolygonArea._transitdirect(self.lon1, lon)
157  self.lat1 = lat
158  self.lon1 = lon
159  self.num += 1
160 
161  # return number, perimeter, area
162  def Compute(self, reverse = False, sign = True):
163  """Compute the properties of the polygon
164 
165  :param reverse: if true then clockwise (instead of
166  counter-clockwise) traversal counts as a positive area
167  :param sign: if true then return a signed result for the area if the
168  polygon is traversed in the "wrong" direction instead of returning
169  the area for the rest of the earth
170  :return: a tuple of number, perimeter (meters), area (meters^2)
171 
172  If the object is a polygon (and not a polygon), the perimeter
173  includes the length of a final edge connecting the current point to
174  the initial point. If the object is a polyline, then area is nan.
175 
176  More points can be added to the polygon after this call.
177 
178  """
179  if self.polyline: area = Math.nan
180  if self.num < 2:
181  perimeter = 0.0
182  if not self.polyline: area = 0.0
183  return self.num, perimeter, area
184 
185  if self.polyline:
186  perimeter = self._perimetersum.Sum()
187  return self.num, perimeter, area
188 
189  _, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
190  self.lat1, self.lon1, self._lat0, self._lon0, self._mask)
191  perimeter = self._perimetersum.Sum(s12)
192  tempsum = Accumulator(self._areasum)
193  tempsum.Add(S12)
194  crossings = self._crossings + PolygonArea._transit(self.lon1, self._lon0)
195  if crossings & 1:
196  tempsum.Add( (1 if tempsum.Sum() < 0 else -1) * self.area0/2 )
197  # area is with the clockwise sense. If !reverse convert to
198  # counter-clockwise convention.
199  if not reverse: tempsum.Negate()
200  # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
201  if sign:
202  if tempsum.Sum() > self.area0/2:
203  tempsum.Add( -self.area0 )
204  elif tempsum.Sum() <= -self.area0/2:
205  tempsum.Add( self.area0 )
206  else:
207  if tempsum.Sum() >= self.area0:
208  tempsum.Add( -self.area0 )
209  elif tempsum.Sum() < 0:
210  tempsum.Add( self.area0 )
211 
212  area = 0.0 + tempsum.Sum()
213  return self.num, perimeter, area
214 
215  # return number, perimeter, area
216  def TestPoint(self, lat, lon, reverse = False, sign = True):
217  """Compute the properties for a tentative additional vertex
218 
219  :param lat: the latitude of the point in degrees
220  :param lon: the longitude of the point in degrees
221  :param reverse: if true then clockwise (instead of
222  counter-clockwise) traversal counts as a positive area
223  :param sign: if true then return a signed result for the area if the
224  polygon is traversed in the "wrong" direction instead of returning
225  the area for the rest of the earth
226  :return: a tuple of number, perimeter (meters), area (meters^2)
227 
228  """
229  if self.polyline: area = Math.nan
230  if self.num == 0:
231  perimeter = 0.0
232  if not self.polyline: area = 0.0
233  return 1, perimeter, area
234 
235  perimeter = self._perimetersum.Sum()
236  tempsum = 0.0 if self.polyline else self._areasum.Sum()
237  crossings = self._crossings; num = self.num + 1
238  for i in ([0] if self.polyline else [0, 1]):
239  _, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
240  self.lat1 if i == 0 else lat, self.lon1 if i == 0 else lon,
241  self._lat0 if i != 0 else lat, self._lon0 if i != 0 else lon,
242  self._mask)
243  perimeter += s12
244  if not self.polyline:
245  tempsum += S12
246  crossings += PolygonArea._transit(self.lon1 if i == 0 else lon,
247  self._lon0 if i != 0 else lon)
248 
249  if self.polyline:
250  return num, perimeter, area
251 
252  if crossings & 1:
253  tempsum += (1 if tempsum < 0 else -1) * self.area0/2
254  # area is with the clockwise sense. If !reverse convert to
255  # counter-clockwise convention.
256  if not reverse: tempsum *= -1
257  # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
258  if sign:
259  if tempsum > self.area0/2:
260  tempsum -= self.area0
261  elif tempsum <= -self.area0/2:
262  tempsum += self.area0
263  else:
264  if tempsum >= self.area0:
265  tempsum -= self.area0
266  elif tempsum < 0:
267  tempsum += self.area0
268 
269  area = 0.0 + tempsum
270  return num, perimeter, area
271 
272  # return num, perimeter, area
273  def TestEdge(self, azi, s, reverse = False, sign = True):
274  """Compute the properties for a tentative additional edge
275 
276  :param azi: the azimuth at the current the point in degrees
277  :param s: the length of the edge in meters
278  :param reverse: if true then clockwise (instead of
279  counter-clockwise) traversal counts as a positive area
280  :param sign: if true then return a signed result for the area if the
281  polygon is traversed in the "wrong" direction instead of returning
282  the area for the rest of the earth
283  :return: a tuple of number, perimeter (meters), area (meters^2)
284 
285  """
286 
287  if self.num == 0: # we don't have a starting point!
288  return 0, Math.nan, Math.nan
289  num = self.num + 1
290  perimeter = self._perimetersum.Sum() + s
291  if self.polyline:
292  return num, perimeter, Math.nan
293 
294  tempsum = self._areasum.Sum()
295  crossings = self._crossings
296  _, lat, lon, _, _, _, _, _, S12 = self.earth._GenDirect(
297  self.lat1, self.lon1, azi, False, s, self._mask)
298  tempsum += S12
299  crossings += PolygonArea._transitdirect(self.lon1, lon)
300  _, s12, _, _, _, _, _, _, _, S12 = self.earth._GenInverse(
301  lat, lon, self._lat0, self._lon0, self._mask)
302  perimeter += s12
303  tempsum += S12
304  crossings += PolygonArea._transit(lon, self._lon0)
305 
306  if crossings & 1:
307  tempsum += (1 if tempsum < 0 else -1) * self.area0/2
308  # area is with the clockwise sense. If !reverse convert to
309  # counter-clockwise convention.
310  if not reverse: tempsum *= -1
311  # If sign put area in (-area0/2, area0/2], else put area in [0, area0)
312  if sign:
313  if tempsum > self.area0/2:
314  tempsum -= self.area0
315  elif tempsum <= -self.area0/2:
316  tempsum += self.area0
317  else:
318  if tempsum >= self.area0:
319  tempsum -= self.area0
320  elif tempsum < 0:
321  tempsum += self.area0
322 
323  area = 0.0 + tempsum
324  return num, perimeter, area
def Compute(self, reverse=False, sign=True)
Definition: polygonarea.py:162
def TestPoint(self, lat, lon, reverse=False, sign=True)
Definition: polygonarea.py:216
def TestEdge(self, azi, s, reverse=False, sign=True)
Definition: polygonarea.py:273
def __init__(self, earth, polyline=False)
Definition: polygonarea.py:78


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:14