基于python的气象区域站画图，考虑到区域站站点密集，使用线性差值比较卡顿，而且结果为变态矩阵（from scipy.interpolate import Rbf），建议使用（from scipy.interpolate import griddata）
1、使用的降水资料，基于cimiss的区域站资料，时间为上篇帖子所写的玛利亚台风时间。
直接上代码：
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.interpolate import Rbf
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as colors
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] #SimHei是黑体的意思
data = pd.read_csv('./2018.csv', header=0)
data_new=data.replace(to_replace=data["PRE_Time_0808"].max(),value=0)
temp=data_new[data_new["Day"]==5]
temp["PRE_Time_0808"].max()
rain_data = temp['PRE_Time_0808'].values
lon = temp['Lon'].values
lat = temp['Lat'].values
#线性差值
#numcols, numrows = 100,100
#olon = np.linspace(lon.min()-0.9,lon.max()+0.8, numcols)
#olat = np.linspace(lat.min()-0.5,lat.max()+0.5, numrows)
#olon,olat = np.meshgrid(olon,olat)
#func = Rbf(lon, lat, rain_data,function='linear')
#rain_data_new = func(olon, olat)
#https://my.oschina.net/u/2306127/blog/600628，
#linear' Linear interpolation (default) 双线性插值
#'cubic' Cubic interpolation 双三次插值
#'natural' Natural neighbor interpolation 自然邻近插值
#'nearest' Nearest neighbor interpolation 最近邻近插值
from scipy.interpolate import griddata
points=np.dstack((lon,lat))
points.resize((len(lat),2))
grid_x, grid_y = np.mgrid[lon.min()-0.9:lon.max()+0.8: 50j,lat.min()-0.9:lat.max()+0.8:50j]
rain_data_new= griddata(points, rain_data, (grid_x, grid_y), method='nearest')
rain_data_new
fig = plt.figure(figsize=(16,9))
plt.rc('font',size=15,weight='bold')
ax = fig.add_subplot(111)
lon_leftup=lon.min()-0.9;lat_leftup=lat.max()+1
lon_rightdown=lon.max()+0.8;lat_rightdown=lat.min()-1
res=0.1
lon=np.arange(lon_leftup,lon_rightdown,res)
lat=np.arange(lat_rightdown,lat_leftup,res)
m = Basemap(projection='cyl', llcrnrlat=lat_rightdown, urcrnrlat=lat_leftup, llcrnrlon=lon_leftup, urcrnrlon=lon_rightdown, resolution='l')
parallels = np.arange(20,90,10) #纬线
m.drawparallels(parallels,labels=[True,False,False,False],linewidth=0.2,dashes=[1,4])
meridians = np.arange(0,200,10) #经线
m.drawmeridians(meridians,labels=[False,False,False,True],linewidth=0.2,dashes=[1,4])
m.readshapefile('shp/cnhimap', 'cnhimap')
m.readshapefile('shp/ne_50m_coastline', 'ne_50m_coastline')
x,y = m(grid_x,grid_y)
xx,yy = m(lon,lat)
cdict = [(151 / 255, 250 / 255, 151 / 255), (49/ 255, 204 / 255, 49/ 255), (126/ 255, 191 / 255, 237 / 255)
, (0 / 255, 0 / 255, 255 / 255), (237/ 255, 0 / 255, 237 / 255)]
my_cmap = colors.ListedColormap(cdict,'pre3h')
my_cmap.set_under('w')
my_cmap.set_over((135/ 255, 25 / 255, 25 / 255))
lev=np.array([0.1,10,25,50,100,250])
norm3 = mpl.colors.BoundaryNorm(lev, my_cmap.N)
cf = m.contourf(x,y,rain_data_new, cmap=my_cmap,levels=lev,norm=norm3,extend='both')
cbar = m.colorbar(cf,location='right',format='%d',size=0.3,label='mm')
#cb=fig.colorbar(cf,ax=ax,pad=0.07,shrink=0.7,aspect=25,orientation='horizontal')
#for i in range(0,len(xx)):
# plt.text(xx[i],yy[i],data['站名'][i],va='center',fontsize=15)
# print (data['站名'][i])
#白化
import shapefile
from matplotlib.path import Path
from matplotlib.patches import PathPatch
def shp2clip(originfig,ax,m,shpfile,region):
sf = shapefile.Reader(shpfile)
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
#ID_0为第一维度，NAME_1为第5维度
if shape_rec.record[4] in region: ####这里需要找到和region匹配的唯一标识符，record[]中必有一项是对应的。
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
vertices.append(m(pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform=ax.transData)
for contour in originfig.collections:
contour.set_clip_path(clip)
return clip
#clip = shp2clip(cf,ax,m,'D:/shp_map/dongting7',[0])
#clip = shp2clip(cf,ax,m,'D:/shp_map/dongting','Hubei,Hunan,Hunan,Hunan')
#clip = shp2clip(cf,ax,m,'D:/shp_map/湖南省/hunannew','Hunan')
clip = shp2clip(cf,ax,m,'D:/shp_map/湖南省/chinanew','Hunan')
plt.savefig('test1.png', bbox_inches='tight',dpi=300)
plt.show()
白化要注意修改 if shape_rec.record[4] in region中列索引，不同的shp文件列索引不一样，建议使用arcgis打开shp看下属性，我的shp文件如下
白化说明2.png
这样便可以任意挑选你需要白化的省份了，大功告成。
test1.png
