按照maskout和帖子的plot脚本，我画出来的图是空白的，这是为什么？

Nozigaka46

按照maskout和帖子的plot脚本，https://bbs.06climate.com/forum.php?mod=viewthread&tid=93995&extra=&page=1我画出来的图是空白的，这是为什么？


import netCDF4 as nc      
import numpy as np
from wrf import (smooth2d, interplevel, getvar, to_np, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords)
import sys
sys.path.append(r'C:\Users\a652887227\Desktop\新建文件夹\代码\地图')
import maskout

# 读取数据
ncfile_s = nc.Dataset(r'D:\Meteorology Data\GFS Data\output\2021.08.21-2021.08.23\1\wrfout_d02_2021-08-21_00-00-00')
ncfile_e = nc.Dataset(r'D:\Meteorology Data\GFS Data\output\2021.08.21-2021.08.23\1\wrfout_d02_2021-08-22_00-00-00')
# 获取变量并做平滑处理
rainc_s = smooth2d(getvar(ncfile_s,"RAINC"),3)
rainnc_s = smooth2d(getvar(ncfile_s,"RAINNC"),3)
rainsh_s = smooth2d(getvar(ncfile_s,"RAINSH"),3)
rainc_e = smooth2d(getvar(ncfile_e,"RAINC"),3)
rainnc_e = smooth2d(getvar(ncfile_e,"RAINNC"),3)
rainsh_e =smooth2d(getvar(ncfile_e,"RAINSH"),3)

prei = (rainc_e + rainnc_e + rainsh_e) - (rainc_s - rainnc_s - rainsh_s)
#降雨是个累积量，而不是瞬时值，在wrfout中如果需要得到模拟起止日期内的降雨，需要将最后一个时刻的降雨量减去第一个时刻的降雨量。

# 获取关注对象的经纬点
lats, lons = latlon_coords(prei)

# 匹配兰伯特投影                                            
cart_proj = get_cartopy(rainnc_e)                       

# =========================================================================
# ==========================画图部分========================================
#import cartopy.feature as cfeature
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as ccrs
import cmaps
import matplotlib.colors as colors
import cartopy.io.shapereader as shpreader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
plt.rcParams['font.sans-serif']=['SimHei']  #显示中文
plt.rcParams['axes.unicode_minus']=False#用来正常显示负号

plotcrs = ccrs.PlateCarree()
fig = plt.figure(1, figsize=(16,9))
ax = plt.axes(projection=cart_proj)
proj =  cart_proj

zone=[97., 109., 25., 35.]
ax.set_extent(zone, ccrs.PlateCarree())            #设置出图的经纬度范围

# 设置底图特征
china = shpreader.Reader(r'C:\Users\a652887227\Desktop\新建文件夹\代码\地图\四川省2\四川省.shp').geometries()
ax.add_geometries(china, plotcrs, facecolor='none', edgecolor='black', zorder=1)


gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, x_inline=False, y_inline=False, linewidth=1.2 ,color='k', alpha=0.5, linestyle='--')
gl.top_labels = False      #关闭顶端标签
gl.right_labels = False    #关闭右侧标签
gl.rotate_labels = False   #禁止标签旋转
gl.xformatter = LONGITUDE_FORMATTER                                               #x轴设为经度格式
gl.yformatter = LATITUDE_FORMATTER                                                #y轴设为纬度格式


############################################################################################################################
clevs = [0, 25, 50, 100, 250, ]
cdict = ['#fff', '#63B7FF','#0000FE', '#FF00FC', '#850042']
my_cmap = colors.ListedColormap(cdict)
norm = colors.BoundaryNorm(boundaries=np.array(clevs), ncolors=len(clevs)-1)

cs1 = plt.contourf(to_np(lons), to_np(lats), to_np(prei),levels=clevs, norm=norm,transform=ccrs.PlateCarree(), extend='both',cmap=my_cmap)
cb = plt.colorbar(cs1, orientation='vertical', extendfrac='auto',pad=0.06, shrink=0.6)




cb.set_label('Precipitation (mm)')   
plt.title('GFS 24h forcast_Prei\n',fontsize=12,loc='left')
plt.title('8.21 08:00:00 到 8.22 08:00:00',fontsize=12,loc='right')

clip1 = maskout.shp2clip(cs1, ax, r'C:\Users\a652887227\Desktop\新建文件夹\代码\地图\四川省2\四川省.shp', 510400, proj=proj)
plt.show()