Valar Dohaeris

Prepbufr 格式数据介绍及解码

Prepbufr的背景及介绍

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

NCEP Prepbufr介绍：https://www.emc.ncep.noaa.gov/mmb/data_processing/prepbufr.doc/document.htm
GSI同化系统中Prepbufr介绍：
https://dtcenter.org/sites/default/files/community-code/gsi/downloads/BUFR/presentations/gsi-tutorial-bufr/L6-06292011-BUFRandPrepBUFR-Ruifang_Li.pdf
http://www.dtcenter.org/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.noaa.gov/mmb/data_processing/prepbufr.doc/table_2.htm

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

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

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

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

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

NCEP obsproc 源代码地址: https://www.nco.ncep.noaa.gov/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.noaa.gov/emc/pages/infrastructure/bufrlib.php

# 解压
tar xvf BUFRLIB_v11-3-0.tar
# 编译说明在解压后的README_BUFRLIB中
# 首先设置编译器
export CC=gcc
export FC=gfortan
# 编译C程序
$CC -O3 -DUNDERSCORE -w  -c `./getdefflags_C.sh` *.c
# 编译Fortran
$FC -O3 -DUNDERSCORE -fno-second-underscore -w  -c `./getdefflags_F.sh` modv*.F moda*.F `ls -1 *.F *.f | grep -v "mod[av]_"`
# 打包为静态库
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.noaa.gov/

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

下面以GFS的观测数据为例做解码：https://nomads.ncep.noaa.gov/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并行下载
lftp -e 'pget -c -n 20 "https://nomads.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.20210130/00/gfs.t00z.prepbufr.nr"'
# 重命名
mv gfs.t00z.prepbufr.nr gfs.2021013000.prepbufr

2. 使用`prepbufr_decode_all_evn.exe`

# 链接
ln -s gfs.2021013000.prepbufr prepbufr
# 运行解码程序，这里也可以直接改prepbufr_decode_all_evn.f90进行个性化输入输出
./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为解码后的结果。

# 将输出的文件名写入files.txt
echo Surface2021013000.obs > files.txt
# 转little r
./exe/runob2lit_imd_obs.x files.txt 2021013000

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

4. 使用`py-ncepbufr`

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

#!/bin/env python
from __future__ import print_function
import argparse
from datetime import datetime, timedelta
import ncepbufr
import os
import sys

parser = argparse.ArgumentParser()

parser.add_argument("--file", "-f", type=str, required=True)

args = parser.parse_args()

hdstr='SID XOB YOB DHR TYP ELV SAID T29'
obstr='POB QOB TOB ZOB UOB VOB PWO MXGS HOVI CAT PRSS TDO PMO'
qcstr='PQM QQM TQM ZQM WQM NUL PWQ PMQ'
oestr='POE QOE TOE NUL WOE NUL PWE'

# read prepbufr file.

bufr = ncepbufr.open(args.file)
while bufr.advance() == 0: # loop over messages.
    ref_time = datetime.strptime(str(bufr.msg_date), "%Y%m%d%H")
    while bufr.load_subset() == 0: # loop over subsets in message.
        hdr = bufr.read_subset(hdstr).squeeze()
	try:
          station_id = hdr[0].tostring()
        except Exception as e:
          station_id = "unknown"
        obs = bufr.read_subset(obstr)
        nlevs = obs.shape[-1]
        oer = bufr.read_subset(oestr)
        qcf = bufr.read_subset(qcstr)
        time = ref_time + timedelta(seconds=hdr[3]*3600)
        print("%s %8.2f %8.2f %s %3.3d %2.2d" % (station_id,hdr[1],hdr[2],time,int(hdr[4]),nlevs))
        print(obs)
bufr.close()

运行解码程序：

1	python read_prepbufr.py -f gfs.2021013000.prepbufr

通过修改以上程序，可以获得个性化输出，比如输出POB QOB TOB即为气压，湿度，温度等。

此外，py-ncepbufr也提供了一些小工具，如：py-ncepbufr/utils/{prepbufr2nc,nc2prepbufr}等。

4. 使用`MET`画观测分布

MET是模式检验工具，全称The Model Evaluation Tools，the Developmental Testbed Center所开发。

在线文档：https://dtcenter.github.io/MET/latest/Users_Guide/index.html

MET有现成的docker镜像，可pull下来直接食用。还有相应的Python脚本程序metplus。

1	docker pull dtcenter/met

Prepbufr转为netcdf

1	pb2nc gfs.2021013000.prepbufr gdas_2021013000.nc PB2NCConfig

画地面站分布

1	plot_point_obs -msg_typ ADPSFC gdas_2021013000.nc gdas_2021013000_adpsfc.ps

全球地面站分布

画探空站分布

1	plot_point_obs -msg_typ ADPUPA gdas_2021013000.nc gdas_2021013000_adpupa.ps

全球探空站分布

写在最后：

以上解码Prepbufr我们可以获得完整的全球GTS交换数据。对于国内常常拿不到的常规观测或探空，通过解码Prepbufr获取也不失为一种“正当途径”。

背景

