使用wgrib2实现GRIB2文件要素场抽取

meteorology/data
目录

本文介绍如何使用 wgrib2 和 grep 从 GRIB2 文件中抽取指定要素场生成新的 GRIB2 文件。

wgrib2 是美国 NCEP 开发的 GRIB2 格式文件处理工具。

需求

CEMC 正在建立 CMA-MESO 1KM 逐小时循环预报系统，并面向 2021chengdu 开发预报产品。 因为逐小时循环系统一天运行 24 次，所以一个预报时次的所有产品需要在一个小时内完成传输。 前期已生成一套 GRIB2 产品（grib2-orig），每个文件 340 MB，一个时次 28 个文件一共 9.5 GB。 假设下载带宽为 3MB/s，那么下载一个时次全部数据需要 56 分钟，对于每个小时都运行的系统来说基本没有容错空间。 考虑到网络传输带宽有限，需要对数据进行精简，提供尽可能少的要素场，减少单个数据文件的大小。 这就需要从已有 GRIB2 产品文件中抽取指定要素场生成新的 GRIB2 产品文件（grib2-bccd）。

grib2-orig 产品示例 rmf.hgra.2023042300006.grb2 文件的 wgrib2 输出如下所示，每行代表一个要素场，一共 256 个要素场。

1:0:d=2023042300:APCP:surface:0-6 hour acc fcst:
2:2462709:d=2023042300:TMP:surface:6 hour fcst:
3:4321648:d=2023042300:HGT:surface:6 hour fcst:
4:6400691:d=2023042300:SPFH:2 m above ground:6 hour fcst:
5:7338630:d=2023042300:TMP:2 m above ground:6 hour fcst:
...
251:358186619:d=2023042300:var discipline=0 center=38 local_table=1 parmcat=7 parm=225:surface:6 hour fcst:
252:358387277:d=2023042300:PWAT:entire atmosphere:6 hour fcst:
253:360263004:d=2023042300:VIS:surface:6 hour fcst:
254:361389756:d=2023042300:GUST:10 m above ground:6 hour fcst:
255:363531777:d=2023042300:CDCTOP:surface:6 hour fcst:
256:365300689:d=2023042300:CDCB:surface:6 hour fcst:

根据要素表格对数据进行抽取，要素表格如下所示：

序号要素要素名单位层次数层次清单
1东西风UGRDm/s21100, 150,200, 250,300,350, 400, 450, 500, 550, 600, 650, 700, 750, 800, 850, 900, 925, 950, 975, 1000
2南北风VGRDm/s21同上
310米UUGRDm/s110 m above ground
410米VVGRDm/s110 m above ground
………………
2210米阵风风速GUSTm/s110m
23等高面u风UGRDm930，50，70，100，120，140，160，180，200（等高面）
24等高面v风VGRDm930，50，70，100，120，140，160，180，200（等高面）

抽取后的 grib2-bccd 产品示例 bccd_rmf.hgra.2023042300006.grb2 文件的 wgrib2 输出如下所示，一共 80 个要素场。 文件大小为 134 MB，比源文件减少一半以上，有效缩短传输时间。

