Coverage for pygeodesy/geodesicx/gxarea.py: 95%
213 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2024-06-10 14:08 -0400
1# -*- coding: utf-8 -*-
3u'''Slightly enhanced versions of classes U{PolygonArea
4<https://GeographicLib.SourceForge.io/1.52/python/code.html#
5module-geographiclib.polygonarea>} and C{Accumulator} from
6I{Karney}'s Python U{geographiclib
7<https://GeographicLib.SourceForge.io/1.52/python/index.html>}.
9Class L{GeodesicAreaExact} is intended to work with instances
10of class L{GeodesicExact} and of I{wrapped} class C{Geodesic},
11see module L{pygeodesy.karney}.
13Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023)
14and licensed under the MIT/X11 License. For more information, see the
15U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.
16'''
17# make sure int/int division yields float quotient
18from __future__ import division as _; del _ # PYCHOK semicolon
20from pygeodesy.basics import isodd, unsigned0
21from pygeodesy.constants import NAN, _0_0, _0_5, _720_0
22# from pygeodesy.interns import _COMMASPACE_ # from .lazily
23from pygeodesy.karney import Area3Tuple, _diff182, GeodesicError, \
24 _norm180, _remainder, _sum2_
25from pygeodesy.lazily import _ALL_DOCS, printf, _COMMASPACE_
26from pygeodesy.named import ADict, callername, _NamedBase, pairs
27from pygeodesy.props import Property, Property_RO, property_RO
28# from pygeodesy.streprs import pairs # from .named
30from math import fmod as _fmod
32__all__ = ()
33__version__ = '24.05.19'
36class GeodesicAreaExact(_NamedBase):
37 '''Area and perimeter of a geodesic polygon, an enhanced
38 version of I{Karney}'s Python class U{PolygonArea
39 <https://GeographicLib.SourceForge.io/html/python/
40 code.html#module-geographiclib.polygonarea>} using
41 the more accurate surface area.
43 @note: The name of this class C{*Exact} is a misnomer, see
44 I{Karney}'s comments at C++ attribute U{GeodesicExact._c2
45 <https://GeographicLib.SourceForge.io/C++/doc/
46 GeodesicExact_8cpp_source.html>}.
47 '''
48 _Area = None
49 _g_gX = None # Exact or not
50 _lat0 = _lon0 = \
51 _lat1 = _lon1 = NAN
52 _mask = 0
53 _num = 0
54 _Peri = None
55 _verbose = False
56 _xings = 0
58 def __init__(self, geodesic, polyline=False, **name):
59 '''New L{GeodesicAreaExact} instance.
61 @arg geodesic: A geodesic (L{GeodesicExact}, I{wrapped}
62 C{Geodesic} or L{GeodesicSolve}).
63 @kwarg polyline: If C{True}, compute the perimeter only,
64 otherwise area and perimeter (C{bool}).
65 @kwarg name: Optional C{B{name}=NN} (C{str}).
67 @raise GeodesicError: Invalid B{C{geodesic}}.
68 '''
69 try: # results returned as L{GDict}
70 if not (callable(geodesic._GDictDirect) and
71 callable(geodesic._GDictInverse)):
72 raise TypeError()
73 except (AttributeError, TypeError):
74 raise GeodesicError(geodesic=geodesic)
76 self._g_gX = g = geodesic
77 # use the class-level Caps since the values
78 # differ between GeodesicExact and Geodesic
79 self._mask = g.DISTANCE | g.LATITUDE | g.LONGITUDE
80 self._Peri = _Accumulator(name='_Peri')
81 if not polyline: # perimeter and area
82 self._mask |= g.AREA | g.LONG_UNROLL
83 self._Area = _Accumulator(name='_Area')
84 if g.debug: # PYCHOK no cover
85 self.verbose = True # debug as verbosity
86 if name:
87 self.name = name
89 def AddEdge(self, azi, s):
90 '''Add another polygon edge.
92 @arg azi: Azimuth at the current point (compass
93 C{degrees360}).
94 @arg s: Length of the edge (C{meter}).
95 '''
96 if self.num < 1:
97 raise GeodesicError(num=self.num)
98 r = self._Direct(azi, s)
99 p = self._Peri.Add(s)
100 if self._Area:
101 a = self._Area.Add(r.S12)
102 self._xings += r.xing
103 else:
104 a = NAN
105 self._lat1 = r.lat2
106 self._lon1 = r.lon2
107 self._num += 1
108 if self.verbose: # PYCHOK no cover
109 self._print(self.num, p, a, r, lat1=r.lat2, lon1=r.lon2,
110 azi=azi, s=s)
111 return self.num
113 def AddPoint(self, lat, lon):
114 '''Add another polygon point.
116 @arg lat: Latitude of the point (C{degrees}).
117 @arg lon: Longitude of the point (C{degrees}).
118 '''
119 if self.num > 0:
120 r = self._Inverse(self.lat1, self.lon1, lat, lon)
121 s = r.s12
122 p = self._Peri.Add(s)
123 if self._Area:
124 a = self._Area.Add(r.S12)
125 self._xings += r.xing
126 else:
127 a = NAN
128 else:
129 self._lat0 = lat
130 self._lon0 = lon
131 a = p = s = _0_0
132 r = None
133 self._lat1 = lat
134 self._lon1 = lon
135 self._num += 1
136 if self.verbose: # PYCHOK no cover
137 self._print(self.num, p, a, r, lat1=lat, lon1=lon, s=s)
138 return self.num
140 @Property_RO
141 def area0x(self):
142 '''Get the ellipsoid's surface area (C{meter} I{squared}),
143 more accurate for very I{oblate} ellipsoids.
144 '''
145 return self.ellipsoid.areax # not .area!
147 area0 = area0x # for C{geographiclib} compatibility
149 def Compute(self, reverse=False, sign=True):
150 '''Compute the accumulated perimeter and area.
152 @kwarg reverse: If C{True}, clockwise traversal counts as a
153 positive area instead of counter-clockwise
154 (C{bool}).
155 @kwarg sign: If C{True}, return a signed result for the area if
156 the polygon is traversed in the "wrong" direction
157 instead of returning the area for the rest of the
158 earth.
160 @return: L{Area3Tuple}C{(number, perimeter, area)} with the
161 number of points, the perimeter in C{meter} and the
162 area in C{meter**2}. The perimeter includes the
163 length of a final edge, connecting the current to
164 the initial point, if this polygon was initialized
165 with C{polyline=False}. For perimeter only, i.e.
166 C{polyline=True}, area is C{NAN}.
168 @note: Arbitrarily complex polygons are allowed. In the case
169 of self-intersecting polygons, the area is accumulated
170 "algebraically". E.g., the areas of the 2 loops in a
171 I{figure-8} polygon will partially cancel.
173 @note: More points and edges can be added after this call.
174 '''
175 r, n = None, self.num
176 if n < 2:
177 p = _0_0
178 a = NAN if self.polyline else p
179 elif self._Area:
180 r = self._Inverse(self.lat1, self.lon1, self.lat0, self.lon0)
181 a = self._reduced(r.S12, reverse, sign, r.xing)
182 p = self._Peri.Sum(r.s12)
183 else:
184 p = self._Peri.Sum()
185 a = NAN
186 if self.verbose: # PYCHOK no cover
187 self._print(n, p, a, r, lat0=self.lat0, lon0=self.lon0)
188 return Area3Tuple(n, p, a)
190 def _Direct(self, azi, s):
191 '''(INTERNAL) Edge helper.
192 '''
193 lon1 = self.lon1
194 r = self._g_gX._GDictDirect(self.lat1, lon1, azi, False, s, self._mask)
195 if self._Area: # aka transitDirect
196 # Count crossings of prime meridian exactly as
197 # int(ceil(lon2 / 360)) - int(ceil(lon1 / 360))
198 # Since we only need the parity of the result we
199 # can use std::remquo but this is buggy with g++
200 # 4.8.3 and requires C++11. So instead we do:
201 lon1 = _fmod( lon1, _720_0) # r.lon1
202 lon2 = _fmod(r.lon2, _720_0)
203 # int(True) == 1, int(False) == 0
204 r.set_(xing=int(lon2 > 360 or -360 < lon2 <= 0) -
205 int(lon1 > 360 or -360 < lon1 <= 0))
206 return r
208 @Property_RO
209 def ellipsoid(self):
210 '''Get this area's ellipsoid (C{Ellipsoid[2]}).
211 '''
212 return self._g_gX.ellipsoid
214 @Property_RO
215 def geodesic(self):
216 '''Get this area's geodesic object (C{Geodesic[Exact]}).
217 '''
218 return self._g_gX
220 earth = geodesic # for C{geographiclib} compatibility
222 def _Inverse(self, lat1, lon1, lat2, lon2):
223 '''(INTERNAL) Point helper.
224 '''
225 r = self._g_gX._GDictInverse(lat1, lon1, lat2, lon2, self._mask)
226 if self._Area: # aka transit
227 # count crossings of prime meridian as +1 or -1
228 # if in east or west direction, otherwise 0
229 lon1 = _norm180(lon1)
230 lon2 = _norm180(lon2)
231 lon12, _ = _diff182(lon1, lon2)
232 r.set_(xing=int(lon12 > 0 and lon1 <= 0 and lon2 > 0) or
233 -int(lon12 < 0 and lon2 <= 0 and lon1 > 0))
234 return r
236 @property_RO
237 def lat0(self):
238 '''Get the first point's latitude (C{degrees}).
239 '''
240 return self._lat0
242 @property_RO
243 def lat1(self):
244 '''Get the most recent point's latitude (C{degrees}).
245 '''
246 return self._lat1
248 @property_RO
249 def lon0(self):
250 '''Get the first point's longitude (C{degrees}).
251 '''
252 return self._lon0
254 @property_RO
255 def lon1(self):
256 '''Get the most recent point's longitude (C{degrees}).
257 '''
258 return self._lon1
260 @property_RO
261 def num(self):
262 '''Get the current number of points (C{int}).
263 '''
264 return self._num
266 @Property_RO
267 def polyline(self):
268 '''Is this perimeter only (C{bool}), area NAN?
269 '''
270 return self._Area is None
272 def _print(self, n, p, a, r, **kwds): # PYCHOK no cover
273 '''(INTERNAL) Print a verbose line.
274 '''
275 d = ADict(p=p, s12=r.s12 if r else NAN, **kwds)
276 if self._Area:
277 d.set_(a=a, S12=r.S12 if r else NAN)
278 t = _COMMASPACE_.join(pairs(d, prec=10))
279 printf('%s %s: %s (%s)', self.named2, n, t, callername(up=2))
281 def _reduced(self, S12, reverse, sign, xing):
282 '''(INTERNAL) Accumulate and reduce area to allowed range.
283 '''
284 a0 = self.area0x
285 A = _Accumulator(self._Area)
286 _ = A.Add(S12)
287 a = A.Remainder(a0) # clockwise
288 if isodd(self._xings + xing):
289 a = A.Add((a0 if a < 0 else -a0) * _0_5)
290 if not reverse:
291 a = A.Negate() # counter-clockwise
292 # (-area0x/2, area0x/2] if sign else [0, area0x)
293 a0_ = a0 if sign else (a0 * _0_5)
294 if a > a0_:
295 a = A.Add(-a0)
296 elif a <= -a0_:
297 a = A.Add( a0)
298 return unsigned0(a)
300 def Reset(self):
301 '''Reset this polygon to empty.
302 '''
303 if self._Area:
304 self._Area.Reset()
305 self._Peri.Reset()
306 self._lat0 = self._lon0 = \
307 self._lat1 = self._lon1 = NAN
308 self._num = self._xings = n = 0
309 if self.verbose: # PYCHOK no cover
310 printf('%s %s: (%s)', self.named2, n, self.Reset.__name__)
311 return n
313 Clear = Reset
315 def TestEdge(self, azi, s, reverse=False, sign=True):
316 '''Compute the properties for a tentative, additional edge
318 @arg azi: Azimuth at the current the point (compass C{degrees}).
319 @arg s: Length of the edge (C{meter}).
320 @kwarg reverse: If C{True}, clockwise traversal counts as a
321 positive area instead of counter-clockwise
322 (C{bool}).
323 @kwarg sign: If C{True}, return a signed result for the area if
324 the polygon is traversed in the "wrong" direction
325 instead of returning the area for the rest of the
326 earth.
328 @return: L{Area3Tuple}C{(number, perimeter, area)}.
330 @raise GeodesicError: No points.
331 '''
332 n = self.num + 1
333 p = self._Peri.Sum(s)
334 if self.polyline:
335 a, r = NAN, None
336 elif n < 2:
337 raise GeodesicError(num=self.num)
338 else:
339 d = self._Direct(azi, s)
340 r = self._Inverse(d.lat2, d.lon2, self.lat0, self.lon0)
341 a = self._reduced(d.S12 + r.S12, reverse, sign, d.xing + r.xing)
342 p += r.s12
343 if self.verbose: # PYCHOK no cover
344 self._print(n, p, a, r, azi=azi, s=s)
345 return Area3Tuple(n, p, a)
347 def TestPoint(self, lat, lon, reverse=False, sign=True):
348 '''Compute the properties for a tentative, additional vertex
350 @arg lat: Latitude of the point (C{degrees}).
351 @arg lon: Longitude of the point (C{degrees}).
352 @kwarg reverse: If C{True}, clockwise traversal counts as a
353 positive area instead of counter-clockwise
354 (C{bool}).
355 @kwarg sign: If C{True}, return a signed result for the area if
356 the polygon is traversed in the "wrong" direction
357 instead of returning the area for the rest of the
358 earth.
360 @return: L{Area3Tuple}C{(number, perimeter, area)}.
361 '''
362 r, n = None, self.num + 1
363 if n < 2:
364 p = _0_0
365 a = NAN if self.polyline else p
366 else:
367 i = self._Inverse(self.lat1, self.lon1, lat, lon)
368 p = self._Peri.Sum(i.s12)
369 if self._Area:
370 r = self._Inverse(lat, lon, self.lat0, self.lon0)
371 a = self._reduced(i.S12 + r.S12, reverse, sign, i.xing + r.xing)
372 p += r.s12
373 else:
374 a = NAN
375 if self.verbose: # PYCHOK no cover
376 self._print(n, p, a, r, lat=lat, lon=lon)
377 return Area3Tuple(n, p, a)
379 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature
380 '''Return this C{GeodesicExactArea} as string.
382 @kwarg prec: The C{float} precision, number of decimal digits (0..9).
383 Trailing zero decimals are stripped for B{C{prec}} values
384 of 1 and above, but kept for negative B{C{prec}} values.
385 @kwarg sep: Separator to join (C{str}).
387 @return: Area items (C{str}).
388 '''
389 n, p, a = self.Compute()
390 d = dict(geodesic=self.geodesic, num=n, area=a,
391 perimeter=p, polyline=self.polyline)
392 return sep.join(pairs(d, prec=prec))
394 @Property
395 def verbose(self):
396 '''Get the C{verbose} option (C{bool}).
397 '''
398 return self._verbose
400 @verbose.setter # PYCHOK setter!
401 def verbose(self, verbose): # PYCHOK no cover
402 '''Set the C{verbose} option (C{bool}) to print
403 a message after each method invokation.
404 '''
405 self._verbose = bool(verbose)
408class PolygonArea(GeodesicAreaExact):
409 '''For C{geographiclib} compatibility, sub-class of L{GeodesicAreaExact}.
410 '''
411 def __init__(self, earth, polyline=False, **name):
412 '''New L{PolygonArea} instance.
414 @arg earth: A geodesic (L{GeodesicExact}, I{wrapped}
415 C{Geodesic} or L{GeodesicSolve}).
416 @kwarg polyline: If C{True} perimeter only, otherwise
417 area and perimeter (C{bool}).
418 @kwarg name: Optional C{B{name}=NN} (C{str}).
420 @raise GeodesicError: Invalid B{C{earth}}.
421 '''
422 GeodesicAreaExact.__init__(self, earth, polyline=polyline, **name)
425class _Accumulator(_NamedBase):
426 '''Like C{math.fsum}, but allowing a running sum.
428 Original from I{Karney}'s U{geographiclib
429 <https://PyPI.org/project/geographiclib>}C{.accumulator},
430 enhanced to return the current sum by most methods.
431 '''
432 _n = 0 # len()
433 _s = _t = _0_0
435 def __init__(self, y=0, **name):
436 '''New L{_Accumulator}.
438 @kwarg y: Initial value (C{scalar}).
439 @kwarg name: Optional C{B{name}=NN} (C{str}).
440 '''
441 if isinstance(y, _Accumulator):
442 self._s, self._t, self._n = y._s, y._t, 1
443 elif y:
444 self._s, self._n = float(y), 1
445 if name:
446 self.name = name
448 def Add(self, y):
449 '''Add a value.
451 @return: Current C{sum}.
452 '''
453 self._n += 1
454 self._s, self._t = _sum2_(self._s, self._t, y)
455 return self._s # current .Sum()
457 def Negate(self):
458 '''Negate sum.
460 @return: Current C{sum}.
461 '''
462 self._s = s = -self._s
463 self._t = -self._t
464 return s # current .Sum()
466 @property_RO
467 def num(self):
468 '''Get the current number of C{Add}itions (C{int}).
469 '''
470 return self._n
472 def Remainder(self, y):
473 '''Remainder on division by B{C{y}}.
475 @return: Remainder of C{sum} / B{C{y}}.
476 '''
477 self._s = _remainder(self._s, y)
478# self._t = _remainder(self._t, y)
479 self._n = -1
480 return self.Add(_0_0)
482 def Reset(self, y=0):
483 '''Set value from argument.
484 '''
485 self._n, self._s, self._t = 0, float(y), _0_0
487 Set = Reset
489 def Sum(self, y=0):
490 '''Return C{sum + B{y}}.
492 @note: B{C{y}} is included in the returned
493 result, but I{not} accumulated.
494 '''
495 if y:
496 s = _Accumulator(self, name='_Sum')
497 s.Add(y)
498 else:
499 s = self
500 return s._s
502 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature
503 '''Return this C{_Accumulator} as string.
505 @kwarg prec: The C{float} precision, number of decimal digits (0..9).
506 Trailing zero decimals are stripped for B{C{prec}} values
507 of 1 and above, but kept for negative B{C{prec}} values.
508 @kwarg sep: Separator to join (C{str}).
510 @return: Accumulator (C{str}).
511 '''
512 d = dict(n=self.num, s=self._s, t=self._t)
513 return sep.join(pairs(d, prec=prec))
516__all__ += _ALL_DOCS(GeodesicAreaExact, PolygonArea)
518# **) MIT License
519#
520# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.
521#
522# Permission is hereby granted, free of charge, to any person obtaining a
523# copy of this software and associated documentation files (the "Software"),
524# to deal in the Software without restriction, including without limitation
525# the rights to use, copy, modify, merge, publish, distribute, sublicense,
526# and/or sell copies of the Software, and to permit persons to whom the
527# Software is furnished to do so, subject to the following conditions:
528#
529# The above copyright notice and this permission notice shall be included
530# in all copies or substantial portions of the Software.
531#
532# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
533# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
534# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
535# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
536# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
537# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
538# OTHER DEALINGS IN THE SOFTWARE.