WRF-Python库插值到指定离地高度层并绘图

任务

从wrfout文件中提取o3变量，并将其数据结果插值到想要的离地高度层上(示例中是1km、3km、5km、10km)，进行可视化。下面提供示例代码，其中用虚线----框住的部分是插值的关键代码。

代码

import numpy as np    
import xarray as xr   
from netCDF4 import Dataset
from wrf import to_np, getvar, get_cartopy, cartopy_xlim,cartopy_ylim, latlon_coords,interpz3d,destagger

import matplotlib as mpl
import matplotlib.pyplot as plt  
from matplotlib.cm import get_cmap
from matplotlib.colors import ListedColormap  
import cartopy.crs as ccrs       
import cartopy.feature as cfeature 
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter, LatitudeLocator, LongitudeLocator)
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
import matplotlib.patches as mpatches
import matplotlib.transforms as mtransforms
from matplotlib import ticker as mticker
import cmaps

def mapart(ax): 
    ax.coastlines('50m', color='k', lw=0.3)                
    ax.set_xlim(cartopy_xlim(o3))
    ax.set_ylim(cartopy_ylim(o3))
    ax.set_yticks([])

title='O3 concentrations above surface'
rowlabels=[ '10km','5km', '3km','1km']      

#------------------------------------------------------
ncfile = Dataset('wrfout_d01_2018-08-02_00:00:00')

ph=getvar(ncfile, "PH",timeidx=0)[:,10:140,10:140]
phb=getvar(ncfile, "PHB",timeidx=0)[:,10:140,10:140]
hgt=getvar(ncfile, "HGT",timeidx=0)[10:140,10:140]
o3 = getvar(ncfile, "o3",timeidx=0)[:,10:140,10:140]

P=ph+phb
P = destagger(P,0,meta=True)
gmp=P/9.81-hgt

z_list=[10000.,5000.,3000.,1000.] #unit:m

o3_z = interpz3d(o3,gmp,np.array(z_list))
#--------------------------------------------------------

lats, lons = latlon_coords(o3)
cart_proj = get_cartopy(o3)

fig,ax=plt.subplots(4,1,figsize=(6,6),dpi=300,subplot_kw={'projection': cart_proj})

for i in range(len(z_list)):
    
    mapart(ax[i]) 
    ac0=ax[i].contourf(to_np(lons), to_np(lats), 1000.0*to_np(o3_z[i]), levels=np.arange(12,54,3),
             transform=ccrs.PlateCarree(),
             cmap=cmaps.BlGrYeOrReVi200,extend='both')
    
    ax[i].set_ylabel(rowlabels[i],size=10)     
 
plt.suptitle(title,fontsize='small',fontweight='heavy')
plt.savefig('test.png')
plt.show()