时间剖面图？so easy

核心函数：mpcalc.divergence

前言

在本文中，我们将利用WRFOUT数据进行处理和分析，并生成直观明了的时间剖面图。你将能够清楚地看到水汽通量散度随着时间和高度的变化趋势，从而更好地理解大气中水汽的传播与运动机制

1. 导入库与读取数据

代码语言：javascript复制

from wrf import uvmet, to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords,ALL_TIMES,xy_to_ll
import numpy as np
from netCDF4 import Dataset
import xarray as xr
from metpy.units import units
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
from matplotlib.colors import from_levels_and_colors
import cmaps
from glob import glob
import metpy.calc as mpcalc
import metpy.constants as constants
import os
# 定义 WRF 文件夹路径和文件名前缀
wrfout_path = "/home/mw/input/wrfout3385/"
filename_prefix = "wrfout_d02_"

wrf_files = sorted([os.path.join(wrfout_path, f) for f in os.listdir(wrfout_path) if f.startswith(filename_prefix)])
wrf_list = [Dataset(f) for f in wrf_files]

代码语言：javascript复制

2. 提取变量

代码语言：javascript复制

lon = getvar(wrf_list, 'lon',timeidx=ALL_TIMES, method='cat')
lat = getvar(wrf_list, 'lat', timeidx=ALL_TIMES, method='cat')
u = getvar(wrf_list, 'ua', timeidx=ALL_TIMES, method='cat')
v = getvar(wrf_list, 'va', timeidx=ALL_TIMES, method='cat')
#在wrf中q单位是kg/kg
q = getvar(wrf_list, 'QVAPOR',timeidx=ALL_TIMES, method='cat')*1000
# 提取WRF模拟的经纬度数组
lats, lons = latlon_coords(u)
time = u.Time
# 提取WRF模拟的地图投影
wrf_proj = get_cartopy(u)

代码语言：javascript复制

3. 循环计算逐个时次水汽通量

代码语言：javascript复制

%time
p = getvar(wrf_list, 'pressure', timeidx=ALL_TIMES, method='cat')
levels = [1000,925,850,700,600,500,400,300]
u_interp = interplevel(u, p, levels)
v_interp = interplevel(v, p, levels)
q_interp = interplevel(q, p, levels) ** units('g/kg')
lev = u_interp.level
dx, dy = mpcalc.lat_lon_grid_deltas(lon, lat)
# 计算水汽通量散度
qu = u_interp *q_interp/constants.g
qv = v_interp *q_interp/constants.g
q_flux_divergence_all = np.zeros((time.shape[0],lev.shape[0],lat.shape[2],lon.shape[1]))

代码语言：javascript复制

CPU times: user 2 µs, sys: 0 ns, total: 2 µs
Wall time: 5.48 µs

代码语言：javascript复制

lon.shape,lat.shape,q_flux_divergence_all.shape

代码语言：javascript复制

((9, 90, 90), (9, 90, 90), (9, 8, 90, 90))

代码语言：javascript复制

%time
for j in range(time.shape[0]):
    for i in range(lev.shape[0]):
        q_flux_divergence_all[j,i] = mpcalc.divergence(u = to_np(qu[j,i]),v = to_np(qv[j,i]),dx = to_np(dx[j]) ,dy = to_np(dy[j]))

#将 q_flux_divergence_all 中的 NaN 值替换为 0
q_flux_divergence_all = np.nan_to_num(q_flux_divergence_all, nan=0)
q_flux_divergence_all.shape

代码语言：javascript复制

CPU times: user 2 µs, sys: 0 ns, total: 2 µs
Wall time: 5.72 µs

Out[7]:

代码语言：javascript复制

(9, 8, 90, 90)

metpy的mpcalc.divergence中，x轴与y轴的默认排序是x_dim=-1, y_dim=-2，即(：，：，Y，X)，还请留意

4. 绘图部

代码语言：javascript复制

qfd = q_flux_divergence_all[:,:,35,35]
latlon= xy_to_ll(wrf_list[0], 35, 35)
print(latlon)
fig, ax = plt.subplots(figsize=(15, 15))
ax.set_title('lev-time', loc='center', fontsize=20)
ax.set_yscale('symlog')
ax.set_yticks([1000, 925, 850, 700, 600, 500, 400, 300])
ax.set_yticklabels(['1000', '925', '850', '700', '600', '500', '400', '300'], fontsize=16)
ax.invert_yaxis()
ax.set_ylabel('Level (hPa)', fontsize=18)
ax.set_xlabel('TIME', fontsize=18)
c = ax.contourf(time,lev ,qfd.T,levels=np.arange(-50,55,5), extend='both',zorder=0, cmap=plt.cm.bwr)
position = fig.add_axes([0.95, 0.16, 0.05, 0.65])
position.tick_params(axis='both', which='major', labelsize=14, width=2, length=6)
position.set_ylabel('Q Flux Divergence', fontsize=16)
fig.colorbar(c,cax=position,orientation='vertical',format='%d')
plt.show()

代码语言：javascript复制

<xarray.DataArray 'latlon' (lat_lon: 2)>
array([ 30.16124186, 104.7080441 ])
Coordinates:
    xy_coord  object CoordPair(x=35, y=35)
  * lat_lon   (lat_lon) <U3 'lat' 'lon'

小结

可自行探索对应站点经纬度的高度时间剖面，wrfpython有换算xy与经纬度的函数
可对某一纬度进行平均后再绘图分析
优化方向可以是计算速度的提升，例如使用dask或者向量化，懂的同学可手动优化

完整文件与代码在此

排序优化 import time 函数

0 人点赞