从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()