Prepbufr introduction and decoding

奇妙_光明 2021-09-16

展开全文

Prepbufr 格式数据介绍及解码

Prepbufr的背景及介绍

PrepBUFR是NCEP在WMO BUFR格式的基础上，对码表内容进行扩充的数据格式，不仅包含观测数据，还包含观测误差、背景场等信息。该格式目前是NCEP业务模式与同化系统（GSI）所使用的常规观测资料数据格式。RAP/HRRR同样也使用了这套观测数据格式。国内也有引进RAP/HRRR区域数值预报系统，如武汉暴雨研究所，福建省局等。中国气象局数值预报中心的GRAPES模式中观测是文本文件，没有直接使用Prepbufr。

NCEP Prepbufr介绍：https://www.emc.ncep./mmb/data_processing/prepbufr.doc/document.htm
GSI同化系统中Prepbufr介绍：
https:///sites/default/files/community-code/gsi/downloads/BUFR/presentations/gsi-tutorial-bufr/L6-06292011-BUFRandPrepBUFR-Ruifang_Li.pdf
http://www./sites/default/files/community-code/gsi/downloads/BUFR/BUFR_PrepBUFR_User_Guide_v1.pdf

NCEP数据处理系统流程

Prepbufr由美国NCEP的质量控制系统（obsproc）生成，而质控控制系统也是结合数值模式来运行，对比和背景场的差异。当然也可以单独运行，但质量控制效果就差很多，背景场本质代表过去的观测对当前的影响。总体分为以下步骤：

PREPOBS_PREPDATA
总体质量控制
GLERLADJ
对水体周边观测的调整
PREVENTS
在气候模式分析准备观测资料时运行
比照背景场，对地面气压进行一些粗略的质量控制检查
CQCBUFR
对探空观测的质控
PROFCQC
对飞机观测的质控
PREPACQC
对风廓线仪观测的质控
OIQCBUFR
卡住值检查
使用所有的观测来执行独立的质控检查

obsproc对不同尺度的分析质量控制过程不一样：如全球预报系统GFS，RAP/HRRR，以及中尺度实况分析RTMA。

GFS全球预报系统：https://www.emc.ncep./mmb/data_processing/prepbufr.doc/table_2.htm

RAP/HRRR区域模式：https://www.emc.ncep./mmb/data_processing/prepbufr.doc/table_5.htm

NAM北美区域模式：https://www.emc.ncep./mmb/data_processing/prepbufr.doc/table_4.htm

RTMA中尺度实况分析：https://www.emc.ncep./mmb/data_processing/prepbufr.doc/table_19.htm
更多信息请访问：https://www.emc.ncep./mmb/data_processing/prepbufr.doc/

质控码表： https://www.emc.ncep./mmb/data_processing/prepbufr.doc/table_7.htm
其中质控码小于2的认为是质量比较好的观测，比如WRFDA中质控码设置的就是2。

关于NCEP质量控制的更多信息请参见上面NCEP的Prepbufr介绍。

NCEP obsproc 源代码地址: https://www.nco.ncep./pmb/codes/nwprod/

安装bufrlib及其他解码程序

NCEP对bufr的编码和解码用bufrlib库，欧洲中心的eccodes也可以对bufr进行操作，但本人没尝试成功。

基于bufrlib库，目前有很多现成的解码程序使用，如：

GSI同化系统中的Fortan解码程序: GSI/util/bufr_tools
NCEP ADP BUFR 转little r程序
美国卫星联合同化中心（JCSDA）开发的Python库：pyncepbufr

1. 编译bufrlib

下载bufrlib源码：https://emc.ncep./emc/pages/infrastructure/bufrlib.php

tar xvf BUFRLIB_v11-3-0.tar


 CC=gcc
 FC=gfortan

 -O3 -DUNDERSCORE -w  -c `./getdefflags_C.sh` *.c

 -O3 -DUNDERSCORE -fno-second-underscore -w  -c `./getdefflags_F.sh` modv*.F moda*.F `ls -1 *.F *.f | grep -v `

ar crv bufrlib.a *.o

