Basemap contourf inside multiple shapes

Question

I am trying to plot interpolated weather data on a map only within the bounds of the towns contained in a shapefile. The following is the unclipped contourf over the Basemap with imported shapefile: Contourf overlaid on Basemap with Shapefile

I have tried clipping to contour collections by iterating through the contour collections like so:

m.readshapefile('data/grense', 'grense',zorder=10,linewidth=1, 
drawbounds=True)

patches   = []
for info, shape in zip(m.grense_info, m.grense):
   patches.append( Polygon(np.array(shape), linestyle=':', fill=False) )

for poly in patches:
   for collection in cs.collections:
      collection.set_clip_path(poly)

This obviously confines the contour to one polygon, i.e. one town, like this: Contourf clipped to one ploygon

Is it possible to create a collection of contour collections which I can then add by using ax.add_collection(new_contour_collection)? Something along the lines of:

for poly in patches:
   for collection in cs.collections:
     contour_collection.append(collection)
ax.add_collection(contour_collection)

Or can I create a single path from the Patchcollection and then clip each of the contour collections using collection.set_clip_patch(patches)?

If you only want visualization, plot your contourf at lower zorder. Then plot filled polygons with higher zorder to hide un-wanted parts of contourf. — swatchai
Very nice solution. Please consider posting the solution as an answer to your own question and accepting it. This way others see that your problem need no further attention (and it may score you some votes). — Thomas Kühn

Stef Marais Stef Marais · Accepted Answer · 2018-12-15T11:50:09

Following swatchai's suggestion and Thomas Kühn's previous answer here, I managed to solve my problem, seen here ,

by doing the following:

#Combine shapefile shape object (m.grense) with map edges
##limits of the map:
x0,x1 = ax.get_xlim()
y0,y1 = ax.get_ylim()
map_edges = np.array([[x0,y0],[x1,y0],[x1,y1],[x0,y1]])

polys = [map_edges] + m.grense

codes = [
    [Path.MOVETO] + [Path.LINETO for p in p[1:]]
    for p in polys
]
polys_lin = [v for p in polys for v in p]
codes_lin = [c for cs in codes for c in cs]
path = Path(polys_lin, codes_lin)

#Important - Set Zorder greater than Contour and less than Map 
borders
patch = PathPatch(path,facecolor='white', lw=0, zorder =2)

##masking the data:
ax.add_patch(patch)

Basemap contourf inside multiple shapes

1 Answers