5

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?

9

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:

enter image description here

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)

enter image description here

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

enter image description here

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.

  • I think there is a step missing here. type(poly) returns matplotlib.patches.Polygon, which doesn't have an area method or property. How do I convert this to a shapely object? – Rich Signell Jun 5 '14 at 14:25
  • Polygon([(i[0], i[1]) for i in zip(x,y)]) is a shapely object, look at The Shapely User Manual – gene Jun 5 '14 at 20:54
  • Grr... was using ipython notebook --pylab and was getting the wrong Polygon constructor (from Matplotlib instead of Shapely) – Rich Signell Jun 5 '14 at 21:09
  • Ipython notebook here: nbviewer.ipython.org/gist/rsignell-usgs/dd0c63ea8755d6317db2 – Rich Signell Jun 5 '14 at 21:18
4

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.