# 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? Jun 5, 2014 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, 2014 at 20:54
• Grr... was using `ipython notebook --pylab` and was getting the wrong `Polygon` constructor (from Matplotlib instead of Shapely) Jun 5, 2014 at 21:09
• Ipython notebook here: nbviewer.ipython.org/gist/rsignell-usgs/dd0c63ea8755d6317db2 Jun 5, 2014 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.

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)
``````