如何在Python中使用Matplotlib计算轮廓内的面积？

Question

我想找个办法来计算一个特定轮廓线内部的面积。
我使用 matplotlib.pyplot 来绘制我的轮廓图。
有没有人有这方面的经验呢？

非常感谢！

Answer 1

这个问题是关于“特定的轮廓线”，还有PandaScience的回答是个不错的解决方案。

不过，如果同一个轮廓线有多个岛屿，这个方法就不能给出正确的值。下面的代码可以处理多个岛屿的情况：

for i in range(len(levels)):
    contour = cs.collections[i]
    vs = contour.get_paths()
    area=0
    for island in vs:
        ai=np.abs(area(island.vertices))
        area+=ai

Answer 2

显然，在@spfrnd的回答中，当r=1,2,3时，得到的结果和根据公式A(r) = pi r^2计算的圆的实际面积差得很远，即使是在一个比较密集的网格上。上面的代码之所以不能正常工作，是因为返回的顶点不完整，这个问题是由

plt.clabel(cs, inline=1, fontsize=10)

生成的clabels导致的，因此，鞋带算法计算出的面积也是错误的。

你可以通过将返回的顶点和轮廓图一起绘制来轻松验证这一点（最好使用更大的增量，或者通过::操作符只保留每N个点）。

N = 10
vs = cs.collections[0].get_paths()[0].vertices                                   
plt.plot(vs[::N, 0], vs[::N, 1], marker="x", alpha=0.5)

想要可视化效果，可以查看这张图片。

去掉clabels后，结果会变得相当准确，即使是在比较粗糙的网格上。将inline=False放在clabel命令中也能解决这个问题。

Answer 3

这是对@spfrnd的答案的一个向量化版本，用来计算面积：

x=contour.vertices[:,0]
y=contour.vertices[:,1]
area=0.5*np.sum(y[:-1]*np.diff(x) - x[:-1]*np.diff(y))
area=np.abs(area)

需要注意的是，你可能需要对面积取绝对值，因为如果轮廓上的点方向相反，计算出来的结果会是负数。

Answer 4

从轮廓集合的 collections 属性中，你可以获取到每个轮廓的路径，这些路径是通过 contour 函数返回的。路径的 vertices 属性里包含了轮廓的有序顶点。

利用这些顶点，你可以近似计算轮廓的积分 0.5*(x*dy-y*dx)，根据格林定理，这可以给你封闭区域的面积。

不过，轮廓必须完全在图表内，否则轮廓会被分成多个不一定相连的路径，这样方法就不管用了。

下面是用来计算半径函数封闭区域面积的方法，也就是 r = (x^2 + y^2)^0.5，对于 r=1.0, r=2.0, r=3.0 的情况。

import numpy as np
import matplotlib.pylab as plt

# Use Green's theorem to compute the area
# enclosed by the given contour.
def area(vs):
    a = 0
    x0,y0 = vs[0]
    for [x1,y1] in vs[1:]:
        dx = x1-x0
        dy = y1-y0
        a += 0.5*(y0*dx - x0*dy)
        x0 = x1
        y0 = y1
    return a

# Generate some test data.
delta = 0.01
x = np.arange(-3.1, 3.1, delta)
y = np.arange(-3.1, 3.1, delta)
X, Y = np.meshgrid(x, y)
r = np.sqrt(X**2 + Y**2)

# Plot the data
levels = [1.0,2.0,3.0]
cs = plt.contour(X,Y,r,levels=levels)
plt.clabel(cs, inline=1, fontsize=10)

# Get one of the contours from the plot.
for i in range(len(levels)):
    contour = cs.collections[i]
    vs = contour.get_paths()[0].vertices
    # Compute area enclosed by vertices.
    a = area(vs)
    print "r = " + str(levels[i]) + ": a =" + str(a)

plt.show()

输出：

r = 1.0: a = 2.83566351207
r = 2.0: a = 11.9922190971
r = 3.0: a = 27.3977413253

如何在Python中使用Matplotlib计算轮廓内的面积？

4 个回答

撰写回答