fix(cvfdutil): polygon area and centroid calculations now use shapely (…

…#2165)

autotest/test_grid.py

@@ -19,7 +19,12 @@
 from flopy.modflow import Modflow, ModflowDis
 from flopy.utils import import_optional_dependency
 from flopy.utils.crs import get_authority_crs
-from flopy.utils.cvfdutil import gridlist_to_disv_gridprops, to_cvfd
+from flopy.utils.cvfdutil import (
+    area_of_polygon,
+    centroid_of_polygon,
+    gridlist_to_disv_gridprops,
+    to_cvfd,
+)
 from flopy.utils.triangle import Triangle
 from flopy.utils.voronoi import VoronoiGrid
 
@@ -946,6 +951,28 @@ def test_tocvfd3():
         assert i == j, f"{i} not equal {j}"
 
 
+@requires_pkg("shapely")
+def test_area_centroid_polygon():
+    pts = [
+        (685053.450097303, 6295544.549730939),
+        (685055.8377391606, 6295545.167682521),
+        (685057.3028430222, 6295542.712221102),
+        (685055.3500302795, 6295540.907246565),
+        (685053.2040466429, 6295542.313082705),
+        (685053.450097303, 6295544.549730939),
+    ]
+    xc, yc = centroid_of_polygon(pts)
+    result = np.array([xc, yc])
+    answer = np.array((685055.1035824707, 6295543.12059913))
+    assert np.allclose(
+        result, answer
+    ), "cvfdutil centroid of polygon incorrect"
+    x, y = list(zip(*pts))
+    result = area_of_polygon(x, y)
+    answer = 11.228131838368032
+    assert np.allclose(result, answer), "cvfdutil area of polygon incorrect"
+
+
 def test_unstructured_grid_shell():
     # constructor with no arguments.  incomplete shell should exist
     g = UnstructuredGrid()

flopy/utils/cvfdutil.py

@@ -1,38 +1,22 @@
 import numpy as np
 
+from .utl_import import import_optional_dependency
+
 
 def area_of_polygon(x, y):
-    """Calculates the signed area of an arbitrary polygon given its vertices
-    https://stackoverflow.com/a/4682656/ (Joe Kington)
-    http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm#2D%20Polygons
-    """
-    area = 0.0
-    for i in range(-1, len(x) - 1):
-        area += x[i] * (y[i + 1] - y[i - 1])
-    return area / 2.0
+    shapely = import_optional_dependency("shapely")
+    from shapely.geometry import Polygon
+
+    pgon = Polygon(zip(x, y))
+    return pgon.area
 
 
 def centroid_of_polygon(points):
-    """
-    https://stackoverflow.com/a/14115494/ (mgamba)
-    """
-    import itertools as IT
-
-    area = area_of_polygon(*zip(*points))
-    result_x = 0
-    result_y = 0
-    N = len(points)
-    points = IT.cycle(points)
-    x1, y1 = next(points)
-    for i in range(N):
-        x0, y0 = x1, y1
-        x1, y1 = next(points)
-        cross = (x0 * y1) - (x1 * y0)
-        result_x += (x0 + x1) * cross
-        result_y += (y0 + y1) * cross
-    result_x /= area * 6.0
-    result_y /= area * 6.0
-    return (result_x, result_y)
+    shapely = import_optional_dependency("shapely")
+    from shapely.geometry import Polygon
+
+    pgon = Polygon(points)
+    return pgon.centroid.x, pgon.centroid.y
 
 
 class Point: