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

编辑:

点很可能不是以任何指定的顺序(顺时针或逆时针),并且可能非常复杂,因为它们是来自一组边界下的shapefile的一组utm坐标

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笔记本中完成.

Answer 2

您可以使用鞋带配方,例如

def PolygonArea(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

# examples
corners = [(2.0, 1.0), (4.0, 5.0), (7.0, 8.0)]

^{这仅适用于简单的多边形}

• 如果您有一个带孔的多边形:计算外圈的面积并减去内圈的面积

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

Answer 3

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. 它在每个点积中产生正值和负值的混合，这将在很大程度上抵消。

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

Answer 4

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], ...]

Answer 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

Answer 6

通过分析马赫迪的回答，我得出结论，大部分时间都花在了做上np.roll()。通过消除滚动的需要，并且仍然使用numpy，与Mahdi的41µs（相比之下，Mahdi的功能在我的机器上平均花费37µs）相比，我使每个循环的执行时间降至4-5µ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()

Answer 7

shapely.geometry.Polygon使用起来比自己计算要快。

\n

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

\n

使用这些代码，并执行%timeit\xef\xbc\x9a

\n

%timeit PolyArea(x,y)\n46.4 \xc2\xb5s \xc2\xb1 2.24 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10000 loops each)\n%timeit Polygon(coords).area\n20.2 \xc2\xb5s \xc2\xb1 414 ns per loop (mean \xc2\xb1 std. dev. of 7 runs, 100000 loops each)\n

\n