Converting Matplotlib contour objects to Shapely objects

Question

I would like to convert a Matplotlib contour collection to a Shapely geometry collection, in order to calculate the area contained within closed contour lines.

it seems that all the information is there in the Matplotlib contour collection to create Shapely polygon collection, based on these stackoverflow answers on extracting data from matplotlib contours: https://stackoverflow.com/questions/5666056/matplotlib-extracting-data-from-contour-lines and https://stackoverflow.com/questions/17051131/matplotlib-extracting-values-from-contour-lines/17056513#17056513

Has someone already done this?

Answer 1

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.

Answer 2

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.

Answer 3

There is also an issue with the answer from nbeuchat.

If the levels you pass to contourf intersect z, get_paths() will return some lines, which cannot be converted to polygons. to_polygons() returns None, shapely allows Polygons to be defined with empty features, and difference is supported for polygons with empty features, which leads to multiple counting, since the difference between an empty polygon and a non-empty polygon is the non-empty polygon.

This code rectifies the problem:

    contours = contourf(lons, lats, z, levels = levels)
    for i, collection in enumerate(contours.collections):
        for path in collection.get_paths():
            if path.to_polygons():
                for npoly, polypoints in enumerate(path.to_polygons()):
                    poly_lons = polypoints[:, 0]
                    poly_lats = polypoints[:, 1]
                    # be careful with the following---check to make sure
                    # your coordinate system expects lat first, lon second
                    poly_init = Polygon([coords for coords in \
                                         zip(poly_lats, poly_lons)])
                    if poly_init.is_valid:
                        poly_clean = poly_init
                    else:
                        poly_clean = poly_init.buffer(0.)
                    if npoly == 0:
                        poly = poly_clean
                    else:
                        poly = poly.difference(poly_clean)