python中纬向平均、经向平均、水平流场、水汽柱浓度计算、坐标刻度朝向、label上下角标示例

import pandas as pd
import numpy as np
import xarray as xr
from wrf import to_np,interpz3d,destagger
import glob 
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cmaps
import warnings

warnings.filterwarnings('ignore')

file_names = glob.glob('./wrfout_d01_*')
file_names = np.sort(file_names)
file_names = file_names[:3]
file_names

#计算TSK的相关数据
TSK_mean_space = []
TSK_mean_time = []
TSK_mean_lat = []
TSK_mean_lon = []
for file in file_names:
    print(file)
    data = xr.open_dataset(file)
    TSK_mean_space_one = np.mean(data['TSK'],axis=0)
    TSK_mean_time_one = (data['TSK'].mean(dim=['south_north','west_east'])).values
    TSK_mean_lat_one = np.mean(data['TSK'],axis=2)
    TSK_mean_lon_one = np.mean(data['TSK'],axis=1)
    TSK_mean_space.append(TSK_mean_space_one)
    TSK_mean_time.append(TSK_mean_time_one)
    TSK_mean_lat.append(TSK_mean_lat_one)
    TSK_mean_lon.append(TSK_mean_lon_one)

#经向平均
TSK_mean_lon_concat = xr.concat([TSK_mean_lon[0],TSK_mean_lon[1],TSK_mean_lon[2]],'Time')

fig,ax=plt.subplots(1,1,figsize=(2,1),dpi=400)

ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0][0,:],np.arange(0,72,1),TSK_mean_lon_concat,levels=np.arange(260,310,2),cmap='rainbow',extend='both')
ax.set_title('Meridional Average TSK',fontdict={'size':4})
ax.set_ylabel("Simulation hours",fontdict={'size':3})
ax.set_xlabel("lon",fontdict={'size':3})
cb = plt.colorbar(c,shrink=0.9)
cb.set_label('TSK [K]',fontdict={'size':3})
cb.set_ticks([260,270,280,290,300])
cb.ax.tick_params(direction="out",length=1,labelsize=3)
plt.savefig('TSK mean lon.png')
plt.show()

#纬向平均
TSK_mean_lat_concat = xr.concat([TSK_mean_lat[0],TSK_mean_lat[1],TSK_mean_lat[2]],'Time')

fig,ax=plt.subplots(1,1,figsize=(2,1),dpi=400)

ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(np.arange(0,72,1),data['XLAT'][0][:,0],TSK_mean_lat_concat.values.T,levels=np.arange(260,310,2),cmap='rainbow',extend='both')
ax.set_title('Zonal Average TSK',fontdict={'size':4})
ax.set_xlabel("Simulation hours",fontdict={'size':3})
ax.set_ylabel("Lat",fontdict={'size':3})
cb = plt.colorbar(c,shrink=0.9)
cb.set_label('TSK [K]',fontdict={'size':3})
cb.set_ticks([260,270,280,290,300])
cb.ax.tick_params(direction="out",length=1,labelsize=3)
plt.savefig('TSK mean lat.png')
plt.show()

#全球平均温度时间折线
fig,ax = plt.subplots(1,1,figsize=(4,2),dpi=200)

time = np.arange(0,72,1)
ax.plot(time, np.array(TSK_mean_time).flatten(), label="TSK_mean_global", color="blue",  linestyle="-.",linewidth=0.7)
ax.set_xticks([0,12,24,36,48,60,72])
ax.set_xlabel("Simulation hours",fontdict={'size':5})
ax.set_ylabel("TSK [K]",fontdict={'size':5})
ax.set_title('Average TSK',fontdict={'size':5})
ax.tick_params(color = 'gray',length=1,labelsize=5)
plt.savefig('time TSK mean global.png')
plt.show()

#ax.legend()

#画空间分布图
fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=400,subplot_kw={'projection': ccrs.PlateCarree()})


lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())
ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0],data['XLAT'][0],(TSK_mean_space[0]+TSK_mean_space[1]+TSK_mean_space[2])/3.0,levels=np.arange(240,320,2),cmap='rainbow',transform=ccrs.PlateCarree(),extend='both')

ax.coastlines('50m', color='k', lw=0.3)
ax.set_title('The Distribution of Average TSK',fontdict={'size':5})
cb = plt.colorbar(c,shrink=0.25)
cb.set_label('TSK [K]',fontdict={'size':5})
cb.set_ticks([240,255,270,285,300,315])
cb.ax.tick_params(direction="out",length=1,labelsize=5)
plt.savefig('space TSK mean.png')
plt.show()

