计算给定（x，y）坐标的多边形面积

Posted 2023-02-14

【中文标题】计算给定（x，y）坐标的多边形面积

【英文标题】：Calculate area of polygon given (x,y) coordinates

【发布时间】：2014-08-19 12:24:25

【问题描述】：

我有一组点，想知道是否有一个函数（为了方便和可能的速度）可以计算一组点所包围的面积。

例如：

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 坐标的任何负值......它们将是没有任何定义函数的多边形。

编辑：

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

【问题讨论】：

谷歌首发：people.virginia.edu/~ll2bf/docs/various/polyarea.html

这里给出了一个更有效的公式：softsurfer.com/Archive/algorithm_0101/…。 Python 实现：***.com/a/4682656/190597.

【参考方案1】：

Shoelace formula 的实现可以在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 notebook 中完成。

【讨论】：

很好的解决方案。我不知道为什么，但是当某些坐标为负时，@Nikos Athanasiou 的“最佳”答案不起作用。列出的另一个解决方案here 也有这个问题。您的解决方案是唯一有效的解决方案。只需检查xxx = np.array([[-100,0],[100,0],[100,150],[-100,150],[-100,0]])

@user989762：但是我使用这两种方法得到了相同的答案！

新手错误：不按顺序（顺时针/逆时针）方式提供积分会产生错误的结果。

你能解释一下你是如何使用点积而不是鞋带论坛所说的叉积吗？

@pstatix：鞋带公式确实可以写成外部产品，但是你可以扩展产品，你会看到有两种类型的术语：正术语和负术语。如果将它们分成两个项，您会看到它们是 x 和 y 的乘积，那么您可以将这些 x 和 y 写为两个向量，它们之间有一个点积。在此处查看proof for a triangle 部分：en.wikipedia.org/wiki/Shoelace_formula

【参考方案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)

---

要从头开始构建它或了解底层算法的工作原理，请查看shoelace formula：

# 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

因为这适用于简单的多边形：

如果你有一个带孔的多边形：计算外环的面积并子跟踪内环的面积

如果你有自相交环：你必须将它们分解成简单的扇区

【讨论】：

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

@user2593236：只要你的多边形边界不交叉（在这种情况下，这就是“简单”的意思），你应该没问题。

@user2593236 Simple 表示没有孔或自相交的凹凸。

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

@diegopso 似乎只有当点在一系列绘图中时才有效。所以它适用于[(0, 0), (0, 1), (1, 1), (1, 0)]

【参考方案3】：

通过分析 Mahdi 的回答，我得出结论，大部分时间都花在了 np.roll() 上。通过消除滚动的需要，并且仍然使用 numpy，与 Mahdi 的 41μs 相比，我将每个循环的执行时间降低到 4-5μ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() 完成计时

【讨论】：

当我定义x 和y 如x_n+1 = x_1 and x_0 = x_n, as well as y_n+1 = y_1 and y_0 = y_n 应用鞋带公式（见en.wikipedia.org/wiki/Shoelace_formula#Definition）时，我得到了这种方法和Mahdi 方法的区别是轻微的，因为这些点是顶点彼此如此接近但存在并且在处理具有较长边的多边形时可能会被放大。

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

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

@Eskapp 很高兴听到一切正常！

@pstatix 如果您查看 Wikipedia 文章中的 Shoelace formula，可以将其可视化为移位点积。我自己没有想出公式，但我确实意识到使用的计算模式直接匹配使用点积（或者更确切地说是两个点积），每个乘积中的一个向量移动。有关更多信息，我刚刚阅读了这篇文章，我为这个答案所做的唯一事情就是提高算法的性能。

【参考方案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) 的位置，例如通过减去上面代码中的质心。这有两个方面的帮助：

它从每个产品中消除了x.mean() * y.mean() 的因子 它会在每个点积中产生正负值的混合，这将在很大程度上抵消。

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

【讨论】：

唯一提供正确结果的解决方案！赞！请参阅我的答案，以获取带有元组列表的稍作修改的版本。

【参考方案5】：

上面的代码有一个错误，因为它在每次迭代中都没有采用绝对值。上面的代码将始终返回零。 （在数学上，它是带符号面积或楔积与实际面积之间的差异 http://en.wikipedia.org/wiki/Exterior_algebra.) 这是一些备用代码。

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

【参考方案6】：

这里有点晚了，但是您是否考虑过简单地使用sympy？

一个简单的代码是：

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

【参考方案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)

【讨论】：

numpy 相当标准，但 shapely 更快

【参考方案8】：

我将这里提供的每个解决方案都与 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

【参考方案9】：

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

示例：

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

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

【讨论】：

你这里没有提到它适用于整数数组

这实际上似乎是最快的，至少对于我测试的简单多边形而言

【参考方案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

【讨论】：

有趣。似乎用 numba 进行 jit 编译会快速且容易。你有这方面的参考吗？

# 给定边数 n 和边长 s，多边形的面积为 # 1/4 n s2 / tan( pi/n) Interactive Python (Rice University, Coursera ) 再次在这里：多边形的面积 (academia.edu/5179705/Exercise_1_How_to_design_programs) 我从那里完成了这个功能......

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

【参考方案11】：

基于

https://www.mathsisfun.com/geometry/area-irregular-polygons.html

def _area_(coords):
    t=0
    for count in range(len(coords)-1):
        y = coords[count+1][1] + coords[count][1]
        x = coords[count+1][0] - coords[count][0]
        z = y * x
        t += z
    return abs(t/2.0)

a=[(5.09,5.8), (1.68,4.9), (1.48,1.38), (4.76,0.1), (7.0,2.83), (5.09,5.8)]
print _area_(a)

诀窍是第一个坐标也应该是最后一个。

【讨论】：

当我尝试具有 15 个顶点的更复杂区域时，它给出了错误的结果。

能否提供坐标？

对不起，这是我的错。我对您的代码进行了几次测试并将结果与​​ CAD 软件进行了比较，我测试了 coords=[(1141.784,893.124), (1521.933,893.124), (1521.933,999.127), (1989.809,999.127), (1989.809,622.633), (2125.054,622.633), (2125.054,326.556), (1372.067,326.556), (1372.067,-60.903), (1872.84,-60.903), (1872.84,52.41), (2015.396,52.41), (2015.396,52.41), , (1090.611,-455.673), (1086.955,436.214), (1141.784,893.124)] 昨天我得到了错误的结果，也许我错过了一些东西，今天它像 PolygonArea 函数一样工作得很好。

我想我评论错了，也许我昨天在这里尝试了另一个功能。

很高兴我能帮上忙