GRIB笔记：使用Matplotlib绘制剖面图
之前在两篇文章中分别介绍使用 nwpc-data 封装的 cfgrib 和 eccodes-python 接口加载绘制剖面图需要的数据。
文章中使用的接口将 GRIB 2 消息加载为 xarray.DataArray 对象，并使用 xarray 集成的绘图接口绘制剖面图。但大部分应用都没有必要增加 xarray 库，虽然笔者强烈建议使用 xarray 库处理要素场数据。
nwpc-data 还提供对 eccodes-python 的另外一套封装函数，返回 eccodes-python 的 GRIB 2 消息对象，即代表 handler 的整数。
本文介绍如何利用 nwpc-data 加载 GRIB 2 消息，并使用 matplotlib 绘制 Zonal Mean 偏差的剖面图。
数据计算
计算某个预报时效的月平均偏差可以通过计算每个预报时次的偏差再求偏差的平均得到。
单个偏差
计算单个时次的偏差需要加载两个要素场并作差。
下面的代码从 fcstdata 和 analdata 两个文件中加载层次为 plev 的 var 要素场，并将偏差累加到 bias 中。
每个要素场都是形状为 [121, 240] 的二维数组。
import eccodes
from nwpc_data.grib.eccodes import load_message_from_file
t = load_message_from_file(
file_path=fcstdata,
parameter=var,
level_type="isobaricInhPa",
level=plev,
)
fcst = eccodes.codes_get_double_array(t, "values")
fcst = fcst.reshape([121, 240])
t = load_message_from_file(
file_path=analdata,
parameter=var,
level_type="isobaricInhPa",
level=plev,
)
anal = eccodes.codes_get_double_array(t, "values")
anal = anal.reshape([121, 240])
bias += fcst - anal所有层次偏差
将所有层次的要素场偏差保存到一个三维数组中。
该数组使用下面的代码初始化，其中 plevs 是层次列表。
import numpy as np
plevs = [1000,850,700,500,400,300,250,100,70,50,20,10]
biasm = np.full((len(plevs), 121, 240), np.nan)当我们获取某个层次所有偏差累加得到的 bias 后，求均值并赋值给该数组。
其中 time_count 是起报时次个数，level_index 是层次序号。
bias = bias / time_count
biasm[level_index,] = bias求纬向平均
沿最后一个维度求平均，得到剖面图的最终数据
data = biasm.mean(axis=2)绘图
坐标信息
需要注意，数据的纬度是按照从 90 到 -90 排列的
xi = np.flip(np.linspace(-90, 90, 121))
yi = plevs
z = data准备 x 轴的标签序列
def format_x(x):
if x<0:
return f"{abs(x):.0f}S"
elif x>0:
return f"{x:.0f}N"
else:
return "EQ"
xlabels = [format_x(x) for x in np.arange(-90, 91, 10)]绘制
绘制填充图
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
fig, ax = plt.subplots(figsize=(12, 6))
cf = ax.contourf(xi, yi, z, levels=10, cmap='bwr')处理 x 轴
ax.set_xticks(np.arange(-90, 91, 10))
ax.xaxis.set_major_formatter(
mticker.FixedFormatter(xlabels)
)处理 y 轴，使用对数坐标，并从大到小排列。
关闭 minor 标签，只显示 plevs 中的值。
ax.invert_yaxis()
ax.set_yscale("log")
ax.set_yticks(plevs)
ax.yaxis.set_major_formatter(mticker.ScalarFormatter())
ax.yaxis.set_minor_locator(mticker.NullLocator())增加标注
ax.set_title("Zonal Mean")
fig.colorbar(cf, ax=ax)绘图结果如下图所示