#计算rain的相关数据
data03  = xr.open_dataset(file_names[2])
RAIN_mean_space = (data03['RAINC'][-1] + data03['RAINNC'][-1])/(24*3.0) #unit:mm/h

QFX_mean_space = []
for file in file_names:
    print(file)
    data = xr.open_dataset(file)
    QFX_mean_space_one = np.mean(data['QFX'],axis=0)
    QFX_mean_space.append(QFX_mean_space_one)#画降水空间分布图

fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=400,subplot_kw={'projection': ccrs.PlateCarree()})


lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())
ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0],data['XLAT'][0],RAIN_mean_space,levels=np.arange(0,1.5,0.02),cmap='rainbow',transform=ccrs.PlateCarree(),extend='both')

ax.coastlines('50m', color='k', lw=0.3)
ax.set_title('The Distribution of Average Rain',fontdict={'size':5})
cb = plt.colorbar(c,shrink=0.25)
cb.set_label('[mm * $h^{-1}$]',fontdict={'size':5})
cb.set_ticks([0.4,0.8,1.2,])
cb.ax.tick_params(direction="out",length=1,labelsize=5)
plt.savefig('space RAIN mean.png')
plt.show()

#画蒸发空间分布图

fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=400,subplot_kw={'projection': ccrs.PlateCarree()})


lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())
ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0],data['XLAT'][0],10000*(QFX_mean_space[0]+QFX_mean_space[1]+QFX_mean_space[2])/3.0,levels=np.arange(0,1.7,0.05),cmap='rainbow',transform=ccrs.PlateCarree(),extend='both')

ax.coastlines('50m', color='k', lw=0.3)
ax.set_title('The Distribution of Average QFX',fontdict={'size':5})
cb = plt.colorbar(c,shrink=0.25)
cb.set_label('[1.0E-4 kg $m^{-2}$ $s^{-1}$]',fontdict={'size':5})
cb.set_ticks([0.4,0.8,1.2])
cb.ax.tick_params(direction="out",length=1,labelsize=5)
plt.savefig('space QFX mean.png')
plt.show()

#计算TOA的相关数据
TOA_mean_space = []
TOA_mean_lat = []
TOA_mean_lon = []
for file in file_names:
    print(file)
    data = xr.open_dataset(file)
    val = data['SWDNT'] - data['SWUPT'] + data['LWDNT'] - data['LWUPT']
    TOA_mean_space_one = np.mean(val,axis=0)
    TOA_mean_lat_one = np.mean(val,axis=2)
    TOA_mean_lon_one = np.mean(val,axis=1)
    TOA_mean_space.append(TOA_mean_space_one)
    TOA_mean_lat.append(TOA_mean_lat_one)
    TOA_mean_lon.append(TOA_mean_lon_one)

#画空间分布图
fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=400,subplot_kw={'projection': ccrs.PlateCarree()})


lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())
ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0],data['XLAT'][0],(TOA_mean_space[0]+TOA_mean_space[1]+TOA_mean_space[2])/3.0,levels=np.arange(-80,195,5),cmap='rainbow',transform=ccrs.PlateCarree(),extend='both')

ax.coastlines('50m', color='k', lw=0.3)
ax.set_title('The Distribution of Average TOA',fontdict={'size':5})
cb = plt.colorbar(c,shrink=0.25)
cb.set_label('TOA [W *$m^{-2}$]',fontdict={'size':5})
cb.set_ticks([-50,0,50,100,150])
cb.ax.tick_params(direction="out",length=1,labelsize=5)
plt.savefig('space TOA mean.png')
plt.show()

#TOA经向平均
TOA_mean_lon_concat = xr.concat([TOA_mean_lon[0],TOA_mean_lon[1],TOA_mean_lon[2]],'Time')
TOA_mean_lon_concat = TOA_mean_lon_concat.mean('Time')

fig,ax = plt.subplots(1,1,figsize=(4,2),dpi=200)

time = np.arange(0,72,1)
ax.plot(data['XLONG'][0][0,:], TOA_mean_lon_concat, label="TOA_mean_lon", color="blue",  linestyle="-.",linewidth=0.7)
#ax.set_xticks([-180,-120,24,36,48,60,72])
ax.set_xlabel("longitude",fontdict={'size':5})
ax.set_ylabel("TOA [W *$m^{-2}$]",fontdict={'size':5})
ax.set_title('Meridional Average TOA',fontdict={'size':5})
ax.tick_params(color = 'gray',length=1,labelsize=5)
plt.savefig('TOA mean lon.png')
plt.show()

