How can I calculate the area within a contour in Python using the Matplotlib?

Question

I am trying to figure out a way to get the area inside a specific contour line?
I use matplotlib.pyplot to create my contours.
Does anyone have experience for this?

Thanks a lot.

Answer 1

From the collections attribute of the contour collection, which is returned by the contour function, you can get the paths describing each contour. The paths' vertices attributes then contain the ordered vertices of the contour.

Using the vertices you can approximate the contour integral 0.5*(x*dy-y*dx), which by application of Green's theorem gives you the area of the enclosed region.

However, the contours must be fully contained in the plot, because otherwise the contours are broken up into multiple, not necessarily connected paths and the method breaks down.

Here's the method used to compute the area enclosed of the radius function, i.e. r = (x^2 + y^2)^0.5, for 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()

Output:

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

Answer 2

A vectorized version of @spfrnd's answer to compute the area:

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)

Note that you may need to take the abs of the area because if the points along the contour are oriented in the opposite direction the result will be negative.

Answer 3

Obviously, the results for r=1,2,3 in @spfrnd's answer are far off the exact areas for circles according to A(r) = pi r^2, even for a rather dense grid. The reason the above code doesn't work properly is that the returned vertices are incomplete due to the clabels generated by

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

and consequently, the shoelace algorithm calculates the wrong area.

You can verify this easily by plotting the returned vertices next to the contour plot (best use a larger delta or keep only every Nth point via the :: operator)

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

For a visualization, see this picture.

Removing the clabels leads to quite accurate results, even for rather coarse grids. Setting inline=False in the clabel command does the job as well.

Answer 4

The question regards a "specific contour line" and spfrnd toghether with PandaScience's response is a good solution.

However, the method will not give a correct value if there are multiple islands for the same contour level. The code below handles multiple islands:

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