pyton画青藏高原及要素散点图

chiqu296

今天周末参考（http://bbs.06climate.com/forum.php?mod=viewthread&tid=38961）用python画了一下青藏高原及西藏地区降水量的散点图，分享之大家可以参考！cmaps（http://bbs.06climate.com/forum.php?mod=viewthread&tid=43521）

1. import cmaps
   
2. import pandas as pd
   
3. import numpy as np
   
4. from netCDF4 import Dataset
   
5. import matplotlib.pyplot as plt
   
6. from mpl_toolkits.basemap import Basemap
   
7. 
   
8. 
   
9. plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文
   
10. plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
    
11. 
    
12. 
    
13. # 打开cvs数据
    
14. data_all = pd.read_csv('rain.dat',header=None, names=['站点','站号','lons','lats','降水量','气温'])
    
15. # 打开高程数据
    
16. deg_data = Dataset('elev.0.25-deg.nc','r')
    
17. data = deg_data.variables['data'][0,:,:]
    
18. # 缺测值maske掉
    
19. data = np.ma.masked_values(data,32767)
    
20. lon = deg_data.variables['lon'][:]
    
21. lat = deg_data.variables['lat'][:]
    
22. # data_units = deg_data['data'].units
    
23. lon, lat = np.meshgrid(lon, lat)
    
24. 
    
25. 
    
26. # 画图
    
27. fig = plt.figure(figsize=(16,16))
    
28. # 设置全局字体等
    
29. plt.rc('font',size=15,weight='bold')
    
30. plt.subplot(121)
    
31. m = Basemap(projection='cyl',llcrnrlat=25,llcrnrlon=70,urcrnrlat=45,urcrnrlon=110)
    
32. # 加载中国地图文件
    
33. m.readshapefile('china_shp/cnhimap','cnhimap.shp',linewidth=1,color ='k')
    
34. # 加载青藏高原边界文件
    
35. m.readshapefile('DBATP/DBATP_Line','DBATP_Line.shp',linewidth = 2,color='r')
    
36. x,y = m(lon,lat)
    
37. levels = np.linspace(3000,np.max(data),50)
    
38. cf = m.contourf(x,y,data,levels=levels,cmap=cmaps.MPL_gist_earth)
    
39. cbar = m.colorbar(cf,location='bottom',size='5%',pad='15%',
    
40.     ticks=np.linspace(3000,np.max(data),10),label='m')
    
41. lon_num = np.arange(70,111,5)
    
42. lon_name = ['70°','75°','80°','85°','90°','95°','100°','105°','110°E']
    
43. lat_num =  np.arange(25,46,5)
    
44. lat_name = ['25°','30°','35°','40°','45°N']
    
45. plt.yticks(lat_num,lat_name)
    
46. plt.xticks(lon_num,lon_name)
    
47. plt.title('(a)')
    
48. plt.subplot(122)
    
49. m = Basemap(projection='cyl',llcrnrlat=25,llcrnrlon=70,urcrnrlat=45,urcrnrlon=110)
    
50. # 加载中国地图文件
    
51. m.readshapefile('china_shp/cnhimap','cnhimap.shp',linewidth=1,color ='black')
    
52. # 加载青藏高原边界文件
    
53. m.readshapefile('DBATP/DBATP_Line','DBATP_Line.shp',linewidth = 2,color='red')
    
54. x,y = m(data_all['lons'],data_all['lats'])
    
55. st = m.scatter(x,y,c=data_all['降水量'],s=data_all['降水量'],marker='o',cmap=cmaps.MPL_BuGn)
    
56. cbar = m.colorbar(st,location='bottom',size='5%',pad='15%',
    
57.                   ticks=np.linspace(np.min(data_all['降水量']),np.max(data_all['降水量']),10),
    
58.                   label='mm',format='%.1f',extend='max')
    
59. 
    
60. lon_num = np.arange(70,111,5)
    
61. lon_name = ['70°','75°','80°','85°','90°','95°','100°','105°','110°E']
    
62. lat_num =  np.arange(25,46,5)
    
63. lat_name = ['25°','30°','35°','40°','45°N']
    
64. plt.yticks(lat_num,lat_name)
    
65. plt.xticks(lon_num,lon_name)
    
66. plt.title('(b)')
    
67. plt.savefig('map_data_test.png', bbox_inches='tight',dpi=300)

复制代码

xigua

楼主厉害，果断收藏

weinihou

谢谢分享

（吐槽下帖子题目的拼写错误）

ybmy001

谢谢分享

（吐槽下帖子题目的拼写错误）

ybmy001

谢谢分享

（吐槽下帖子题目的拼写错误）

ybmy001

谢谢分享

（吐槽下帖子题目的拼写错误）

ybmy001

谢谢分享

（吐槽下帖子题目的拼写错误）

ybmy001

谢谢分享

（吐槽下帖子题目的拼写错误）

ybmy001

谢谢分享

（吐槽下帖子题目的拼写错误）

huangya2009

感谢分享，学习一下