#ax.legend()

#TOA纬向平均
TOA_mean_lat_concat = xr.concat([TOA_mean_lat[0],TOA_mean_lat[1],TOA_mean_lat[2]],'Time')
TOA_mean_lat_concat = TOA_mean_lat_concat.mean('Time')

fig,ax = plt.subplots(1,1,figsize=(4,2),dpi=200)

time = np.arange(0,72,1)
ax.plot(data['XLAT'][0][:,0], TOA_mean_lat_concat, label="TOA_mean_lat", color="blue",  linestyle="-.",linewidth=0.7)
#ax.set_xticks([-180,-120,24,36,48,60,72])
ax.set_xlabel("latitude",fontdict={'size':5})
ax.set_ylabel("TOA [W *$m^{-2}$]",fontdict={'size':5})
ax.set_title('Zonal Average TOA',fontdict={'size':5})
ax.tick_params(color = 'gray',length=1,labelsize=5)
plt.savefig('TOA mean lat.png')
plt.show()

#ax.legend()

#计算TOA长波的相关数据
TOA_LW_mean_space = []
TOA_LW_mean_lat = []
TOA_LW_mean_lon = []
for file in file_names:
    print(file)
    data = xr.open_dataset(file)
    val = data['LWDNT'] - data['LWUPT']
    TOA_LW_mean_space_one = np.mean(val,axis=0)
    TOA_LW_mean_lat_one = np.mean(val,axis=2)
    TOA_LW_mean_lon_one = np.mean(val,axis=1)
    TOA_LW_mean_space.append(TOA_LW_mean_space_one)
    TOA_LW_mean_lat.append(TOA_LW_mean_lat_one)
    TOA_LW_mean_lon.append(TOA_LW_mean_lon_one)
    
#画空间分布图
fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=400,subplot_kw={'projection': ccrs.PlateCarree()})


lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())
ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0],data['XLAT'][0],(TOA_LW_mean_space[0]+TOA_LW_mean_space[1]+TOA_LW_mean_space[2])/3.0,levels=np.arange(-300,-90,10),cmap='rainbow',transform=ccrs.PlateCarree(),extend='both')

ax.coastlines('50m', color='k', lw=0.3)
ax.set_title('The Distribution of Average TOA_LW',fontdict={'size':5})
cb = plt.colorbar(c,shrink=0.25)
cb.set_label('TOA [W *$m^{-2}$]',fontdict={'size':5})
cb.set_ticks([-300,-250,-200,-150,-100])
cb.ax.tick_params(direction="out",length=1,labelsize=5)
plt.savefig('space TOA_LW mean.png')
plt.show()


#TOA经向平均
TOA_LW_mean_lon_concat = xr.concat([TOA_LW_mean_lon[0],TOA_LW_mean_lon[1],TOA_LW_mean_lon[2]],'Time')
TOA_LW_mean_lon_concat = TOA_LW_mean_lon_concat.mean('Time')

fig,ax = plt.subplots(1,1,figsize=(4,2),dpi=200)

time = np.arange(0,72,1)
ax.plot(data['XLONG'][0][0,:], TOA_LW_mean_lon_concat, label="TOA_LW_mean_lon", color="blue",  linestyle="-.",linewidth=0.7)
#ax.set_xticks([-180,-120,24,36,48,60,72])
ax.set_xlabel("longitude",fontdict={'size':5})
ax.set_ylabel("TOA [W *$m^{-2}$]",fontdict={'size':5})
ax.set_title('Meridional Average TOA_LW',fontdict={'size':5})
ax.tick_params(color = 'gray',length=1,labelsize=5)
plt.savefig('TOA_LW mean lon.png')
plt.show()

#ax.legend()

#TOA纬向平均
TOA_LW_mean_lat_concat = xr.concat([TOA_LW_mean_lat[0],TOA_LW_mean_lat[1],TOA_LW_mean_lat[2]],'Time')
TOA_LW_mean_lat_concat = TOA_LW_mean_lat_concat.mean('Time')

fig,ax = plt.subplots(1,1,figsize=(4,2),dpi=200)

time = np.arange(0,72,1)
ax.plot(data['XLAT'][0][:,0], TOA_LW_mean_lat_concat, label="TOA_LW_mean_lat", color="blue",  linestyle="-.",linewidth=0.7)
#ax.set_xticks([-180,-120,24,36,48,60,72])
ax.set_xlabel("latitude",fontdict={'size':5})
ax.set_ylabel("TOA [W *$m^{-2}$]",fontdict={'size':5})
ax.set_title('Zonal Average TOA_LW',fontdict={'size':5})
ax.tick_params(color = 'gray',length=1,labelsize=5)
plt.savefig('TOA_LW mean lat.png')
plt.show()

