# Converting Matplotlib contour objects to Shapely objects

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

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.

• 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

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.