Re: [Matplotlib-users] contour plot confined by shapefile border

[Benjamin R. <ben...@ou...>]

On Sunday, June 19, 2011, Lukmanul Hakim <plg...@ya...> wrote:
> Hello,
>
> I would like to ask for some hints
> and help. I am currently trying to plot a "pseudo contour" over a Basemap. This contour is confined by borders obtained from a shapefile. I can generate the contour, I can retrieve the shapefile and put them on a basemap. What I have not been able to do is to limit the contour plot such that it only fills the area defined by a shapefile (in my case it is bordered by Province Lampung). You can use different shapefile for example, as the point of my question is how to limit contour fill by a shapefile.
>
> Shapefile can be downloaded from http://www.gadm.org/data/shp/IDN_adm.zip
>
> The generated figure can be found here:
> https://lh3.googleusercontent.com/-hMlwe6MAjdc/Tf6HlxYxn1I/AAAAAAAAAEY/exdzSv30ZL4/s640/to_matplotlib_userlists.png
>
> Thanks for the help!
> --
> Lukmanul Hakim
>
> Here's the code:
> ----------------
> from pylab import cm
> import matplotlib.pyplot as plt
> import numpy as np
> from mpl_toolkits.basemap import Basemap
>
> #Define area about Province Lampung
> ulat, llat, ulon, llon = -2.5, -6.5, 107, 103
> m = Basemap(projection='merc', lon_0=0.0, llcrnrlon=llon,
>                      llcrnrlat=llat, urcrnrlon=ulon, urcrnrlat=ulat,
>                      resolution='i')
>
>
> #Read shapefile of Province Lampung
> s = m.readshapefile('E:/Works/UNILA/Research/IDN_adm/LampungMap/LampungMap', 'lampung')
>
> #Define area for contour plot
> llon1,ulon1 = 103,106
> llat1,ulat1 = -6,-3.5
>
> #Generate random data
> nx,ny=5,5
> data2 = np.random.sample((ny,nx))
> x = np.linspace(llon1, ulon1, nx)
> y = np.linspace(llat1, ulat1, ny)
> X,Y = np.meshgrid(x,y)
> px,py = m(X,Y)
> m.contourf(px, py, data2, cmap=cm.jet)
> m.drawcoastlines()#(linewidth=0.5)
> m.drawparallels(np.arange(-6.5,-2.5,1.),labels=[1,0,0,0],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) # draw parallels
> m.drawmeridians(np.arange(102.,108.,1.),labels=[0,0,0,1],color='black',dashes=[1,0],labelstyle='+/-',linewidth=0.2) # draw meridians
>
> plt.show()

Not exactly sure how to do this, but if you can get a true/false mask
of the region in the same shape as the input data, then you can have a
masked array that goes into contourf.  Any area where the mask was
true will be blanked.

The hard part, though is getting the mask.

I hope that helps!
Ben Root