2. 编译GSI中bufr工具

下载源码
https://github.com/NOAA-EMC/GSI/tree/master/util/bufr_tools

编译，以prepbufr_decode_all_evn.f90为例

1	gfortran prepbufr_decode_all_evn.f90 -o prepbufr_decode_all_evn.exe /app/bufrlib/bufrlib.a

3. 编译BUFR转`little r`程序

下载源码：https://github.com/NCAR/rda-bufr-decode-ADPsfc-littler

本人在解码最新的prepbufr时有点问题，需要在src/bufr_sfc2ob.f中第267行再加一个判断，改为：
if(M1(1:1) == 'm’ .or. M1(2:2) == 'm’) then

编译：修改install/install.sh第22行，改为bufrlib.a所在的路径。然后执行install/install.sh
编译成功后，在exe下生成如下可执行文件：bufr_sfc2ob.x bufr_ship2ob.x dumpbufr.x runob2lit_imd_obs.x

3. 安装`py-ncepbufr`

源码地址：https://github.com/JCSDA-internal/py-ncepbufr
git仓库有详细的安装说明，就不再赘述

解码及绘图

下载Prepbufr数据，实时的可在nomads（NOAA模式存档和分发系统）下载：https://nomads.ncep./

历史的可在NCAR RDA网站可下载：https://rda./datasets/ds337.0/

下面以GFS的观测数据为例做解码：https://nomads.ncep./pub/data/nccf/com/gfs/prod/gfs.20210130/00/gfs.t00z.prepbufr.nr
里面包含了，地面站，探空站，浮标船舶，GPS以及卫星风等。
单独的bufr文件如：

文件名	观测类型
gfs.t00z.adpsfc.tm00.bufr_d.nr	地面数据
gfs.t00z.adpupa.tm00.bufr_d	探空数据
gfs.t00z.aircar.tm00.bufr_d.nr	飞机报
gfs.t00z.gpsipw.tm00.bufr_d.nr	GPS
gfs.t00z.satwnd.tm00.bufr_d	卫星风
还有很多其他的观测类型数据，就不一一列出了。

1. 下载`Prepbufr`数据

lftp -e 

mv gfs.t00z.prepbufr.nr gfs.2021013000.prepbufr

2. 使用`prepbufr_decode_all_evn.exe`

ln -s gfs.2021013000.prepbufr prepbufr

./prepbufr_decode_all_evn.exe > gfs.2021013000_qc.txt

gfs.2021013000_qc.txt为解码后的结果，以北京站54511为例：

解码结果，54511地面站

其中，obs为观测值，qcf为质控码。

3. 使用`rda-bufr-decode-ADPsfc`程序包

1
2
3


./exe/bufr_sfc2ob.x gfs.t00z.adpsfc.tm00.bufr_d.nr 2021013000

Surface2021011400.obs为解码后的结果。

 Surface2021013000.obs > files.txt

./exe/runob2lit_imd_obs.x files.txt 2021013000

SURFACE_OBS:2021013000为转换后的little r格式数据。

4. 使用`py-ncepbufr`

读取数据和质控码程序如下：

 __future__  print_function
 argparse
 datetime  datetime, timedelta
 ncepbufr
 os
 sys

parser = argparse.ArgumentParser()

parser.add_argument(, , =, required=)

args = parser.parse_args()

hdstr=
obstr=
qcstr=
oestr=



bufr = ncepbufr.(args.file)
 bufr.advance() == : 
    ref_time = datetime.strptime((bufr.msg_date), )
     bufr.load_subset() == : 
        hdr = bufr.read_subset(hdstr).squeeze()
:
          station_id = hdr[].tostring()
         Exception  e:
          station_id = 
        obs = bufr.read_subset(obstr)
        nlevs = obs.shape[-]
        oer = bufr.read_subset(oestr)
        qcf = bufr.read_subset(qcstr)
        time = ref_time + timedelta(seconds=hdr[]*)
        print( % (station_id,hdr[],hdr[],time,(hdr[]),nlevs))
        print(obs)
bufr.close()

运行解码程序：