根据 (x,y) 坐标计算多边形面积-IT科技

摘要：问题描述：我有一组点，想知道是否有一个函数（为了方便和速度）可以计算一组点所包围的面积。例如：x = np.arange(0,1,0.001) y = np.sqrt(1-x**2) points = zip(x,y) 鉴于points面积应近似等于(pi-2)/4。也许 scipy、matplotlib、...

问题描述：

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

例如：

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 坐标

解决方案 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循环可使该函数比以下函数快约 50 倍PolygonArea：

%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 笔记本中完成的。

解决方案 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

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

如果你有一个带孔的多边形：计算外环的面积并对内环的面积进行子跟踪
如果你有自相交环：你必须将它们分解为简单的扇区

解决方案 3：

通过分析 Mahdi 的答案，我得出结论，大部分时间都花在了做上np.roll()。通过消除 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()

解决方案 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：

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

例子：

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

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

解决方案 6：

上面的代码有错误，因为它没有在每次迭代中取绝对值。上面的代码将始终返回零。（从数学上讲，这是取有符号面积或楔形积与实际面积之间的差值
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

解决方案 7：

有点晚了，但是你有没有考虑过简单地使用sympy？

一个简单的代码是：

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

解决方案 8：

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)

解决方案 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

解决方案 10：

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

import math

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

因为公式是¼ns2/tan(π/n)。给定边数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)

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

解决方案 12：

def find_int_coordinates(n: int, coords: list[list[int]]) -> float:
    rez = 0
    x, y = coords[n - 1]
    for coord in coords:
        rez += (x + coord[0]) * (y - coord[1])
        x, y = coord
    return abs(rez / 2)

解决方案 13：

鞋带公式来自根据多边形周围的连续点计算内部三角形。基于此解释编写代码可能会更有参考价值。多边形的面积就是其两个相邻向量的叉积（行列式）的长度除以二，然后相加，

https://youtu.be/0KjG8Pg6LGk?t=213

import numpy as np, numpy.linalg as lin

def area(pts):
    ps = np.array([0.5 * lin.det(np.vstack((pts[i], pts[i+1]))) for i in range(len(pts)-1)])
    s = np.sum(ps)
    p1,p2 = pts[-1],pts[0] # cycle back, last pt with the first 
    s += 0.5 * lin.det(np.vstack((p1,p2)))
    return np.abs(s)

points = np.array([[0,0],[10,0],[10,10],[0,10]])
area(points) # 100