1:0:d=2023042300:UGRD:1000 mb:6 hour fcst:
2:1994227:d=2023042300:UGRD:975 mb:6 hour fcst:
3:3990785:d=2023042300:UGRD:950 mb:6 hour fcst:
4:5989798:d=2023042300:UGRD:925 mb:6 hour fcst:
5:7995377:d=2023042300:UGRD:900 mb:6 hour fcst:
6:10000325:d=2023042300:UGRD:850 mb:6 hour fcst:
7:11987388:d=2023042300:UGRD:800 mb:6 hour fcst:
8:13938722:d=2023042300:UGRD:750 mb:6 hour fcst:
9:15857046:d=2023042300:UGRD:700 mb:6 hour fcst:
10:17748465:d=2023042300:UGRD:650 mb:6 hour fcst:
11:19580122:d=2023042300:UGRD:600 mb:6 hour fcst:
12:21335505:d=2023042300:UGRD:550 mb:6 hour fcst:
13:22965168:d=2023042300:UGRD:500 mb:6 hour fcst:
14:24502364:d=2023042300:UGRD:450 mb:6 hour fcst:
15:25992306:d=2023042300:UGRD:400 mb:6 hour fcst:
16:27444467:d=2023042300:UGRD:350 mb:6 hour fcst:
17:28905525:d=2023042300:UGRD:300 mb:6 hour fcst:
18:30367368:d=2023042300:UGRD:250 mb:6 hour fcst:
19:31858071:d=2023042300:UGRD:200 mb:6 hour fcst:
20:33316369:d=2023042300:UGRD:150 mb:6 hour fcst:
21:34598102:d=2023042300:UGRD:100 mb:6 hour fcst:
22:36104901:d=2023042300:VGRD:1000 mb:6 hour fcst:
23:38080669:d=2023042300:VGRD:975 mb:6 hour fcst:
24:40059920:d=2023042300:VGRD:950 mb:6 hour fcst:
25:42037693:d=2023042300:VGRD:925 mb:6 hour fcst:
26:44022552:d=2023042300:VGRD:900 mb:6 hour fcst:
27:46013289:d=2023042300:VGRD:850 mb:6 hour fcst:
28:47982656:d=2023042300:VGRD:800 mb:6 hour fcst:
29:49918054:d=2023042300:VGRD:750 mb:6 hour fcst:
30:51805225:d=2023042300:VGRD:700 mb:6 hour fcst:
31:53650037:d=2023042300:VGRD:650 mb:6 hour fcst:
32:55426405:d=2023042300:VGRD:600 mb:6 hour fcst:
33:57125123:d=2023042300:VGRD:550 mb:6 hour fcst:
34:58724694:d=2023042300:VGRD:500 mb:6 hour fcst:
35:60248661:d=2023042300:VGRD:450 mb:6 hour fcst:
36:61719140:d=2023042300:VGRD:400 mb:6 hour fcst:
37:63177060:d=2023042300:VGRD:350 mb:6 hour fcst:
38:64630479:d=2023042300:VGRD:300 mb:6 hour fcst:
39:66057813:d=2023042300:VGRD:250 mb:6 hour fcst:
40:67485994:d=2023042300:VGRD:200 mb:6 hour fcst:
41:68873318:d=2023042300:VGRD:150 mb:6 hour fcst:
42:70198025:d=2023042300:VGRD:100 mb:6 hour fcst:
43:71562029:d=2023042300:UGRD:10 m above ground:6 hour fcst:
44:73573804:d=2023042300:VGRD:10 m above ground:6 hour fcst:
45:75570687:d=2023042300:TMP:2 m above ground:6 hour fcst:
46:77697591:d=2023042300:TMP:surface:6 hour fcst:
47:79556530:d=2023042300:PRMSL:mean sea level:6 hour fcst:
48:81406893:d=2023042300:PRES:surface:6 hour fcst:
49:83475101:d=2023042300:RH:2 m above ground:6 hour fcst:
50:85272968:d=2023042300:APCP:surface:0-6 hour acc fcst:
51:87735677:d=2023042300:HGT:surface:6 hour fcst:
52:89814720:d=2023042300:CAPE:surface:6 hour fcst:
53:90036962:d=2023042300:KX:mean sea level:6 hour fcst:
54:91403569:d=2023042300:var discipline=0 center=38 local_table=1 parmcat=16 parm=224:surface:6 hour fcst:
55:91403750:d=2023042300:BLI:surface:6 hour fcst:
56:93585687:d=2023042300:var discipline=0 center=38 local_table=1 parmcat=1 parm=239:surface:6 hour fcst:
57:93792981:d=2023042300:DPT:2 m above ground:6 hour fcst:
58:95793163:d=2023042300:UPHL:5000-2000 m above ground:0-6 hour max fcst:
59:96411591:d=2023042300:PWAT:entire atmosphere:6 hour fcst:
60:98287318:d=2023042300:var discipline=0 center=38 local_table=1 parmcat=2 parm=237:600-0 m above ground:6 hour fcst:
61:100162117:d=2023042300:VIS:surface:6 hour fcst:
62:101288869:d=2023042300:GUST:10 m above ground:6 hour fcst:
63:103430890:d=2023042300:UGRD:30 m above ground:6 hour fcst:
64:105516703:d=2023042300:UGRD:50 m above ground:6 hour fcst:
65:107602346:d=2023042300:UGRD:70 m above ground:6 hour fcst:
66:109680759:d=2023042300:UGRD:100 m above ground:6 hour fcst:
67:111744612:d=2023042300:UGRD:120 m above ground:6 hour fcst:
68:113798807:d=2023042300:UGRD:140 m above ground:6 hour fcst:
69:115842656:d=2023042300:UGRD:160 m above ground:6 hour fcst:
70:117877950:d=2023042300:UGRD:180 m above ground:6 hour fcst:
71:119904468:d=2023042300:UGRD:200 m above ground:6 hour fcst:
72:121923991:d=2023042300:VGRD:30 m above ground:6 hour fcst:
73:124000886:d=2023042300:VGRD:50 m above ground:6 hour fcst:
74:126076287:d=2023042300:VGRD:70 m above ground:6 hour fcst:
75:128144424:d=2023042300:VGRD:100 m above ground:6 hour fcst:
76:130197248:d=2023042300:VGRD:120 m above ground:6 hour fcst:
77:132241170:d=2023042300:VGRD:140 m above ground:6 hour fcst:
78:134274947:d=2023042300:VGRD:160 m above ground:6 hour fcst:
79:136300472:d=2023042300:VGRD:180 m above ground:6 hour fcst:
80:138318328:d=2023042300:VGRD:200 m above ground:6 hour fcst:

方法

使用 wgrib2 从 GRIB2 文件中抽取要素场可以通过以下三个步骤实现：

  1. 获取 GRIB2 文件的要素索引信息；
  2. 对索引信息进行过滤，选择需要的要素场，生成筛选后的索引信息；
  3. 根据筛选后的索引信息从 GRIB2 文件中抽取要素场。

使用 wgrib2 从 GRIB2 文件中抽取要素场生成新文件

