How to use cropped path for Baseemap polygon

Question

I want to use imshow (for example) to display some data inside the borders of the country (as an example, I chose the USA). The simple example below illustrates what I want:

import numpy as np import matplotlib.pyplot as plt from matplotlib.patches import RegularPolygon data = np.arange(100).reshape(10, 10) fig = plt.figure() ax = fig.add_subplot(111) im = ax.imshow(data) poly = RegularPolygon([ 0.5, 0.5], 6, 0.4, fc='none', ec='k', transform=ax.transAxes) im.set_clip_path(poly) ax.add_patch(poly) ax.axis('off') plt.show() 

Result:

Now I want to do this, but instead of a simple polygon, I want to use a complex US shape. I have created some examples of data contained in the "Z" array, as can be seen from the code below. It is these data that I want to display using the colourmap, but only within the borders of the continental United States.

So far I have tried the following. I get the form file from here , contained in the file "nationp010g.shp.tar.gz", and I use the Basemap module in python to build the USA. Please note that this is the only method I found that gives me the opportunity to get the polygon area I need. If there are alternative methods, I would be interested in them too. Then I create a polygon called "mainpoly", which is almost a polygon that I want to color blue:

Please note that only one body is colored, all other disjoint polygons remain white:

So, the blue colored area is almost what I want, note that there are undesirable borders around Canada because the border actually passes through some lakes, but this is a small problem. The real problem is why my imshow data is not displayed in the USA? Comparing my first and second code examples, I don’t understand why I do not get the clipped imshow in my second example, as I do in the first. Any help would be appreciated in understanding what I am missing.

 import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap as Basemap from matplotlib.patches import Polygon # Lambert Conformal map of lower 48 states. m = Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49, projection='lcc',lat_1=33,lat_2=45,lon_0=-95) shp_info = m.readshapefile('nationp010g/nationp010g', 'borders', drawbounds=True) # draw country boundaries. for nshape,seg in enumerate(m.borders): if nshape == 1873: #This nshape denotes the large continental body of the USA, which we want mainseg = seg mainpoly = Polygon(mainseg,facecolor='blue',edgecolor='k') nx, ny = 10, 10 lons, lats = m.makegrid(nx, ny) # get lat/lons of ny by nx evenly space grid. x, y = m(lons, lats) # compute map proj coordinates. Z = np.zeros((nx,ny)) Z[:] = np.NAN for i in np.arange(len(x)): for j in np.arange(len(y)): Z[i,j] = x[0,i] ax = plt.gca() im = ax.imshow(Z, cmap = plt.get_cmap('coolwarm') ) im.set_clip_path(mainpoly) ax.add_patch(mainpoly) plt.show() 

Update

I understand that the line

 ax.add_patch(mainpoly) 

doesn't even add a polygon shape to the plot. I do not use it correctly? As far as I know, mainpoly was correctly calculated using the Polygon () method. I checked that the coordinate inputs are reasonable:

 plt.plot(mainseg[:,0], mainseg[:,1] ,'.') 

which gives

Answer 1

I have been discussing this issue for so long.
And I found that the NCL language has a function of masking data outside the border.
Here is an example:

http://i5.tietuku.com/bdb1a6c007b82645.png

The contour graph is shown only on the border of China. Click here for the code.

I know that python has a package called PyNCL that supports all NCL code in the Python framework.
But I really want to build such a figure using a basemap. If you find out, post it online. I will study for the first time.

Thanks!

Add 2016-01-16

In a way, I figured it out.
This is my idea and code, and it inspired me to this question that I asked today.

My method:
1. Make a shapefile of the region of interest (for example, US) in shapely.polygon.
2. Check each point of the value to / from the polygon. 3. If the point value is outside the study area, mask it as np.nan

Introduction * polygon xxx was a city in China in the ESRI shapefile format. * fiona, A beautiful package was used here.

 # generate the shapely.polygon shape = fiona.open("xxx.shp") pol = shape.next() geom = shape(pol['geometry']) poly_data = pol["geometry"]["coordinates"][0] poly = Polygon(poly_data) 

It shows how:

http://i4.tietuku.com/2012307faec02634.png

 ### test the value point ### generate the grid network which represented by the grid midpoints. lon_med = np.linspace((xi[0:2].mean()),(xi[-2:].mean()),len(x_grid)) lat_med = np.linspace((yi[0:2].mean()),(yi[-2:].mean()),len(y_grid)) value_test_mean = dsu.mean(axis = 0) value_mask = np.zeros(len(lon_med)*len(lat_med)).reshape(len(lat_med),len(lon_med)) for i in range(0,len(lat_med),1): for j in range(0,len(lon_med),1): points = np.array([lon_med[j],lat_med[i]]) mask = np.array([poly.contains(Point(points[0], points[1]))]) if mask == False: value_mask[i,j] = np.nan if mask == True: value_mask[i,j] = value_test_mean[i,j] # Mask the np.nan value Z_mask = np.ma.masked_where(np.isnan(so2_mask),so2_mask) # plot! fig=plt.figure(figsize=(6,4)) ax=plt.subplot() map = Basemap(llcrnrlon=x_map1,llcrnrlat=y_map1,urcrnrlon=x_map2,urcrnrlat=y_map2) map.drawparallels(np.arange(y_map1+0.1035,y_map2,0.2),labels= [1,0,0,1],size=14,linewidth=0,color= '#FFFFFF') lon_grid = np.linspace(x_map1,x_map2,len(x_grid)) lat_grid = np.linspace(y_map1,y_map2,len(y_grid)) xx,yy = np.meshgrid(lon_grid,lat_grid) pcol =plt.pcolor(xx,yy,Z_mask,cmap = plt.cm.Spectral_r ,alpha =0.75,zorder =2) 

result

http://i4.tietuku.com/c6620c5b6730a5f0.png

http://i4.tietuku.com/a22ad484fee627b9.png

source result

http://i4.tietuku.com/011584fbc36222c9.png