最近重新整理之前的分析结果，发现一些有意思的图，来分享出来。2020年美国飓风 Laura 临近登陆时，NOAA派出飞机投放探空仪观测，很快推特上便有了相应的T-lnP图，如下图。

而我国在2020年，也进行了某某计划，在台风森拉克过境后，即在台风外围投放了数十枚探空仪，获取到了非常密集的温压湿风探空曲线。这是当时的卫星云图(来源：韩国GK2A卫星)，可以看到当时台风主体已经进入越南境内，而释放探空仪的位置在海南岛东南近海。

真彩色云图

红外增强云图

读取数据

由于数据是excel的表格，我们用pandas来读取表单。

from datetime import datetime,timedelta
import pandas as pd
f = '17_54_4.xls'
# 将时间解析为datetime类
dateparse = lambda x: datetime.strptime('2020-08-02 ' + x, '%Y-%m-%d %H:%M:%S')
df = pd.read_excel(f,sheet_name=0,parse_dates=['时间'], date_parser=dateparse)
# 时间         纬度          经度     位势高度     温度     湿度     气压  风速  风向
#                                     m           degC     %       hPa   m/s  deg   
# convert BJT -> UTC
df['时间'] = df['时间'] - timedelta(hours=8)
df[df == -9999] = np.nan
df = df[::-1]

读取结果

绘制T-lnP图

用NCAR开发的Metpy包来画T-lnP图，官方链接为：https://unidata.github.io/MetPy/latest/
官方提供了很多例子，包括画常规天气图和T-lnP，以及计算物理量、诊断量。链接：https://unidata.github.io/MetPy/latest/examples/index.html

根据Simple_Sounding例子，来画我们需要的最简单的T-lnP。注意：美国人用的是斜T-lnP，即旋转30度，而国内是T-lnP，需要改成0度。对应的函数接口为：
https://unidata.github.io/MetPy/latest/api/generated/metpy.plots.SkewT.html#metpy.plots.SkewT

SkewT的函数说明

下面为具体绘图代码

import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots import SkewT, Hodograph
from metpy.units import units
import metpy.interpolate as mpi

# 转换为Metpy的数据类
p = df['气压'].values * units.hPa
T = df['温度'].values * units.degC
rh = df['湿度'].values / 100
Td = mpcalc.dewpoint_rh(T, rh)
wsp = df['风速'].values * units('m/s')
wdr = df['风向'].values * units.degrees
u, v = mpcalc.wind_components(wsp, wdr)

# 开一个画布
fig = plt.figure(figsize=(12, 12))
skew = SkewT(fig, rotation=0)

# 画温度、露点和风
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')
skew.plot_barbs(p[::5], u[::5], v[::5])
skew.ax.set_ylim(1000, 250)
skew.ax.set_xlim(-40, 30)

# 计算对流凝结高度，p[0]为最低层气压
lcl_pressure, lcl_temperature = mpcalc.lcl(p[0], T[0], Td[0])
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')

# 画状态曲线
prof = mpcalc.parcel_profile(p, T[0], Td[0]).to('degC')
skew.plot(p, prof, 'k', linewidth=2)

# 画对流有效位能和抑制能
skew.shade_cin(p, T, prof)
skew.shade_cape(p, T, prof)

# 画0度温度线，图中紫色虚线
skew.ax.axvline(0, color='c', linestyle='--', linewidth=2)

# 干绝热线
skew.plot_dry_adiabats()
# 湿绝热线
skew.plot_moist_adiabats()
# 等比湿线
skew.plot_mixing_lines()

T-lnP图, 世界时: 2020-08-02T10

可以看出台风过后CAPE还是比较大。好了，这期就到这了。

Valar Dohaeris

Prepbufr introduction and decoding

Prepbufr 格式数据介绍及解码

Prepbufr的背景及介绍

安装bufrlib及其他解码程序

1. 编译bufrlib

2. 编译GSI中bufr工具

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

3. 安装`py-ncepbufr`

解码及绘图

1. 下载`Prepbufr`数据

2. 使用`prepbufr_decode_all_evn.exe`

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

4. 使用`py-ncepbufr`

4. 使用`MET`画观测分布

Python processing sounding data and plotting

背景

读取数据

绘制T-lnP图

Container technology in meteorological applications

容器技术

Hello World

Quick Start

Create a new post

Run server

Generate static files

Deploy to remote sites

Prepbufr 格式数据介绍及解码

Prepbufr的背景及介绍

安装bufrlib及其他解码程序

1. 编译bufrlib

2. 编译GSI中bufr工具

3. 编译BUFR转little r程序

3. 安装py-ncepbufr

解码及绘图

1. 下载Prepbufr数据

2. 使用prepbufr_decode_all_evn.exe

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

4. 使用py-ncepbufr

4. 使用MET画观测分布

背景

读取数据

绘制T-lnP图

容器技术

Quick Start

Create a new post

Run server

Generate static files

Deploy to remote sites

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

3. 安装`py-ncepbufr`

1. 下载`Prepbufr`数据

2. 使用`prepbufr_decode_all_evn.exe`

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

4. 使用`py-ncepbufr`

4. 使用`MET`画观测分布