#ax.legend()

#计算TOA长波的相关数据
TOA_SW_mean_space = []
TOA_SW_mean_lat = []
TOA_SW_mean_lon = []
for file in file_names:
    print(file)
    data = xr.open_dataset(file)
    val = data['SWDNT'] - data['SWUPT']
    TOA_SW_mean_space_one = np.mean(val,axis=0)
    TOA_SW_mean_lat_one = np.mean(val,axis=2)
    TOA_SW_mean_lon_one = np.mean(val,axis=1)
    TOA_SW_mean_space.append(TOA_SW_mean_space_one)
    TOA_SW_mean_lat.append(TOA_SW_mean_lat_one)
    TOA_SW_mean_lon.append(TOA_SW_mean_lon_one)
    
#画空间分布图
fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=400,subplot_kw={'projection': ccrs.PlateCarree()})


lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())
ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0],data['XLAT'][0],(TOA_SW_mean_space[0]+TOA_SW_mean_space[1]+TOA_SW_mean_space[2])/3.0,levels=np.arange(100,400,20),cmap='rainbow',transform=ccrs.PlateCarree(),extend='both')

ax.coastlines('50m', color='k', lw=0.3)
ax.set_title('The Distribution of Average TOA_SW',fontdict={'size':5})
cb = plt.colorbar(c,shrink=0.25)
cb.set_label('TOA [W *$m^{-2}$]',fontdict={'size':5})
cb.set_ticks([100,150,200,250,300,350,400])
cb.ax.tick_params(direction="out",length=1,labelsize=5)
plt.savefig('space TOA_SW mean.png')
plt.show()


#TOA经向平均
TOA_SW_mean_lon_concat = xr.concat([TOA_SW_mean_lon[0],TOA_SW_mean_lon[1],TOA_SW_mean_lon[2]],'Time')
TOA_SW_mean_lon_concat = TOA_SW_mean_lon_concat.mean('Time')

fig,ax = plt.subplots(1,1,figsize=(4,2),dpi=200)

time = np.arange(0,72,1)
ax.plot(data['XLONG'][0][0,:], TOA_SW_mean_lon_concat, label="TOA_SW_mean_lon", color="blue",  linestyle="-.",linewidth=0.7)
#ax.set_xticks([-180,-120,24,36,48,60,72])
ax.set_xlabel("longitude",fontdict={'size':5})
ax.set_ylabel("TOA [W *$m^{-2}$]",fontdict={'size':5})
ax.set_title('Meridional Average TOA_SW',fontdict={'size':5})
ax.tick_params(color = 'gray',length=1,labelsize=5)
plt.savefig('TOA_SW mean lon.png')
plt.show()

#ax.legend()

#TOA纬向平均
TOA_SW_mean_lat_concat = xr.concat([TOA_SW_mean_lat[0],TOA_SW_mean_lat[1],TOA_SW_mean_lat[2]],'Time')
TOA_SW_mean_lat_concat = TOA_SW_mean_lat_concat.mean('Time')

fig,ax = plt.subplots(1,1,figsize=(4,2),dpi=200)

time = np.arange(0,72,1)
ax.plot(data['XLAT'][0][:,0], TOA_SW_mean_lat_concat, label="TOA_SW_mean_lat", color="blue",  linestyle="-.",linewidth=0.7)
#ax.set_xticks([-180,-120,24,36,48,60,72])
ax.set_xlabel("latitude",fontdict={'size':5})
ax.set_ylabel("TOA [W *$m^{-2}$]",fontdict={'size':5})
ax.set_title('Zonal Average TOA_SW',fontdict={'size':5})
ax.tick_params(color = 'gray',length=1,labelsize=5)
plt.savefig('TOA_SW mean lat.png')
plt.show()

#ax.legend()

#风800hpq空间分布图
z_list=[8000.,5000.,2000.]
result_U=[]
result_V=[]
for file in file_names:
    print(file)
    ds = xr.open_dataset(file)
    U = destagger(ds['U'],3,meta=True)
    V = destagger(ds['V'],2,meta=True)
    P = ds['PB']+ds['P']
    U_inter = interpz3d(U,P,np.array(z_list))
    V_inter = interpz3d(V,P,np.array(z_list))
    result_U.append(U_inter)
    result_V.append(V_inter)
    

#画流场
fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=200,subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})

lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())

