matplotlib+cartopy画WRF小时降水

风暴之灵

这几天想学学python，画了一个WRF小时降水的图，不同时次的结果存储在同一个PDF文件中。
总体感觉画图不如NCL方便，lambert投影的时候无法添加经纬度坐标，也可能是我刚学，没有弄明白？
代码如下，给初学者大家一起学习交流：

1. from netCDF4 import Dataset
   
2. import xarray as xr
   
3. import numpy as np
   
4. import matplotlib.pyplot as plt
   
5. import cartopy.crs as ccrs
   
6. from matplotlib.cm import get_cmap
   
7. import matplotlib.colors as colors
   
8. from cartopy.io.shapereader import Reader
   
9. import cartopy.feature as cfe
   
10. import cmaps
    
11. from matplotlib.backends.backend_pdf import PdfPages
    
12. from wrf import (get_cartopy, latlon_coords, to_np, cartopy_xlim, cartopy_ylim,
    
13.                  getvar, ALL_TIMES)
    
14.                  
    
15. f = Dataset('wrfout_d01')
    
16. # rainc = np.array(f.variables['RAINC'][:])
    
17. # rainnc = np.array(f.variables['RAINNC'][:])
    
18. 
    
19. rainc = getvar(f, 'RAINC', timeidx=ALL_TIMES)
    
20. rainnc = getvar(f, 'RAINNC', timeidx=ALL_TIMES)
    
21. 
    
22. # f = xr.open_dataset('wrfout_d01')
    
23. # rainc = f['RAINC']
    
24. # rainnc = f['RAINNC']
    
25. 
    
26. tot = rainc.values + rainnc.values
    
27. wrf_time = f.variables['Times']
    
28. 
    
29. # tot = rainc + rainnc
    
30. dim_tot = tot.shape
    
31. h1prep = tot[2:dim_tot[0]:2,:,:] - tot[0:dim_tot[0]-2:2,:,:]
    
32. 
    
33. # lats = f['XLAT'].values[0]
    
34. # lons = f['XLONG'].values[0]
    
35. 
    
36. lats = f.variables['XLAT'][0]
    
37. lons = f.variables['XLONG'][0]
    
38. 
    
39. # lats, lons = latlon_coords(rainc)
    
40. cart_proj = get_cartopy(rainc)
    
41. # print(cart_proj)
    
42. 
    
43. levels = [0.1,0.5,1.,2.,3.,4.,5.,6.,8.,10.,20.,40]
    
44. norm = colors.BoundaryNorm(boundaries=np.array(levels), ncolors=len(levels)-1)
    
45. 
    
46. shp_path = '/Users/james/Documents/GitHub/py_china_shp/'
    
47. reader = Reader(shp_path  +  'Province_9/Province_9.shp')
    
48. provinces = cfe.ShapelyFeature(reader.geometries(), ccrs.PlateCarree(), edgecolor='k', facecolor='none')
    
49. 
    
50. 
    
51. with PdfPages('multipage_pdf.pdf') as pdf:
    
52.     for i in range(3):
    
53.         fig = plt.figure(figsize=(14,8.5))
    
54.         ax = fig.add_subplot(111, projection=cart_proj)
    
55.         print(ax)   
    
56. 
    
57.         ax.add_feature(provinces, linewidth=0.6)
    
58.         ax.coastlines('50m', linewidth=0.8)
    
59.    
    
60.         pcm = ax.contourf(lons, lats, h1prep[i], levels, norm=norm,
    
61.                      transform=ccrs.PlateCarree(), extend='both',
    
62.                      cmap=cmaps.cmorph[1:])
    
63.         #pcm.cmap.set_over('white')
    
64. 
    
65.         pcm.cmap.set_under(cmaps.cmorph.colors[0])  
    
66. 
    
67.         # Set the map bounds
    
68.         # 两种方式设定画图范围
    
69.         # 方式1
    
70.         # xs, ys, _ = cart_proj.transform_points(ccrs.PlateCarree(),np.array([lons[0,0],lons[-1,-1]]),np.array([lats[0,0],lats[-1,-1]])).T
    
71.         # xlimits = xs.tolist()
    
72.         # ylimits = ys.tolist()
    
73.         # ax.set_xlim(xlimits)
    
74.         # ax.set_ylim(ylimits)
    
75.         # print(xlimits)
    
76.         # print(ylimits)
    
77.         # 方式2
    
78.         # ax.set_xlim(cartopy_xlim(rainc))
    
79.         # ax.set_ylim(cartopy_ylim(rainc))  
    
80. 
    
81.         ax.set_title("WRF Prepcipitation {}".format(str(wrf_time[i],'utf-8')))
    
82.         ax.gridlines(linewidth=0.5, color='gray', linestyle='--')   
    
83. 
    
84.         # plt.colorbar(pcm, ax=ax, ticks=levels, extendfrac='auto',
    
85.         #              aspect=20, extendrect=True, format='%3.1f', drawedges=True, shrink=.7)
    
86.         plt.colorbar(pcm, ax=ax, ticks=levels, extendfrac='auto',
    
87.                      aspect=18, format='%3.1f', shrink=.95)
    
88.         pdf.savefig()  # saves the current figure into a pdf page
    
89.         plt.close()
    
90. 
    
91.     # plt.suptitle('Figure title')
    
92.     # # shrink 调整colorbar大小   
    
93. 
    
94.     # plt.show()

复制代码

以下是画出来的图：

暗藏

lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()

火枪裸鞋

楼上正解！

风暴之灵

暗藏 发表于 2020-4-28 17:19
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()

谢谢，我之前用的cartopy0.17版本，更新到0.18版本以后就可以画这个经纬度坐标了。

蝶与鲸

我也更新至18.0，使用兰伯特投影，还是会报错诶。TypeError: Cannot label LambertConformal gridlines. Only PlateCarree gridlines are currently supported.

言墨

蝶与鲸 发表于 2020-5-9 18:04
我也更新至18.0，使用兰伯特投影，还是会报错诶。TypeError: Cannot label LambertConformal gridlines. On ...

我也是  不过我的报错是：
    raise RuntimeError('Cannot handle non-rectangular coordinate '
RuntimeError: Cannot handle non-rectangular coordinate systems.

风暴之灵

风暴之灵 发表于 2020-5-7 10:58
谢谢，我之前用的cartopy0.17版本，更新到0.18版本以后就可以画这个经纬度坐标了。

这个不清楚了，我更新了之后就不会报这个错误了。。

青山见我应如是

风暴之灵 发表于 2020-5-11 11:04
这个不清楚了，我更新了之后就不会报这个错误了。。

请问，能不能把添加坐标的代码发一下

hanf

青山见我应如是 发表于 2020-5-20 19:11
请问，能不能把添加坐标的代码发一下

lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)

青山见我应如是

hanf 发表于 2020-5-22 15:53
lon_formatter = cticker.LongitudeFormatter()
lat_formatter = cticker.LatitudeFormatter()
ax.xaxi ...

如果cticker.LongitudeFormatter()指的是from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter的话，我想问一下 ax.set_xticks(）语句该怎么写 ，只写lon_formatter = cticker.LongitudeFormatter()lat_formatter = cticker.LatitudeFormatter()ax.xaxis.set_major_formatter(lon_formatter)ax.yaxis.set_major_formatter(lat_formatter)不显示坐标啊，但是加上ax.set_yticks(np.linspace(box[2], box[3],5), crs=proj )就报错