I’m in the process of finishing up a web-app that I’ve written using Dash/Plotly and as someone who is not an expert in front end web development, I’ve really enjoyed how easy it is to build an interactive app without getting bogged down by all the javascript/html/css (though knowing some definitely helps!). But recently I’ve run into an issue while figuring how to plot various datasets as a background contour on a number of interactive maps: while Dash/Plotly is great at plotting interactive points on a map, there’s not a simple way to contour a dataset and use it as the base layer (there is the new scattermapbox class, but there is not yet a simple contouring method for maps like, for example, contourf in the basemap library).
But Dash/Plotly is very good at chloropleth (colored shapes) plots, so I figured that if I could convert a matplotlib.pyplot.contourf plot into a shapefile of polygons, I could then load that shapefile in Dash/Plotly as needed and plot it as a chloropleth plot, which seems easier to me than rendering a bitmap image of the contour plot in the background (the only other option I can think of)… and so I spent the morning experimenting with how to do this. There’s a few scattered stackoverflow questions that related to this, but I had to pull together a bunch of pieces to get this working, so I figured a short writeup here might be useful to others out there.
So with that said, here are they python libraries that you’ll need to run my test script:
- shapely (for polygons, used previously here, here and here)
- fiona (for writing a shapefile)
- descartes (for making some plots)
- matplotlib (for the initial contouring of the data)
- numpy (for making up some initial data and some array manipulation)
The full code is available in my learning_shapefiles repo and the new script is contourf_to_shp.py. And some useful stackoverflow answers that helped significantly: Converting Matplotlib contour objects to Shapely objects, matplotlib – extracting values from contour lines.
Creating data to contour
To start off, I wanted to work with a small test data set that would include multiple domains at the same contour level, so I created a lat/lon grid and then superimposed two 2d gaussian curves with different amplitudes and decay rates and plotted those up:
from shapely import geometry
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
import fiona
import os,json
from descartes.patch import PolygonPatch
# create some test data with multiple peaks
lon = np.linspace(0,45,100)
lat = np.linspace(-20,32,90)
long,latg=np.meshgrid(lon,lat)
C1=np.sqrt((long-5.)**2+(latg-25)**2)/30.
C2=np.sqrt((long-30.)**2+(latg-1)**2)/10.
m = 30*np.exp(-C1**2)+20.*np.exp(-C2**2)
# make the contourf plot, storing the resulting ContourSet in cs
plt.figure(figsize=[10,5])
plt.subplot(1,2,1)
Nlevels=10
cs = plt.contourf(lon,lat,m,Nlevels,cmap='gist_heat')
plt.title('contourf figure with Nlevels='+str(Nlevels))
In my actual use-case I’ll be loading a 2D dataset, so most of this will be replaced. But the important bit is that I store the object that results from the contourf call, that cs object contains all the data on the contours.
Storing the contour results
The first step in saving the contoured data is to create a lookup table for the contour levels. The cs object stores each contour in cs.collections and the levels in cs.levels:
# create lookup table for levels lvl_lookup = dict(zip(cs.collections, cs.levels))
Next, I loop over each contour of cs.collections, convert the contour’s shape into a shapely polygon and store it in a list of dictionaries that also saves that contour’s level value:
# loop over collections (and polygons in each collection), store in list for fiona
PolyList=[]
for col in cs.collections:
z=lvl_lookup[col] # the value of this level
for contour_path in col.get_paths():
# create the polygon for this level
for ncp,cp in enumerate(contour_path.to_polygons()):
lons = cp[:,0]
lats = cp[:,1]
new_shape = geometry.Polygon([(i[0], i[1]) for i in zip(lons,lats)])
if ncp == 0:
poly = new_shape # first shape
else:
poly = poly.difference(new_shape) # Remove the holes
PolyList.append({'poly':poly,'props':{'z': z}})
The list PolyList now holds a shapely polygon as well as the z-value (or level) of that polygon.
Writing the data
In order to use these polygons elsewhere, the easiest thing to do is simply save them within a shapefile. I used the Fiona library for writing, which entails writing the geometry (the polygons) and the properties for each polygon that I stored earlier:
# define ESRI schema, write each polygon to the file
outfi=os.path.join(outname,'shaped_contour.shp')
schema = {'geometry': 'Polygon','properties': {'z': 'float'}}
with fiona.collection(outfi, "w", "ESRI Shapefile", schema) as output:
for p in PolyList:
output.write({'properties': p['props'],
'geometry': geometry.mapping(p['poly'])})
In order to reconstruct the contour plot, I need to save the information about the different levels (to get the colormap correct). And while I could load the shapefile and loop over the shapes to calculate min/max values, it saves time to simply write that data to a different file. And rather than write it as metadata of some kind for the shapefile, I just drop it into a json file:
# save the levels and global min/max as a separate json for convenience
Lvls={'levels':cs.levels.tolist(),'min':m.min(),'max':m.max()}
with open(os.path.join(outname,'levels.json'), 'w') as fp:
json.dump(Lvls, fp)
Plotting the polygons
From there, the shapefile and levels data can be re-loaded and plotted to get the following:
On the left is the original contour plot, on the write is the plot of Polygon patches colored using the same colormap. And it’s an almost perfect replica… I’m not sure if it’s my eyes or screen, but the lightest colors seem just a hair different between the two. But that difference is even less noticeable as I use more contours, e.g. 1000:
The final bit of code is mostly self explanatory: looping over the shapes in the shapefile, then plotting colored polygon patches like in previous posts. The only tricky business is making sure the polygon facecolor matches the filled contour color. To do this, I load in the contour level json file that I saved off, and load the matplotlib colormap that I used in the initial call to contourf:
# read in levels, define colormap with open(os.path.join(outname,'levels.json')) as jfile: Lvls=json.load(jfile) levels=np.array(Lvls['levels']) cmap=plt.cm.gist_heat lv_range=[Lvls['min'],Lvls['max']]
And then for each polygon, I pull out the z-value that I saved and find the appropriate RGB value for that level with:
lv=shp['properties']['z'] # this shape's level
clr=cmap((lv - lv_range[0])/(lv_range[1]-lv_range[0]))
The bit of math with lv_range is ensure that the value is between 0 and 1. The clr variable is a tuple that can be given to the descartes PolygonPatch:
patch = PolygonPatch(poly, facecolor=clr, edgecolor=clr)
And that’s pretty much it. I still have some experimenting to do with how many contour levels I need for the datasets I’m working with, and I expect to run into some scaling issues (the full dataset fields I’m working with are 100s of MB), but at least I have a feasible approach!
Chris Havlin 