python 站点资料插值画图及白化

chiqu296

在学习python气象数据可视化过程中，在论坛里学习到了很多，今天尝试用python将站点资料插值画图，效果还可以，做个记录！用到了又是那隻貓的源码python中使用NCL的colormap以及平流层的萝卜的Python完美白化。在此感谢他们的分享！

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

复制代码

黯隐琉璃

跪求楼主大大回复 ！可不可以共享一下您的maskout文件。我按照“平流层的萝卜”分享的包里面的程序maskout.py以及country1.shp就可以画出来。但是一旦换成bou2_4.shp或者其他shp文件，都失败。maskout里面的参数我也改了，但是还是不行。所以，可不可以看一下您的maskout.py是怎么改的？非常感谢！！！！！

lifuxi

好东西，但是好像是因为我的电脑装了2.7和3.6版本的python，所以我一直装不了basemap，想问问有什么解决办法

charleszhou

谢谢分享！学习到了

charleszhou

学习到了很多，很感谢

fei@fei

你好，我降雨数据是excel用你的代码报MemoryError得错，excel里有9万6千多条数据，请教下怎么处理啊

当未来

楼主我看你用的线性插值，有没有试过cressman的

sxpeng3677

画出来图 ，很漂亮啊

yandongna

确实很实用，楼主，有个小问题想请教您，就是在import maskout, import camps总是找不到相关文件，楼主有办法解决吗？

zdqpla

不收费真好