ax.tick_params(color = 'gray',direction='in',labelsize=5,length=0.9)

ax.streamplot(data['XLONG'][0].values,data['XLAT'][0].values,result_U[-1][-1,0,:,:].values ,(xr.concat(result_V,'Time').mean('Time'))[0].values, transform=ccrs.PlateCarree(),linewidth=0.7,arrowsize=0.5,density=1)

ax.coastlines('50m', color='k', lw=0.15)

ax.set_title('Wind (800hpa)',fontdict={'size':10})
plt.savefig('Wind (800hpa).png')
plt.show()

#计算水汽柱浓度

result=[]
for file in file_names:
    #读取数据
    ds = xr.open_dataset(file)
    #提取计算水汽柱浓度的变量
    P = ds['PH']+ds['PHB']
    gmp=P/9.81
    z = gmp[:,1:31,:,:]-gmp[:,0:30,:,:]
    z = z.rename({"bottom_top_stag": "bottom_top"})
    #计算水汽柱浓度
    q_sumone = ((1.0/ds['ALT'])*ds['QVAPOR']*z).sum(axis=1)  #unit:kg/m2
    result.append(q_sumone)
    
q_all = xr.concat(result,dim='Time')
q_mean = q_all.mean(axis=0)

#水汽柱浓度空间分布图
fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=400,subplot_kw={'projection': ccrs.PlateCarree()})


lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())
ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0],data['XLAT'][0],q_mean,levels=np.arange(0,65,2),cmap='rainbow',transform=ccrs.PlateCarree(),extend='both')

ax.coastlines('50m', color='k', lw=0.3)
ax.set_title('Qvapor Column Concentration',fontdict={'size':5})
cb = plt.colorbar(c,shrink=0.25)
cb.set_label('[kg/m2]',fontdict={'size':5})
cb.ax.tick_params(direction="out",length=1,labelsize=5)
cb.set_ticks([0,10,20,30,40,50,60])
plt.savefig('qvapor column concentration.png')
plt.show()

#水汽800hpq空间分布图
z_list=[8000.,]
result_q=[]
for file in file_names:
    print(file)
    ds = xr.open_dataset(file)
    qvapor = ds['QVAPOR']
    P = ds['PB']+ds['P']
    qvapor_800 = interpz3d(qvapor,P,np.array(z_list))
    result_q.append(qvapor_800)
    

fig,ax=plt.subplots(1,1,figsize=(4,4),dpi=400,subplot_kw={'projection': ccrs.PlateCarree()})


lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_extent([-180,180,-20,75],crs = ccrs.PlateCarree()) 
ax.set_xticks(np.arange(-180,181,30), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(-20,76,15), crs=ccrs.PlateCarree())
ax.tick_params(color = 'gray',direction='in',length=1,labelsize=3)

c = ax.contourf(data['XLONG'][0],data['XLAT'][0],1.0e6*(xr.concat(result_q,'Time').mean('Time')[0]),levels=np.arange(1.5,4,0.1),cmap='rainbow',transform=ccrs.PlateCarree(),extend='both')

ax.coastlines('50m', color='k', lw=0.3)
ax.set_title('The Distribution of Avarage Qvapor (800 hpa)',fontdict={'size':5})
cb = plt.colorbar(c,shrink=0.25)
cb.set_label('[1.0E-6 kg kg-1]',fontdict={'size':5})
cb.ax.tick_params(direction="out",length=1,labelsize=5)
cb.set_ticks([1.5,2.0,2.5,3.0,3.5,4.0])
plt.savefig('qvapor 800 hpa')
plt.show()

'''
result=[]
for file in file_names:
    #读取数据
    ds = xr.open_dataset(file)
    #提取计算水汽柱浓度的变量
    P = ds['P']+ds['PB']
    g = 9.81
    P_top = (ds['P_TOP'].values[0])*np.ones((1,90,180))
    P_TOP = xr.DataArray(P_top, ds['PSFC'].coords,ds['PSFC'].dims)
    P = (xr.concat([ds['PSFC'].rename({"Time": "bottom_top"}),((P[0,0:29,:,:]+P[0,1:30,:,:])/2.0),P_TOP.rename({"Time": "bottom_top"})],'bottom_top'))
    P = P[0:30,:,:] -P[1:31,:,:]
    #计算水汽柱浓度
    q_sumone = ((ds['QVAPOR'][0]/g)*P).sum(axis=0)  #unit:kg/m2
    result.append(q_sumone)

q_all = xr.concat(result,dim='Time')
q_mean = q_all.mean(axis=0)
'''