时间剖面图？so easy
核心函数：mpcalc.divergence
前言
在本文中，我们将利用WRFOUT数据进行处理和分析，并生成直观明了的时间剖面图。你将能够清楚地看到水汽通量散度随着时间和高度的变化趋势，从而更好地理解大气中水汽的传播与运动机制
1. 导入库与读取数据
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()
wrf_list =
2. 提取变量
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)
3. 循环计算逐个时次水汽通量
%time
p = getvar(wrf_list, 'pressure', timeidx=ALL_TIMES, method='cat')
levels =
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,lev.shape,lat.shape,lon.shape))
CPU times: user 2 µs, sys: 0 ns, total: 2 µs
Wall time: 5.48 µs
lon.shape,lat.shape,q_flux_divergence_all.shape
((9, 90, 90), (9, 90, 90), (9, 8, 90, 90))
%time
for j in range(time.shape):
for i in range(lev.shape):
q_flux_divergence_all = mpcalc.divergence(u = to_np(qu),v = to_np(qv),dx = to_np(dx) ,dy = to_np(dy))
#将 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
CPU times: user 2 µs, sys: 0 ns, total: 2 µs
Wall time: 5.72 µs
Out:
(9, 8, 90, 90)
metpy的mpcalc.divergence中，x轴与y轴的默认排序是x_dim=-1, y_dim=-2，即(：，：，Y，X)，还请留意
4. 绘图部
qfd = q_flux_divergence_all[:,:,35,35]
latlon= xy_to_ll(wrf_list, 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()
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()
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()
<xarray.DataArray 'latlon' (lat_lon: 2)>
array([ 30.16124186, 104.7080441 ])
Coordinates:
xy_coordobject CoordPair(x=35, y=35)
* lat_lon (lat_lon) <U3 'lat' 'lon'
小结
1、可自行探索对应站点经纬度的高度时间剖面，wrfpython有换算xy与经纬度的函数
2、可对某一纬度进行平均后再绘图分析
3、优化方向可以是计算速度的提升，例如使用dask或者向量化，懂的同学可手动优化
完整代码与文件在此
文章来源于微信公众号：气ython风雨
页:
[1]