Piggy_Packages V2026.1 帮助文档（九）模式评估

获取Piggy_Packages还没有Piggy_Packages的同学请参考这篇帖子获取Piggy_Packages V2026.1 帮助文档一开箱即用模型评估工具MET是一种常用的对WRF预报结果进行评估的工具。今天我们来学习一项它的基础用法。有关MET的更多用法请参阅 https://dtcenter.org/software-tools/model-evaluation-tools-met/documentation首先我们回到Home目录创建一个名称为Case7的目录并进入。cd mkdir -p Case7 cd Case7准备数据将我们之前下载的观测资料和对应时刻的wrfout拷贝过来。cp ~/Case5/WRFDA/var/da/prepbufr.gdas.20260330.t00z.nr.48h . cp ~/Case1/WRF-Real/run/wrfout_d01_2026-03-30_00:00:00 .解码观测资料首先我们需要从prepbufr文件中解析并提取指定范围的观测资料。新建配置文件并编辑touch PB2NCConfig open -e PB2NCConfig粘贴以下内容// 观测时间窗口obs_window {beg -900;end 900;};// 消息类型空列表表示保留所有消息类型message_type [];// 站号空列表表示保留所有观测站station_id [];// 垂直层级范围1-255 通常涵盖所有层级level_range {beg 1;end 255;};// 数据层级类别空列表表示保留所有类别地面、强制层、显著层等level_category [];// 观测变量空列表表示保留文件中存在的所有 BUFR 变量obs_bufr_var [];// 质量标记阈值// 通常 2 表示保留较好的数据设置得越大限制越少。quality_mark_thresh 2;// 事件堆栈标志event_stack_flag TOP;// 空间遮罩mask {llpnt [ {name Region_Box;lat_thresh 10.048.0;lon_thresh 84.0136.0;} ];}// 仰角/海拔范围设置为极宽范围以包含所有高度elevation_range {beg -1000;end 100000;};// 报告类型限制均设置为空列表pb_report_type [];in_report_type [];instrument_type [];运行pb2nc:pb2nc prepbufr.gdas.20260330.t00z.nr.48h prepbufr.gdas.20260330.t00z.nc PB2NCConfig运行需要较长时间。请耐心等待。运行评估评估2米温度然后我们开始评估。我们想要检验地面2米温度。新建并编辑配置文件touch PointStatConfig open -e PointStatConfig粘贴以下内容model WRF_Model; // 输出到 STAT 文件中的模型标识符desc China_Evaluation; // 任务描述fcst {field [{ name T2; level [ Z2 ]; }];}obs {message_type [ ADPSFC, MSONET ];field [{ name TMP; level [ Z2 ]; }];}obs_window {beg -900; // 预报有效时间前 15 分钟end 900; // 预报有效时间后 15 分钟};interp {field BOTH;type [{method [ BILIN ]; // 双线性插值width [ 2 ]; // 2x2 邻域}];};output_flag {fho NONE;ctc NONE;cts BOTH; // 输出分类统计cnt BOTH; // 输出连续变量统计sl112 BOTH; // 输出标量部分和用于后续汇总mpr STAT; // 将匹配对写入 .stat 文件};然后开始检验point_stat wrfout_d01_2026-03-30_00:00:00 prepbufr.gdas.20260330.t00z.nc PointStatConfig检验后生成统计数据DEBUG 1: Output file: ./point_stat_240000L_20260330_000000V.statDEBUG 1: Output file: ./point_stat_240000L_20260330_000000V_cts.txtDEBUG 1: Output file: ./point_stat_240000L_20260330_000000V_cnt.txt我们可以通过python脚本来美化展示评估结果touch statis.py open -e statis.py粘贴下方内容#!/usr/bin/env python3 MET point_stat CNT 评估结果可视化脚本 用法: python plot_met_cnt.py [文件路径] 默认文件: point_stat_240000L_20260330_000000V_cnt.txt import sys import os try: from rich.console import Console from rich.table import Table from rich.panel import Panel from rich.columns import Columns from rich import box from rich.text import Text HAS_RICH True except ImportError: HAS_RICH False print(提示: 安装 rich 库可获得更佳显示效果 → pip install rich\n) # ───────────────────────────────────────────── # 解析函数 # ───────────────────────────────────────────── def parse_cnt_file(filepath): 解析 MET point_stat CNT 文件返回行列表每行为字典 records [] with open(filepath, r) as f: lines [l.rstrip(\n) for l in f if l.strip()] # 第一行为表头 headers lines[0].split() for line in lines[1:]: values line.split() if len(values) len(headers): values [NA] * (len(headers) - len(values)) records.append(dict(zip(headers, values))) return records def fmt(val, decimals4): 格式化数值NA → 斜体灰色rich模式或 - if val in (NA, N/A, , None): return None try: return f{float(val):.{decimals}f} except ValueError: return val # ───────────────────────────────────────────── # Rich 版本输出 # ───────────────────────────────────────────── def rich_display(records): console Console() for rec in records: mask rec.get(VX_MASK, UNKNOWN) total rec.get(TOTAL, ?) fcst_var rec.get(FCST_VAR, ?) fcst_lead rec.get(FCST_LEAD, ?) valid_time rec.get(FCST_VALID_END, ?) model rec.get(MODEL, ?) desc rec.get(DESC, ?) # ── 标题面板 ────────────────────────────────── title ( f[bold cyan]{model}[/] | {desc} | 区域: [bold yellow]{mask}[/] | f预报时效: {fcst_lead}s | 有效时间: {valid_time} | 变量: [bold]{fcst_var}[/] | 样本数: {total} ) console.print(Panel(title, expandTrue, border_stylecyan)) # ── 1. 均值与标准差 ─────────────────────────── t1 Table(title[bold]均值 标准差[/], boxbox.SIMPLE_HEAVY, show_headerTrue, header_stylebold magenta) t1.add_column(统计量, styledim, width18) t1.add_column(值, justifyright) t1.add_column(置信下限, justifyright) t1.add_column(置信上限, justifyright) rows1 [ (FBAR预报均值, FBAR, FBAR_NCL, FBAR_NCU), (OBAR观测均值, OBAR, OBAR_NCL, OBAR_NCU), (FSTDEV预报标准差,FSTDEV, FSTDEV_NCL, FSTDEV_NCU), (OSTDEV观测标准差,OSTDEV, OSTDEV_NCL, OSTDEV_NCU), ] for label, k, kl, ku in rows1: v, vl, vu rec.get(k), rec.get(kl), rec.get(ku) t1.add_row(label, fmt(v) or [dim]-[/], fmt(vl) or [dim]-[/], fmt(vu) or [dim]-[/]) console.print(t1) # ── 2. 偏差与误差 ──────────────────────────── t2 Table(title[bold]偏差 误差[/], boxbox.SIMPLE_HEAVY, show_headerTrue, header_stylebold magenta) t2.add_column(统计量, styledim, width22) t2.add_column(值, justifyright) t2.add_column(置信下限, justifyright) t2.add_column(置信上限, justifyright) rows2 [ (ME平均误差/偏差, ME, ME_NCL, ME_NCU), (MAE平均绝对误差, MAE, MAE_BCL, MAE_BCU), (MSE均方误差, MSE, MSE_BCL, MSE_BCU), (RMSE均方根误差, RMSE, RMSE_BCL, RMSE_BCU), (BCMSE偏差修正MSE, BCMSE, BCMSE_BCL, BCMSE_BCU), (ESTDEV误差标准差, ESTDEV, ESTDEV_NCL, ESTDEV_NCU), (MBIAS乘性偏差, MBIAS, MBIAS_BCL, MBIAS_BCU), ] for label, k, kl, ku in rows2: v, vl, vu rec.get(k), rec.get(kl), rec.get(ku) # 高亮ME正负 vfmt fmt(v) if vfmt and k ME: color red if float(v) 0 else blue vfmt f[{color}]{vfmt}[/] t2.add_row(label, vfmt or [dim]-[/], fmt(vl) or [dim]-[/], fmt(vu) or [dim]-[/]) console.print(t2) # ── 3. 相关系数 ────────────────────────────── t3 Table(title[bold]相关系数[/], boxbox.SIMPLE_HEAVY, show_headerTrue, header_stylebold magenta) t3.add_column(统计量, styledim, width25) t3.add_column(值, justifyright) t3.add_column(置信下限, justifyright) t3.add_column(置信上限, justifyright) rows3 [ (PR_CORRPearson相关, PR_CORR, PR_CORR_NCL, PR_CORR_NCU), (ANOM_CORR距平相关/ACC, ANOM_CORR, ANOM_CORR_NCL, ANOM_CORR_NCU), (SP_CORRSpearman等级相关, SP_CORR, None, None), (KT_CORRKendall τ相关, KT_CORR, None, None), ] for label, k, kl, ku in rows3: v rec.get(k) vl rec.get(kl) if kl else None vu rec.get(ku) if ku else None vfmt fmt(v) if vfmt: fv float(v) color green if fv 0.9 else (yellow if fv 0.7 else white) vfmt f[{color}]{vfmt}[/] t3.add_row(label, vfmt or [dim]-[/], fmt(vl) or [dim]-[/], fmt(vu) or [dim]-[/]) console.print(t3) # ── 4. 误差分位数 ──────────────────────────── t4 Table(title[bold]误差分位数[/], boxbox.SIMPLE_HEAVY, show_headerTrue, header_stylebold magenta) t4.add_column(分位数, styledim, width18) t4.add_column(值, justifyright) t4.add_column(BCL, justifyright) t4.add_column(BCU, justifyright) rows4 [ (E1010th percentile, E10, E10_BCL, E10_BCU), (E2525th percentile, E25, E25_BCL, E25_BCU), (E50中位数/50th, E50, E50_BCL, E50_BCU), (E7575th percentile, E75, E75_BCL, E75_BCU), (E9090th percentile, E90, E90_BCL, E90_BCU), (EIQR四分位距, EIQR,EIQR_BCL,EIQR_BCU), (MAD绝对中位差, MAD, MAD_BCL, MAD_BCU), ] for label, k, kl, ku in rows4: v, vl, vu rec.get(k), rec.get(kl), rec.get(ku) t4.add_row(label, fmt(v) or [dim]-[/], fmt(vl) or [dim]-[/], fmt(vu) or [dim]-[/]) console.print(t4) # ── 5. 技巧评分 ────────────────────────────── t5 Table(title[bold]技巧评分[/], boxbox.SIMPLE_HEAVY, show_headerTrue, header_stylebold magenta) t5.add_column(统计量, styledim, width25) t5.add_column(值, justifyright) t5.add_column(BCL, justifyright) t5.add_column(BCU, justifyright) rows5 [ (MSESSMSE技巧评分, MSESS, MSESS_BCL, MSESS_BCU), (ME2平均误差平方, ME2, ME2_BCL, ME2_BCU), (SI散度指数, SI, SI_BCL, SI_BCU), (RMSFA预报均方根距平, RMSFA, RMSFA_BCL, RMSFA_BCU), (RMSOA观测均方根距平, RMSOA, RMSOA_BCL, RMSOA_BCU), ] for label, k, kl, ku in rows5: v, vl, vu rec.get(k), rec.get(kl), rec.get(ku) t5.add_row(label, fmt(v) or [dim]-[/], fmt(vl) or [dim]-[/], fmt(vu) or [dim]-[/]) console.print(t5) console.print() # 行间距 # ───────────────────────────────────────────── # 纯文本版本无 rich # ───────────────────────────────────────────── def plain_display(records): def row(label, val, loNone, hiNone, w28): val_s f{float(val):.4f} if val not in (None, NA) else - lo_s f{float(lo):.4f} if lo not in (None, NA) else - hi_s f{float(hi):.4f} if hi not in (None, NA) else - print(f {label:{w}} {val_s:10} {lo_s:10} {hi_s:10}) def header(title): print(f\n{─*65}) print(f {title}) print(f{─*65}) print(f {统计量:28} {值:10} {置信/下限:10} {置信/上限:10}) print(f {-*28} {-*10} {-*10} {-*10}) for rec in records: print(\n ═*65) print(f 模型: {rec.get(MODEL)} 区域: {rec.get(VX_MASK)}) print(f 变量: {rec.get(FCST_VAR)} 时效: {rec.get(FCST_LEAD)}s 样本数: {rec.get(TOTAL)}) print(═*65) header(均值 标准差) for lbl, k, kl, ku in [ (FBAR 预报均值, FBAR, FBAR_NCL, FBAR_NCU), (OBAR 观测均值, OBAR, OBAR_NCL, OBAR_NCU), (FSTDEV 预报标准差,FSTDEV,FSTDEV_NCL,FSTDEV_NCU), (OSTDEV 观测标准差,OSTDEV,OSTDEV_NCL,OSTDEV_NCU), ]: row(lbl, rec.get(k), rec.get(kl), rec.get(ku)) header(偏差 误差) for lbl, k, kl, ku in [ (ME 平均误差, ME, ME_NCL, ME_NCU), (MAE 平均绝对误差,MAE, MAE_BCL, MAE_BCU), (MSE 均方误差, MSE, MSE_BCL, MSE_BCU), (RMSE 均方根误差,RMSE, RMSE_BCL, RMSE_BCU), (BCMSE 偏差修正MSE,BCMSE,BCMSE_BCL,BCMSE_BCU), (MBIAS 乘性偏差, MBIAS,MBIAS_BCL,MBIAS_BCU), ]: row(lbl, rec.get(k), rec.get(kl), rec.get(ku)) header(相关系数) for lbl, k, kl, ku in [ (PR_CORR Pearson相关, PR_CORR, PR_CORR_NCL, PR_CORR_NCU), (ANOM_CORR 距平相关ACC, ANOM_CORR, ANOM_CORR_NCL, ANOM_CORR_NCU), (SP_CORR Spearman, SP_CORR, None, None), (KT_CORR Kendall τ, KT_CORR, None, None), ]: row(lbl, rec.get(k), rec.get(kl) if kl else None, rec.get(ku) if ku else None) header(误差分位数) for lbl, k, kl, ku in [ (E10 10th percentile,E10,E10_BCL,E10_BCU), (E25 25th percentile,E25,E25_BCL,E25_BCU), (E50 中位数, E50,E50_BCL,E50_BCU), (E75 75th percentile,E75,E75_BCL,E75_BCU), (E90 90th percentile,E90,E90_BCL,E90_BCU), (EIQR 四分位距, EIQR,EIQR_BCL,EIQR_BCU), (MAD 绝对中位差, MAD,MAD_BCL,MAD_BCU), ]: row(lbl, rec.get(k), rec.get(kl), rec.get(ku)) header(技巧评分) for lbl, k, kl, ku in [ (MSESS MSE技巧评分,MSESS,MSESS_BCL,MSESS_BCU), (ME2 平均误差平方,ME2, ME2_BCL, ME2_BCU), (SI 散度指数, SI, SI_BCL, SI_BCU), ]: row(lbl, rec.get(k), rec.get(kl), rec.get(ku)) print() # ───────────────────────────────────────────── # 主程序 # ───────────────────────────────────────────── def main(): filepath sys.argv[1] if len(sys.argv) 1 else point_stat_240000L_20260330_000000V_cnt.txt if not os.path.exists(filepath): print(f[错误] 文件不存在: {filepath}) sys.exit(1) records parse_cnt_file(filepath) print(f\n共解析 {len(records)} 条 CNT 评估记录\n) if HAS_RICH: rich_display(records) else: plain_display(records) if __name__ __main__: main()运行python statis.py point_stat_240000L_20260330_000000V_cnt.txt评估500hPa温度编辑配置文件open -e PointStatConfig粘贴以下内容