Python处理.nc数据可视化 指定区域三小时变温图 WRF模式

有错误欢迎提出指正！
两个月前大气环境模式上机作业，部分库的下载有点点麻烦（好像是cartopy），不能直接用anaconda下载。该代码里面有部分内容是画累积降水的，没有删除，但是对运行这个没有影响o(￣▽￣)ブ

import netCDF4 as nc
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pandas as pd
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import numpy as np
from matplotlib.patches import Rectangle
from cartopy.io.shapereader import Reader
import matplotlib

plt.rcParams['font.family'] = ['sans-serif']
plt.rcParams['font.sans-serif'] = ['SimHei']
#读取文件
file = r'H:\wrfout.nc'
dataset =nc.Dataset(file,'r') #读取nc文件
Rain=dataset.variables['RAINNC'][:,::-1] #需要先找出nc文件里面的变量名再赋值
T2=dataset.variables['T2'][:,::-1]
lon=dataset.variables['XLONG'][:]
lat=dataset.variables['XLAT'][:,::-1]
#print(dataset.variables['T2'])
lon_need=np.array(lon[0][115:130,117:132])#这个区域是我根据降水分布图找到的降水量特大区，然后自己指定的该区域
lat_need=np.array(lat[0][115:130,117:132])
fig=plt.figure(figsize=(12,10),dpi=550)
#绘制3小时气温
fig.suptitle('整个模拟区域逐3小时气温',fontsize=8)
fig.suptitle('指定区域逐3小时气温(UTC)',fontsize=8)

levels=np.arange(20,30,1)
for i in range(1,10):
    ax=fig.add_subplot(2,5,i,projection=ccrs.PlateCarree())

    a=12+3*(i-1)
    if a>24:
        ax.set_title('22日{}时'.format(a-24),fontsize=3,loc='left',y=0.8)
    else:
        ax.set_title('21日{}时'.format(a),fontsize=3,loc='left',y=0.8)

    #地理信息
    ax.add_feature(cfeature.COASTLINE,lw=0.1)
    ax.add_feature(cfeature.LAND.with_scale('50m'))
    shp_path_1=r'E:\WRF\区划\省.shp' 
    reader_1=Reader(shp_path_1)
    enshicity_1 = cfeature.ShapelyFeature(reader_1.geometries(), crs=ccrs.PlateCarree(), edgecolor='k', facecolor='none')
    ax.add_feature(enshicity_1,lw=0.1)
    ##经纬度范围
    #extent=[105,120,20,30]
    extent=[108,111,24,26]
    ax.set_extent(extent)

    #绘制经纬度
    gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.2, color='k', alpha=0.5, linestyle='--')
    gl.top_labels = False ##关闭上侧坐标显示
    gl.right_labels = False ##关闭右侧坐标显示
    if i !=1 and i!=6 :
        gl.left_labels=False

    gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式
    gl.yformatter = LATITUDE_FORMATTER
    gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1], 1))
    gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3], 1))
    gl.xlabel_style={'size':3}
    gl.ylabel_style={'size':3}
    T = T2-273.15
    b = ax.contourf(lon[0],lat[0],T[i*3-3],cmap='Spectral_r',levels=levels,extend='both')
    x = ax.plot(lon_need[0], lat_need[0], color='b', linewidth='0.5', linestyle='--')
    y = ax.plot(lon_need[0], lat_need[14], color='b', linewidth='0.5', linestyle='--')
    m = ax.plot(lon_need[:, 0], lat_need[:, 0], color='b', linewidth='0.5', linestyle='--')
    n = ax.plot(lon_need[:, 14], lat_need[:, 0], color='b', linewidth='0.5', linestyle='--')

position=fig.add_axes([0.3,0.175,0.4,0.02])
cb=plt.colorbar(b, cax=position,extend='both',orientation='horizontal')
font = {'family' : 'serif',
        'color'  : 'darkred',
        'weight' : 'normal',
        'size'   : 6,
        }
cb.ax.tick_params(labelsize=5)
plt.subplots_adjust(left=0.110, bottom=0.205, right=0.900, top=0.835, wspace=0.165, hspace=0)

plt.show()

欢迎各位同行的指正与讨论~