Converting Matplotlib contour objects to Shapely objects

If I use your first example matplotlib - extracting data from contour lines

import matplotlib.pyplot as plt
x = [1,2,3,4]
y = [1,2,3,4]
m = [[15,14,13,12],[14,12,10,8],[13,10,7,4],[12,8,4,0]]
cs = plt.contour(x,y,m)

The result is:

The number of elements (lines) is given by:

len(cs.collection)
7

and the result you want is the area of one of the polygons (with contourf(): 7 polygons)

In fact, the xy list determined by:

p = cs.collections[0].get_paths()[0]
v = p.vertices
x = v[:,0]
y = v[:,1]

are the coordinates of the exterior LinearRings of the coloured Polygons. So

from shapely.geometry import polygon
for i in range(len(cs.collections)):
    p = cs.collections[i].get_paths()[0]
    v = p.vertices
    x = v[:,0]
    y = v[:,1]
    poly = Polygon([(i[0], i[1]) for i in zip(x,y)])
    print i, poly
    0 POLYGON ((4 3.5, 4 4, 3.5 4, 4 3.5))
    1 POLYGON ((4 3, 4 3, 4 3.5, 3.5 4, 3 4, 3 4, 3 4, 4 3, 4 3))
    2 POLYGON ((4 2.5, 4 3, 4 3, 3 4, 3 4, 2.5 4, 3 3.333333333333333, 3.333333333333333 3, 4 2.5))
    3 POLYGON ((4 2, 4 2, 4 2.5, 3.333333333333333 3, 3 3.333333333333333, 2.5 4, 2 4, 2 4, 2 4, 2.666666666666667 3, 3 2.666666666666667, 4 2, 4 2))
    4 POLYGON ((3 2, 4 1.5, 4 2, 4 2, 3 2.666666666666667, 2.666666666666667 3, 2 4, 2 4, 1.5 4, 2 3, 2 3, 3 2, 3 2))
    5 POLYGON ((4 1, 4 1, 4 1.5, 3 2, 3 2, 2 3, 2 3, 1.5 4, 1 4, 1 4, 1.333333333333333 3, 2 2, 2 2, 3 1.333333333333333, 4 1))
    6 POLYGON ((2 1, 2 1, 3 1, 4 1, 3 1.333333333333333, 2 2, 2 2, 1.333333333333333 3, 1 4, 1 3, 1 2, 1 2, 2 1))
    7 POLYGON ((1 1, 2 1, 1 2, 1 1))

Plot of the Polygon 4

and the result is given by poly.area

But there are other solutions as in matplotlib - users: pyplot: Extract contourset without plotting or stackoverflow : Python: find contour lines from matplotlib.pyplot.contour() with the undocumented module matplotlib._cntr without plotting anything.

Issue with accepted answer:

To complete the accepted answer, one should note that the method will fail if either of these is true:

1. There is more than one polygon for a given level
2. There are "holes" in the polygon (in this case, the accepted answer would work but would create an invalid Polygon which can be problematic down the line)

Code:

The following code would fix both problems at once:

from shapely import geometry
for col in cs.collections:
    # Loop through all polygons that have the same intensity level
    for contour_path in col.get_paths(): 
        # Create the polygon for this intensity level
        # The first polygon in the path is the main one, the following ones are "holes"
        for ncp,cp in enumerate(contour_path.to_polygons()):
            x = cp[:,0]
            y = cp[:,1]
            new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(x,y)])
            if ncp == 0:
                poly = new_shape
            else:
                # Remove the holes if there are any
                poly = poly.difference(new_shape)
                # Can also be left out if you want to include all rings

        # do something with polygon
        print poly 

Explanations:

1. If more than one polygon exist with the same intensity level, .get_paths() will contain more than one item. Therefore, looping on .get_paths() allows one not to miss any polygon.
2. If there are holes, the vertices property returns all of the points in the polygons, regardless if there are on the exterior or interior. Therefore, one should create a polygon with the exterior and remove all the polygons inside. Using .to_polygons() allows to get all the polygons (exterior and interior), the first one being the exterior one. With the difference function, you can remove all the holes.