实现

下面具体介绍三个步骤的实现代码。

生成要素清单

将 wgrib2 输出的要素清单打印到要素清单文本文件 wgrib2_list.txt 中

wgrib2 rmf.hgra.2023042300006.grb2 > wgrib2_list.txt

生成要素清单文件 wgrib2_list.txt 内容如下：

1:0:d=2023042300:APCP:surface:0-6 hour acc fcst:
2:2462709:d=2023042300:TMP:surface:6 hour fcst:
3:4321648:d=2023042300:HGT:surface:6 hour fcst:
4:6400691:d=2023042300:SPFH:2 m above ground:6 hour fcst:
5:7338630:d=2023042300:TMP:2 m above ground:6 hour fcst:
...

筛选要素

使用 grep 命令将符合条件的要素场条目写入到文件 wgrib2_selected.txt 中。 按照要素表格顺序逐一筛选要素场。

for var in  ":UGRD:[0-9]* mb:" \
            ":VGRD:[0-9]* mb:" \
            ":UGRD:10 m above ground:" \
            ":VGRD:10 m above ground:" \
            ":TMP:2 m above ground:"             \
            ":TMP:surface:"      \
            ":PRMSL:"      \
            ":PRES:"      \
            ":RH:2 m above ground:"      \
            ":APCP:surface:"      \
            ":HGT:surface:"      \
            ":CAPE:"      \
            ":KX:"      \
            ":var discipline=0 center=38 local_table=1 parmcat=16 parm=224:surface:"   \
            ":BLI:"      \
            ":var discipline=0 center=38 local_table=1 parmcat=1 parm=239:surface:"    \
            ":DPT:2 m above ground:"   \
            ":UPHL:5000-2000 m above ground:"  \
            "PWAT"   \
            ":var discipline=0 center=38 local_table=1 parmcat=2 parm=237:600-0 m above ground:" \
            ":VIS:"  \
            ":GUST:10 m above ground:"  \
;do
    grep -E "${var}" wgrib2_list.txt >> wgrib2_selected.txt
done

level_list="30 50 70 100 120 140 160 180 200"

for level in ${level_list}; 
do
    grep -E ":UGRD:${level} m above ground:" wgrib2_list.txt >> wgrib2_selected.txt
done

for level in ${level_list}; 
do
    grep -E ":VGRD:${level} m above ground:" wgrib2_list.txt >> wgrib2_selected.txt
done

筛选后生成的要素清单文件 wgrib2_selected.txt 内容如下，一共 80 个要素。

80:121764774:d=2023042300:UGRD:1000 mb:6 hour fcst:
81:123759001:d=2023042300:UGRD:975 mb:6 hour fcst:
...
36:63555002:d=2023042300:VGRD:180 m above ground:6 hour fcst:
37:65572858:d=2023042300:VGRD:200 m above ground:6 hour fcst:

清单文件里的要素场顺序决定生成的 GRIB2 文件的要素场顺序，以上代码按照要素表格排序。 如果想要与输入数据 grib2-orig 保持一样的顺序，可以对要素清单每一行开头的数字（即要素场序号）进行排序。

sort -t ';' -n wgrib2_selected.txt

1:0:d=2023042300:APCP:surface:0-6 hour acc fcst:
2:2462709:d=2023042300:TMP:surface:6 hour fcst:
3:4321648:d=2023042300:HGT:surface:6 hour fcst:
5:7338630:d=2023042300:TMP:2 m above ground:6 hour fcst:
6:9465534:d=2023042300:UGRD:10 m above ground:6 hour fcst:
7:11477309:d=2023042300:VGRD:10 m above ground:6 hour fcst:
13:21157567:d=2023042300:var discipline=0 center=38 local_table=1 parmcat=16 parm=224:surface:6 hour fcst:
...

抽取要素场

使用 wgrib2 根据要素清单抽取要素场，-i_file 指定要素清单文件，-grib 设置输出 GRIB 文件的路径。

wgrib2 rmf.hgra.2023042300006.grb2 \
  -i_file wgrib2_selected.txt \
  -grib bccd_rmf.hgra.2023042300006.grb2

总结

生成索引-筛选索引-抽取要素 是处理 GRIB2 文件的常用方法。 无论是 NCEP 的 wgrib2 还是 CEMC 的 gribpost 都支持基于清单的要素场抽取。

wgrib2 的 -match 参数也可以抽取一个或一类要素场。 但在抽取多个要素场时，应尽量避免多次执行 wgrib2 -match ... -append ... 命令逐步生成 GRIB2 文件的方式，减少 wgrib2 命令执行次数，避免重复遍历 GRIB2 文件。 详细说明请参考 wgrib2 -i 参数文档。

笔者正在使用 cemc-oper/reki 库实现类似的功能，后面会介绍如何使用 reki 库实现基于清单的要素场抽取。

参考

wgrib2 官方网站

wgrib2 -i/-i_file 参数文档

参考博文：

Dask应用：并行抽取站点数据

Dask应用：从多个文件中抽取GRIB2要素场

Dask应用：批量转换GRIB2文件中的要素场