Python 计算给定 (x,y) 坐标的多边形面积
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/24467972/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Calculate area of polygon given (x,y) coordinates
提问by pbreach
I have a set of points and would like to know if there is a function (for the sake of convenience and probably speed) that can calculate the area enclosed by a set of points.
我有一组点，想知道是否有一个函数（为了方便和速度）可以计算一组点所包围的面积。
for example:
例如：
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)
points = zip(x,y)
given pointsthe area should be approximately equal to (pi-2)/4. Maybe there is something from scipy, matplotlib, numpy, shapely, etc. to do this? I won't be encountering any negative values for either the x or y coordinates... and they will be polygons without any defined function.
鉴于points面积应约等于(pi-2)/4。也许有来自 scipy、matplotlib、numpy、shapely 等的东西来做到这一点？我不会遇到 x 或 y 坐标的任何负值……它们将是没有任何定义函数的多边形。
EDIT:
编辑：
points will most likely not be in any specified order (clockwise or counterclockwise) and may be quite complex as they are a set of utm coordinates from a shapefile under a set of boundaries
点很可能没有任何指定的顺序（顺时针或逆时针），并且可能非常复杂，因为它们是来自一组边界下的 shapefile 的一组 utm 坐标
采纳答案by Mahdi
Implementation of Shoelace formulacould be done in Numpy. Assuming these vertices:
Shoelace 公式的实现可以在Numpy. 假设这些顶点：
import numpy as np
x = np.arange(0,1,0.001)
y = np.sqrt(1-x**2)
We can redefine the function in numpy to find the area:
我们可以在 numpy 中重新定义函数来查找区域：
def PolyArea(x,y):
return 0.5*np.abs(np.dot(x,np.roll(y,1))-np.dot(y,np.roll(x,1)))
And getting results:
并得到结果：
print PolyArea(x,y)
# 0.26353377782163534
Avoiding forloop makes this function ~50X faster than PolygonArea:
避免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.
Timing is done in Jupyter notebook.
计时是在 Jupyter notebook 中完成的。
回答by Nikos Athanasiou
You can use the shoelace formula, eg
您可以使用鞋带公式，例如
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)]
This only works for simple polygons
这仅适用于简单的多边形
If you have a polygon with holes: Calculate the area of the outer ring and subtrack the areas of the inner rings
If you have self-intersecting rings: You have to decompose them into simple sectors
如果您有一个带孔的多边形：计算外环的面积并子跟踪内环的面积
如果您有自相交环：您必须将它们分解为简单的扇区
回答by Chris Judge
There's an error in the code above as it doesn't take absolute values on each iteration. The above code will always return zero. (Mathematically, it's the difference between taking signed area or wedge product and the actual area http://en.wikipedia.org/wiki/Exterior_algebra.) Here's some alternate code.
上面的代码中有一个错误，因为它在每次迭代中都没有取绝对值。上面的代码将始终返回零。（从数学上讲，这是带符号面积或楔形乘积与实际面积之间的区别 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
回答by Bizarre
This is much simpler, for regular polygons:
对于正多边形，这要简单得多：
import math
def area_polygon(n, s):
return 0.25 * n * s**2 / math.tan(math.pi/n)
since the formula is ? n s2 / tan(π/n). Given the number of sides, n, and the length of each side, s
因为公式是？n s2 / tan(π/n)。给定边数 n 和每边的长度 s
回答by Takis Tsiberis
Based on
基于
https://www.mathsisfun.com/geometry/area-irregular-polygons.html
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)
The trick is that the first coordinate should also be last.
诀窍是第一个坐标也应该是最后一个。
回答by maxb
By analysis of Mahdi's answer, I concluded that the majority of time was spent doing np.roll(). By removing the need of the roll, and still using numpy, I got the execution time down to 4-5μs per loop compared to Mahdi's 41μs (for comparison Mahdi's function took an average of 37μs on my machine).
通过分析 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)
By calculating the correctional term, and then slicing the arrays, there is no need to roll or create a new array.
通过计算校正项，然后对阵列进行切片，无需滚动或创建新阵列。
Benchmarks:
基准：
10000 iterations
PolyArea(x,y): 37.075μs per loop
polygon_area(x,y): 4.665μs per loop
Timing was done using the timemodule and time.clock()
计时是使用time模块完成的，并且time.clock()
回答by Trenton
maxb's answer gives good performance but can easily lead to loss of precision when coordinate values or the number of points are large. This can be mitigated with a simple coordinate shift:
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)
For example, a common geographic reference system is UTM, which might have (x,y) coordinates of (488685.984, 7133035.984). The product of those two values is 3485814708748.448. You can see that this single product is already at the edge of precision (it has the same number of decimal places as the inputs). Adding just a few of these products, let alone thousands, will result in loss of precision.
例如，常见的地理参考系统是 UTM，它的 (x,y) 坐标可能为(488685.984, 7133035.984)。这两个值的乘积是3485814708748.448。可以看到这个单品已经处于精度的边缘（它的小数位数与输入相同）。添加这些产品中的少数几个，更不用说数千个，都会导致精度损失。
A simple way to mitigate this is to shift the polygon from large positive coordinates to something closer to (0,0), for example by subtracting the centroid as in the code above. This helps in two ways:
缓解这种情况的一种简单方法是将多边形从大的正坐标移动到更接近 (0,0) 的位置，例如通过减去上面代码中的质心。这有两个方面的帮助：
- It eliminates a factor of
x.mean() * y.mean()from each product - It produces a mix of positive and negative values within each dot product, which will largely cancel.
- 它消除了
x.mean() * y.mean()每个产品中的一个因素 - 它在每个点积中产生正值和负值的混合，这将在很大程度上抵消。
The coordinate shift does not alter the total area, it just makes the calculation more numerically stable.
坐标偏移不会改变总面积，只会使计算在数值上更加稳定。
