计算由(x,y)坐标给出的多边形面积

Question

我有一组点，想知道是否有一种函数（为了方便和速度），可以计算由这组点围成的区域面积。

例如：

x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)

points = zip(x,y)

给定points，面积应该大致等于(pi-2)/4。也许可以用scipy、matplotlib、numpy、shapely等库来实现这个功能？我不会遇到x或y坐标为负数的情况……并且它们将是没有定义函数的多边形。

编辑：

points很可能不会按任何指定顺序（顺时针或逆时针）排列，并且可能非常复杂，因为它们是一个shapefile下一组边界内的utm坐标集合。

- pbreach

Answer 1

可以在Numpy中实现鞋带公式。假设这些顶点：

import numpy as np
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)

我们可以在numpy中重新定义函数来计算面积：

def PolyArea(x,y):
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))

获得结果：

print PolyArea(x,y)
# 0.26353377782163534

避免使用for循环，使得这个函数比PolygonArea函数快约50倍：

%timeit PolyArea(x,y)
# 10000 loops, best of 3: 42 µs per loop
%timeit PolygonArea(zip(x,y))
# 100 loops, best of 3: 2.09 ms per loop.

计时是在Jupyter笔记本中完成的。

- Mahdi

---

很棒的解决方案。我不确定为什么，但是@Nikos Athanasiou的“顶部”答案在一些坐标为负数时不起作用。此外，这里列出的另一个解决方案也有这个问题。您的解决方案是唯一有效的。只需使用xxx = np.array([[-100,0],[100,0],[100,150],[-100,150],[-100,0]])进行检查即可。 - user989762

使用两种方法得到了相同的答案！ - Mahdi

新手错误：未按顺序（顺时针/逆时针）提供点会导致错误的结果。 - emredog

你能解释一下你是如何使用点积而不是叉积作为鞋带公式的吗？ - pstatix

@pstatix：确实，鞋带公式可以用外积来表示，但你可以展开这个积，然后你会发现有两种类型的项：正项和负项。如果你将它们分成两个项，你会看到它们是x和y的乘积，然后你可以将这些x和y写成两个向量，并在它们之间进行点积。在这里查看“三角形证明”部分的证明：https://en.wikipedia.org/wiki/Shoelace_formula - Mahdi

Answer 2

覆盖所有可能情况的最优解决方案是使用一个几何学包，例如shapely、scikit-geometry或pygeos。它们都在底层使用C++几何学包。第一个可以通过pip轻松安装:

pip install shapely

而且简单易用：

from shapely.geometry import Polygon
pgon = Polygon(zip(x, y)) # Assuming the OP's x,y coordinates

print(pgon.area)

---

想要从头开始构建或了解底层算法的工作原理，请查看鞋带公式：

# e.g. corners = [(2.0, 1.0), (4.0, 5.0), (7.0, 8.0)]
def Area(corners):
    n = len(corners) # of corners
    area = 0.0
    for i in range(n):
        j = (i + 1) % n
        area += corners[i][0] * corners[j][1]
        area -= corners[j][0] * corners[i][1]
    area = abs(area) / 2.0
    return area

由于这个方法适用于简单的多边形：

• 如果您有一个带洞的多边形：计算外环的面积并减去内环的面积。

• 如果您有自相交的环：您需要将它们分解为简单的部分。

- Nikos Athanasiou

---

我的可能是非常复杂的多边形。这些点是在一组边界下从一个shapefile中选择的utm坐标。 - pbreach

只要你的多边形边界不相交（这就是此上下文中“简单”的意思），那么你应该没有问题。 - Mark Dickinson

“Simple” 的意思是凸多边形或凹多边形且没有空洞或自交。 - Nikos Athanasiou

哦，好的，我现在明白了。解释得很好。 - pbreach

我尝试了非常简单的坐标 [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)]，但是它给出了0.0的面积。您知道有什么限制吗？我还试图将其移出原点，但结果相同。 - diegopso

