import cmaps
import maskout
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import Rbf
from mpl_toolkits.basemap import Basemap
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
data = pd.read_csv('../rain.dat', header=None,
names=['站名','站号','lon','lat','降水','气温'] )
# 插值
lon = data['lon']
lat = data['lat']
rain_data = data['降水']
olon = np.linspace(78,100,88)
olat = np.linspace(26,38,88)
olon,olat = np.meshgrid(olon,olat)
# 插值处理
func = Rbf(lon, lat, rain_data,function='linear')
rain_data_new = func(olon, olat)
# 画图
fig = plt.figure(figsize=(16,9))
plt.rc('font',size=15,weight='bold')
ax = fig.add_subplot(111)
m = Basemap(projection='cyl',llcrnrlat=26,llcrnrlon=78,urcrnrlat=38,urcrnrlon=100)
m.readshapefile('../tibet_shp/xizang_all','xizang_all.shp', linewidth=1, color='k')
m.readshapefile('../tibet_shp/river_1','river_1.shp',linewidth=1,color ='b')
m.readshapefile('../tibet_shp/river_2','river_2.shp',linewidth=0.8,color ='b')
m.readshapefile('../tibet_shp/river_3','river_3.shp',linewidth=0.6,color ='b')
m.readshapefile('../tibet_shp/xzlake','xzlake.shp',linewidth=1,color ='b')
x,y = m(olon,olat)
xx,yy = m(lon,lat)
levels = np.linspace(0,np.max(rain_data_new),50)
cf = m.contourf(x,y,rain_data_new, levels=levels, cmap=cmaps.CBR_wet)
cbar = m.colorbar(cf,location='right',format='%d',size=0.3,
ticks=np.linspace(0,np.max(rain_data_new),10),label='毫米')
st = m.scatter(xx-0.1,yy,c='k',s=10,marker='o')
for i in range(0,len(xx)):
plt.text(xx[i],yy[i],data['站名'][i],va='center',fontsize=10)
lon_num = np.arange(78,101,2)
lon_label = ['78°','80°','82°','84°','86°','88°','90°','92°','94°','96°','98°','100°E']
lat_num = np.arange(26,39,2)
lat_label = ['26°','28°','30°','32°','34°','36°','38°N']
plt.yticks(lat_num,lat_label)
plt.xticks(lon_num,lon_label)
plt.title('测试图')
# 白化
clip = maskout.shp2clip(cf,ax,m,'shapefile/bou2_4p',[540000])
plt.savefig('test.png', bbox_inches='tight',dpi=300)