BUG: Fix handling of polygon holes when calculating area in Geod (#686)

pyproj4 · Jul 22, 2020 · 1781260 · 1781260
1 parent 9222735
commit 1781260
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 12 deletions.
diff --git a/pyproj/geod.py b/pyproj/geod.py
@@ -530,9 +530,11 @@ def geometry_area_perimeter(
 
         .. note:: lats should be in the range [-90 deg, 90 deg].
 
-        .. warning:: The area returned is signed with counter-clockwise traversal being
-                     treated as positive. You can use `shapely.ops.orient` to modify the
-                     orientation.
+        .. warning:: The area returned is signed with counter-clockwise (CCW) traversal
+                     being treated as positive. For polygons, holes should use the
+                     opposite traversal to the exterior (if the exterior is CCW, the
+                     holes/interiors should be CW). You can use `shapely.ops.orient` to
+                     modify the orientation.
 
         If it is a Polygon, it will return the area and exterior perimeter.
         It will subtract the area of the interior holes.
@@ -548,13 +550,13 @@ def geometry_area_perimeter(
         >>> poly_area, poly_perimeter = geod.geometry_area_perimeter(
         ...     Polygon(
         ...         LineString([
-        ...             Point(1, 1), Point(1, 10), Point(10, 10), Point(10, 1)
+        ...             Point(1, 1), Point(10, 1), Point(10, 10), Point(1, 10)
         ...         ]),
         ...         holes=[LineString([Point(1, 2), Point(3, 4), Point(5, 2)])],
         ...     )
         ... )
         >>> f"{poly_area:.3f} {poly_perimeter:.3f}"
-        '-944373881400.339 3979008.036'
+        '944373881400.339 3979008.036'
 
 
         Parameters
@@ -583,7 +585,7 @@ def geometry_area_perimeter(
             # subtract area of holes
             for hole in geometry.interiors:
                 area, _ = self.geometry_area_perimeter(hole, radians=radians)
-                total_area -= area
+                total_area += area
             return total_area, total_perimeter
         # multi geometries
         elif hasattr(geometry, "geoms"):

diff --git a/test/test_geod.py b/test/test_geod.py
@@ -21,6 +21,7 @@
         Point,
         Polygon,
     )
+    from shapely.geometry.polygon import orient
 
     SHAPELY_LOADED = True
 except ImportError:
@@ -380,13 +381,20 @@ def test_geometry_area_perimeter__polygon__radians():
 @skip_shapely
 def test_geometry_area_perimeter__polygon__holes():
     geod = Geod(ellps="WGS84")
+
+    polygon = Polygon(
+        LineString([Point(1, 1), Point(1, 10), Point(10, 10), Point(10, 1)]),
+        holes=[LineString([Point(1, 2), Point(3, 4), Point(5, 2)])],
+    )
+
     assert_almost_equal(
-        geod.geometry_area_perimeter(
-            Polygon(
-                LineString([Point(1, 1), Point(1, 10), Point(10, 10), Point(10, 1)]),
-                holes=[LineString([Point(1, 2), Point(3, 4), Point(5, 2)])],
-            )
-        ),
+        geod.geometry_area_perimeter(orient(polygon, 1)),
+        (944373881400.3394, 3979008.0359657984),
+        decimal=2,
+    )
+
+    assert_almost_equal(
+        geod.geometry_area_perimeter(orient(polygon, -1)),
         (-944373881400.3394, 3979008.0359657984),
         decimal=2,
     )