前言
数据同化是通过插值、变分、滤波等操作去改变数值模式初始场从而改变预报结果，有时候为了观察同化的效果需要对同化前后的初始场和预报场增量进行可视化。本次实验的数值模式方案与前一篇文章一致资料同化简单评估，但时间不同，同化的数据是风机上的观测资料。
一、初始场（wrf_input）增量可视化
水平增量
选取了表面压强、最底层温度和经纬向风速进行增量的可视化
import numpy as np
from netCDF4 import Dataset
import cartopy.crs as ccrs
import matplotlib
matplotlib.use('Agg') # 使用Agg后端进行文件保存
import matplotlib.pyplot as plt
import cartopy.mpl.ticker as cticker
from matplotlib.colors import LinearSegmentedColormap
import os
from wrf import to_np, getvar, get_cartopy, latlon_coords
# 数据路径
f1 = Dataset("/Users/xychen/assimilation/GSI/GSI-Docker/2019080112/bkg/wrfinput_d01")
f2 = Dataset("/Users/xychen/assimilation/GSI/GSI-Docker/real_case_/wrf_inout")
save_path = "/Users/xychen/pythonpainting/input_delta/"
os.makedirs(save_path, exist_ok=True)
def get_symmetric_levels(data, num_levels=10):
"""生成对称的levels并确保0值在白色上"""
range_val = max(abs(np.min(data)), abs(np.max(data)))
return np.linspace(-range_val, range_val, num_levels)
def plot_delta(data1, data2, var_name, levels, cmap, ax, title):
"""绘制差值图像"""
delta = data2 - data1
if delta.ndim > 2:
delta = delta[0, :, :] # 选择第一个时间步和第一个层
c = ax.contourf(to_np(lons), to_np(lats), delta, levels=levels, cmap=cmap, extend='both', transform=ccrs.PlateCarree())
ax.set_title(title, fontsize=12)
fig.colorbar(c, ax=ax, label=var_name)
# 画图
fig = plt.figure(figsize=(16, 12))
lats, lons = latlon_coords(getvar(f1, "tk", 0))
cart_proj = get_cartopy(getvar(f1, "tk", 0))
# 对称levels
variables = {
"slp": ("slp", "delt_p"),
"T2": ("tk", "delt_T"),
"u": ("ua", "delt_U"),
"v": ("va", "delt_V")
}
levels_dict = {}
for var, (var_name, title) in variables.items():
data1 = getvar(f1, var_name, 0).data
data2 = getvar(f2, var_name, 0).data
levels_dict[var] = get_symmetric_levels(data2 - data1)
ax = fig.add_subplot(2, 2, list(variables.keys()).index(var) + 1, projection=cart_proj)
ax.coastlines('50m', linewidth=0.8)
ax.gridlines(color="black", linestyle="dotted")
ax.set_xticks(np.arange(111.5, 124.5, 2), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(21, 32, 2), crs=ccrs.PlateCarree())
ax.set_xlabel('longitude', fontsize=10)
ax.set_ylabel('latitude', fontsize=10)
ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())
myCmap = LinearSegmentedColormap.from_list("", ["#0000FF", "#FFFFFF", "#FF0000"])
plot_delta(data1, data2, var_name, levels_dict[var], myCmap, ax, title)
# 保存图片
plt.savefig(os.path.join(save_path, "input_delta_plot.png"), dpi=800)
plt.close()
print('ok了，老铁')
效果如下
可以通过风场、气温和压强的变化看出GSI同化时较好的物理约束
此次同化同化了海表20m、70m、100m的风机资料（风场、压强、温度），为了观察垂直剖面上的初始场变化，所以要画剖面的增量图
import os
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from wrf import getvar, vertcross, CoordPair, to_np, interpline
from matplotlib.colors import TwoSlopeNorm
# 读取 f1 和 f2 数据集
f1 = Dataset("/Users/xychen/assimilation/GSI/GSI-Docker/2019080112/bkg/wrfinput_d01")
f2 = Dataset("/Users/xychen/assimilation/GSI/GSI-Docker/real_case_/wrf_inout")
# 提取变量
variables = ["ua", "va", "tk", "rh"] # U风速, V风速, 温度, 相对湿度
# 定义剖面线
start_point = CoordPair(lat=25.8, lon=111.5)
end_point = CoordPair(lat=25.8, lon=124.5)
# 保存路径
save_path = "/Users/xychen/pythonpainting/input_delta/"
os.makedirs(save_path, exist_ok=True)
def get_symmetric_levels(data):
"""自动生成对称的levels并确保0值在白色上"""
min_val, max_val = np.min(data), np.max(data)
range_val = max(abs(min_val), abs(max_val))
return np.linspace(-range_val, range_val, 21)
# 计算并绘制差值剖面图
for var_name in variables:
var1 = getvar(f1, var_name, timeidx=0)
var2 = getvar(f2, var_name, timeidx=0)
# 计算差值
delta_var = var2 - var1
# 计算沿剖面线的垂直截面
cross_delta = vertcross(delta_var, getvar(f1, "z", timeidx=0, units="m"), wrfin=f1, start_point=start_point, end_point=end_point, latlon=True, meta=True)
# 获取经度和高度的坐标
lons = to_np(interpline(getvar(f1, "lon", timeidx=0), wrfin=f1, start_point=start_point, end_point=end_point))
heights = to_np(cross_delta.coords["vertical"])
# 自动设置对称的levels
levels = get_symmetric_levels(to_np(cross_delta))
# 设置 colorbar 使其在0值处为白色
norm = TwoSlopeNorm(vmin=levels.min(), vcenter=0, vmax=levels.max())
# 绘制剖面图
plt.figure(figsize=(10, 6))
contour = plt.contourf(lons, heights, to_np(cross_delta), cmap="bwr", extend='both', norm=norm)
plt.colorbar(contour, label=f"{var_name} Difference")
plt.title(f"{var_name} Difference Profile at Latitude 25.8°")
plt.xlabel("Longitude")
plt.ylabel("Height (m)")
# 保存图片
plt.savefig(os.path.join(save_path, f"{var_name}_delta_profile.png"))
plt.close()
print('图画好了')
效果如下（这里就只展示温度的增量剖面图）
二、预报场（wrf_output）增量可视化
预报结果增量可视化
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
from matplotlib.colors import LinearSegmentedColormap
from wrf import getvar, get_cartopy, latlon_coords
import os
def create_cmap():
"""创建对称颜色图谱"""
return LinearSegmentedColormap.from_list("", ["#0000FF", "#FFFFFF", "#FF0000"])
def plot_contour(ax, lons, lats, data, levels, cmap, title, colorbar_label):
"""绘制等值线填充图"""
c = ax.contourf(lons, lats, data, levels=levels, cmap=cmap, extend='both', transform=ccrs.PlateCarree())
ax.set_title(title, fontsize=12)
colorbar = plt.colorbar(c, ax=ax, orientation='vertical', pad=0.05, label=colorbar_label)
return colorbar
def set_map_features(ax):
"""设置地图特征"""
ax.coastlines('50m', linewidth=0.8)
ax.gridlines(color="black", linestyle="dotted")
ax.set_xticks(np.arange(115, 129, 2), crs=ccrs.PlateCarree())
ax.set_yticks(np.arange(27, 39, 2), crs=ccrs.PlateCarree())
ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())
def main():
f1 = Dataset("/Users/xychen/assimilation/GSI/GSI-Docker/2024062812/bkg/wrfout_d01_2024-06-28_12:00:00")
f2 = Dataset("/Users/xychen/WRF/run/wrfout_d01_2024-06-28_12:00:00")
save_path = "/Users/xychen/pythonpainting/output_delta_pic/"
os.makedirs(save_path, exist_ok=True)
cmap = create_cmap()
for i in range(0, 49, 12):
day = (i + 12) // 24 + 28
h = (i + 12) % 24
timestr = f"{day}.{h:02d}"
name = f"{timestr}delta.png"
# 读取数据
T = getvar(f1, "tk", i)
slp1 = getvar(f1, "slp", i, units="Pa").data
slp2 = getvar(f2, "slp", i, units="Pa").data
T2m1 = getvar(f1, "tk", i)[0, :, :].data
T2m2 = getvar(f2, "tk", i)[0, :, :].data
u1 = getvar(f1, "ua", i, units="m s-1")[0, :, :].data
u2 = getvar(f2, "ua", i, units="m s-1")[0, :, :].data
v1 = getvar(f1, "va", i, units="m s-1")[0, :, :].data
v2 = getvar(f2, "va", i, units="m s-1")[0, :, :].data
# 计算差值
delt_slp = slp2 - slp1
delt_T2 = T2m2 - T2m1
delt_u = u2 - u1
delt_v = v2 - v1
# 获取地图投影
lats, lons = latlon_coords(T)
cart_proj = get_cartopy(T)
# 创建图形
fig = plt.figure(figsize=(16, 12))
# 绘制每个图形
ax1 = fig.add_subplot(2, 2, 1, projection=cart_proj)
set_map_features(ax1)
plot_contour(ax1, lons, lats, delt_slp, np.linspace(np.min(delt_slp), np.max(delt_slp), 21), cmap, f"{timestr}delt_slp", 'delt_slp')
ax2 = fig.add_subplot(2, 2, 2, projection=cart_proj)
set_map_features(ax2)
plot_contour(ax2, lons, lats, delt_T2, np.arange(-2, 2.1, 0.2), cmap, f"{timestr}delt_T", 'delt_T2')
ax3 = fig.add_subplot(2, 2, 3, projection=cart_proj)
set_map_features(ax3)
plot_contour(ax3, lons, lats, delt_u, np.arange(-3.5, 4, 0.5), cmap, f"{timestr}delt_U", 'delt_U')
ax4 = fig.add_subplot(2, 2, 4, projection=cart_proj)
set_map_features(ax4)
plot_contour(ax4, lons, lats, delt_v, np.arange(-3.5, 4, 0.5), cmap, f"{timestr}delt_V", 'delt_V')
# 保存图片
plt.savefig(os.path.join(save_path, name), dpi=400)
plt.close()
print('图画好了')
if __name__ == "__main__":
main()
效果图如下
三、附录
1.wrf-python
WRF-Python 是一个用于处理和分析WRF（Weather Research and Forecasting）模型输出数据的Python库。它提供了一系列工具和函数，方便用户从WRF模式输出的NetCDF文件中提取数据、进行后处理和可视化。WRF-Python特别适合气象学家和研究人员用来分析模拟结果、生成统计数据和制作图表。通过集成NetCDF4、Numpy、Matplotlib等Python库，WRF-Python能够处理复杂的数据操作任务，并帮助用户进行深入的气象数据分析。
安装后有很多好用的函数，例如提取向量的getvar()（配有向量表格）、wrf.interplevel()函数用于将三维场数据插值到指定的水平面（例如850 hPa温度）上，采用的是简单且快速的线性插值方法。wrf.vertcross()函数可以将三维场数据插值到垂直平面上，这个垂直平面是通过用户指定的水平线来确定的（即横截面）。wrf.interpline()函数用于将二维场数据插值到用户指定的线上。wrf.vinterp()函数则可以将三维场数据插值到用户指定的“表面”层（例如等θ-e层），虽然这个方法比wrf.interplevel()更智能，但速度相对较慢。
扫一扫
5060
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
