计算给定(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中实现Shoelace formula。假设这些顶点：

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中完成。

Answer 2

通过分析马赫迪的回答，我得出结论，大部分时间都花在了做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()完成的

Answer 3

上面的代码中有一个错误，因为它没有在每次迭代中获取绝对值。上面的代码总是返回零。(从数学上讲，这是取有符号面积或楔形乘积与实际面积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