GSI同化前后WRF模式增量可视化（python）_gsi同化模型

一、初始场（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()更智能，但速度相对较慢。 
 考虑到上述带星号（ * ）的项目，请评估不同的DA方法在估算真相方面的表现。
DAPPER可以通过各种典型的测试用例和统计数据对DA方法进行数值研究。 它（a）复制了文献中报告的数字基准结果，并且（b）促进了比较研究，从而提高了（a）结果的可靠性和（b）相关性。 例如，此图是由examples/basic_3.py生成的，并且是的图5.7的复制品。
 DAPPER是（c）用Python编写的开放源代码，并且（d）注重可读性； 这促进了基础科学的（c）复制和（d）传播，并使其易于适应和扩展。 它还带有一系列的诊断和统计信息，以及实时绘图（与同化一致）
 导师将分发以协助练习， 并对每个部分进行总结后再进行总结。
本地工作说明
您也可以在自己的（Linux / Windows / Mac）计算机上运行这些笔记本。 这比在线运行它们要快一些。
 先决条件：Python> = 3.6。 如果您不是python专家： 1a。 通过安装Python。 1b。 使用运行以下命令。 1c。 （可选） 。 如果安装（以下）失败，请首先尝试执行步骤1c。
 安装： 在终端中运行以下命令（不包括$符号）： $ git clone https://github.com/nansencenter/DA-tutori
函数或变量 'm_proj' 无法识别。
出错 gebco0704 (第 10 行)
m_proj('mercator','lon',[160 200],'lat',[50 65]);
是不是m_proj是您自己写的某个函数吗？如果是的话方便参考下吗，十分感谢
                用Matlab画GEBCO与ETOPO1的海拔分布
                    阿翔_ocean: 
                    你下错区域了吧，海拔全大于0导致的全变nan值了吧
                用Matlab画GEBCO与ETOPO1的海拔分布
                    Leo_warfare: 
                    您好，博主，我画出来的二维海深图是纯绿色的，并没有岸线的显示，是因为没有导入岸线数据么？