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)
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= #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>)
</i>ac0=ax.contourf(to_np(lons), to_np(lats), 1000.0*to_np(o3_z), levels=np.arange(12,54,3),
transform=ccrs.PlateCarree(),
cmap=cmaps.BlGrYeOrReVi200,extend='both')
ax.set_ylabel(rowlabels,size=10)
plt.suptitle(title,fontsize='small',fontweight='heavy')
plt.savefig('test.png')
plt.show()
结果
文章来源于微信公众号：自学气象人
页:
[1]