@diegopso 看起来这个只对点在一系列的图形中有效。所以对于 [(0, 0), (0, 1), (1, 1), (1, 0)] 这个可以使用。 - Pithikos

Answer 3

通过对Mahdi答案的分析，我得出结论：大部分时间都花费在执行np.roll()上。通过消除滚动的需要，仍然使用numpy，我将执行时间降低到每个循环4-5µs，相比之下，Mahdi的函数需要41µs（值得比较的是，Mahdi的函数在我的机器上平均需要37µs）。

def polygon_area(x,y):
    correction = x[-1] * y[0] - y[-1]* x[0]
    main_area = np.dot(x[:-1], y[1:]) - np.dot(y[:-1], x[1:])
    return 0.5*np.abs(main_area + correction)

通过计算校正项，然后切片数组，无需滚动或创建新的数组。

基准测试：

10000 iterations
PolyArea(x,y): 37.075µs per loop
polygon_area(x,y): 4.665µs per loop

使用time模块和time.clock()函数来计时。

- maxb

---

当我定义x和y时，我使用的方法和Mahdi的方法有所不同，例如 x_{n + 1} = x_1和x_0 = x_n，y_{n + 1} = y_1和y_0 = y_n，这是应用鞋带公式所必需的(请参阅https://en.wikipedia.org/wiki/Shoelace_formula#Definition)。差异很小，因为顶点非常靠近，但存在并且在处理边长更长的多边形时可能会被放大。 - Eskapp

当然，像任何实现一样，会存在浮点误差。你能提供一个完整的例子来说明这种差异吗？如果你需要更高的精度，可以使用任意精度算术。 - maxb

我的错，我对修正术语感到困惑，并认为我在跟踪代码中的错误时可能会出现一些差异。在比较计算多边形面积的不同实现进行了更多测试后，似乎完美地工作了。您的解决方案具有速度优势，并且易于阅读！ - Eskapp

@Eskapp很高兴听到一切都正常工作！ - maxb

你能解释一下为什么在某些情况下使用点积代替叉积吗？ - pstatix

@pstatix 如果你看一下维基百科上的鞋带公式，它可以被视为一个移位的点积。我并不是自己想出这个公式，但我确实意识到了计算模式直接匹配使用点积（或者说两个点积），每个乘积中的一个向量都被移动了。如果需要更多信息，我建议阅读文章，对于这个答案，我所做的唯一事情就是提高算法的性能。 - maxb

Answer 4

maxb的答案提供了良好的性能，但在坐标值或点数较大时很容易导致精度损失。可以通过简单的坐标移动来缓解这种情况：

def polygon_area(x,y):
    # coordinate shift
    x_ = x - x.mean()
    y_ = y - y.mean()
    # everything else is the same as maxb's code
    correction = x_[-1] * y_[0] - y_[-1]* x_[0]
    main_area = np.dot(x_[:-1], y_[1:]) - np.dot(y_[:-1], x_[1:])
    return 0.5*np.abs(main_area + correction)

例如，常见的地理参考系统是UTM，其可能具有（x，y）坐标（488685.984，7133035.984）。这两个值的乘积为3485814708748.448。您可以看到，这个单一乘积已经处于精度的边缘（它具有与输入相同的小数位数）。仅添加其中的一些乘积，更不用说成千上万个，都将导致精度损失。

缓解这种情况的简单方法是将多边形从大正数坐标移动到距离（0,0）更近的位置，例如通过减去质心，如上面的代码所示。这种方法有两个好处：

1. 它消除了每个乘积中的x.mean() * y.mean()一个因子
2. 它在每个点乘积中产生了一组正负值，这些值将在很大程度上相互抵消。

坐标偏移不会改变总面积，只会使计算更加数值稳定。

- Trenton

---

唯一提供正确结果的解决方案！赞！请查看我的答案，其中包含一个稍微修改过的版本，它接受元组列表作为参数。 - HumbleBee

Answer 5

在OpenCV中，cv2.contourArea()提供了另一种方法。

例如：

points = np.array([[0,0],[10,0],[10,10],[0,10]])
area = cv2.contourArea(points)
print(area) # 100.0

参数（在上面的例子中是points）是一个包含整数的numpy数组，其表示多边形的顶点：[[x1,y1],[x2,y2], ...]。

- Toru Kikuchi

---

您在这里没有提到它适用于整数数组。 - Eugene W.

这似乎是最快的，至少对于我测试过的简单多边形来说。 - crazjo

你可以使用cv2.contourArea(np.around(np.array([[pt] for pt in points])).astype(np.int32))。 - Coddy

Answer 6

上面的代码存在一个错误，因为它没有在每次迭代中使用绝对值。上面的代码将始终返回零。（数学上，这是取有符号面积或楔积与实际面积的差异。参见外代数。）这里是一些替代代码。

def area(vertices):
    n = len(vertices) # of corners
    a = 0.0
    for i in range(n):
        j = (i + 1) % n
        a += abs(vertices[i][0] * vertices[j][1]-vertices[j][0] * vertices[i][1])
    result = a / 2.0
    return result

- Chris Judge

Answer 7

使用shapely.geometry.Polygon比自行计算更快。

from shapely.geometry import Polygon
import numpy as np
def PolyArea(x,y):
    return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
coords = np.random.rand(6, 2)
x, y = coords[:, 0], coords[:, 1]

使用这些代码，并执行%timeit：

%timeit PolyArea(x,y)
46.4 µs ± 2.24 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit Polygon(coords).area
20.2 µs ± 414 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

- C.K.

---

numpy是相当标准的，但shapely更快一些。 - Jiadong

Answer 8

可能有点晚了，但您考虑过是否可以简单地使用 sympy 吗？

一个简单的代码如下：

from sympy import Polygon
a = Polygon((0, 0), (2, 0), (2, 2), (0, 2)).area
print(a)

- Alex Jiang

Answer 9

我将这里提供的每个解决方案与Shapely的面积方法结果进行了比较，它们具有正确的整数部分，但小数部分不同。只有@Trenton的解决方案提供了正确的结果。

现在对@Trenton的答案进行改进，以将坐标处理为元组列表，我得出了以下结论：

import numpy as np

def polygon_area(coords):
    # get x and y in vectors
    x = [point[0] for point in coords]
    y = [point[1] for point in coords]
    # shift coordinates
    x_ = x - np.mean(x)
    y_ = y - np.mean(y)
    # calculate area
    correction = x_[-1] * y_[0] - y_[-1] * x_[0]
    main_area = np.dot(x_[:-1], y_[1:]) - np.dot(y_[:-1], x_[1:])
    return 0.5 * np.abs(main_area + correction)

#### Example output
coords = [(385495.19520441635, 6466826.196947694), (385496.1951836388, 6466826.196947694), (385496.1951836388, 6466825.196929455), (385495.19520441635, 6466825.196929455), (385495.19520441635, 6466826.196947694)]

Shapely's area method:  0.9999974610685296
@Trenton's area method: 0.9999974610685296

- HumbleBee

Answer 10

对于正多边形，这要简单得多：

import math

def area_polygon(n, s):
    return 0.25 * n * s**2 / math.tan(math.pi/n)

由于公式为 ¼ n s2 / tan(π/n)。 给定边数n和每条边的长度s

- DrSocket

---

有趣。看起来这个可以很快很容易地使用numba进行jit编译。你有参考资料吗？ - pbreach

给定多边形的边数n和每边长度s，多边形的面积为1/4ns²/tan(π/n)。交互式Python（Rice大学，Coursera） 同样在这里： 多边形的面积 （http://www.academia.edu/5179705/Exercise_1_How_to_design_programs）我从那个链接中得到了这个函数... - DrSocket

这是针对一个“正”多边形的问题，它是这个问题的一种特殊但非常有限的情况。所有的边必须具有相同的长度（这也需要计算）。如果您解释了什么是“n”和“s”，那么可能会更加明显。 - pbreach