利用python解码GDAS的prepbufr格式数据并画出数据来源站点分布

当我们使用GSI同化系统来对大气数据进行同化时，都要使用bufr/prepbufr格式的数据，本文以解码GDAS的prepbufr格式数据来可视化所用的数据信息。prepbufr数据的解码编译一般是使用的Fortran语言，但python也有相应的包可以使用
可以参考github库ncepbufr

GDAS以及prepbufr数据简介

1. GDAS介绍
   GDAS（Global Data Assimilation System，全球资料同化系统）GDAS是由美国国家环境预报中心（NCEP）开发和维护的一个全球数值天气预报系统。GDAS的主要功能是将全球观测数据（如卫星、雷达、地面站、气球探测等）与数值天气预报模型结合，通过同化技术生成一个最优的初始条件，这些初始条件将用于驱动数值天气预报模型，以提高预报的准确性。
   GDAS通常会以6小时为一个周期进行同化和更新，生成全球范围内的分析场，这些分析场可以用来初始化数值天气预报模型。GDAS的输出产品不仅为NCEP的全球预报系统（GFS）提供初始条件，还为其他气象研究和应用提供了高质量的初始数据。质量控制后的数据
   2.prepbufr 数据介绍
   PreBUFR格式的数据是一种用于存储气象观测数据的临时二进制格式。它采用高效的方式存储各种类型的观测数据，如卫星、雷达和地面观测数据。这种格式在存储时会将数据压缩，以减少文件大小，同时保留数据的精度和完整性。PreBUFR文件的设计目标是为了在数据同化和预报系统中快速读取和处理数据。在处理完成后，PreBUFR数据通常会被转换为标准的BUFR格式。
   prepbufr数据使用指南（英文）

prebufr结构

我在阅读很多关于prepbufr的文件后发现其答题结构如下
prebufr数据以单个事件储存，每个有四类词条（每个不一定全部有对应的所有描述），其中hdstr是头文件相当于对于一个事情的定义（时间地点（经纬、高度）观测类型），SAID、T29可以没有，但前面个的站点名称、经度、纬度、时间间距（距离所同化时间的间距h，时间在GDAS数据名字里有）、prepbufr类型、海拔必须全面。

# BUFR 文件配置
hdstr='SID XOB YOB DHR TYP ELV SAID T29'                        #{站点名 经度 纬度 时间间距（h） prebufr类型 海拔（m） 卫星标识 输入类型}
obstr='POB QOB TOB ZOB UOB VOB PWO MXGS HOVI CAT PRSS TDO PMO'  #观测{压强 比湿 温度 高度 U V 总降水 未知 未知 数据类型 表面气压 露点温度 平均海平面压力}
qcstr='PQM QQM TQM ZQM WQM NUL PWQ PMQ'                         #质量{压强 比湿 温度 高度 风 未知 总降水 平均海平面压力}
oestr='POE QOE TOE NUL WOE NUL PWE'                             #误差{压强 比湿 温度 未知 风 未知 降水}

站点可视化

将prebufr数据主要信息导入新的text文件

程序如下

from __future__ import print_function
import ncepbufr #github里有下载教程，链接在上方
import sys

def save_print_output_to_file(file_path, *args, **kwargs):
    with open(file_path, 'a') as f:  # 使用 'a' 模式以追加模式打开文件
        sys.stdout = f  # 重定向标准输出流到文件
        print(*args, **kwargs)
        sys.stdout = sys.__stdout__  # 恢复标准输出流
# 假设 hdstr、obstr、qcstr、oestr 是你的字符串
# 假设 file_path 是你希望保存的文件路径
file_path = '/bufr_work/prebufr_data/prebufr_output.txt'#保存生成的text文件地址

# 在循环外部打开文件
with open(file_path, 'w') as f:
    pass  # 确保文件是空的，或者你可以选择清除文件内容
6

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("/prebufr_diy/2019080112.prebufr.nr")#需要解码的prepbufr数据地址（我这里是自己编译的prebufr不是GDAS的）
bufr.print_table() # print embedded table
bufr.dump_table("/obs/gdas.t12z.prepbufr.table") # dump table to file(解码需要的辅助文件)
while bufr.advance() == 0: # loop over messages.
    print(bufr.msg_counter, bufr.msg_type, bufr.msg_date, bufr.receipt_time)
    #bufr.read_subset(obstr) # should raise subset not loaded error
    while bufr.load_subset() == 0: # loop over subsets in message.
        hdr = bufr.read_subset(hdstr).squeeze()
        station_id = hdr[0].tostring()
        obs = bufr.read_subset(obstr)
        nlevs = obs.shape[-1]
        oer = bufr.read_subset(oestr)
        qcf = bufr.read_subset(qcstr)
        #print('station_id, lon, lat, time, station_type, levels =',\
        #station_id,hdr[1],hdr[2],hdr[3],int(hdr[4]),nlevs)
        # 假设 file_path 是你希望保存的文件路径
        # 将 print 输出保存到文件中
        save_print_output_to_file(file_path,'station_id,lon,lat, time,prebufr_type,input_type =',station_id,hdr[1],hdr[2],hdr[3],hdr[4],hdr[-1])

    # stop after first 2 messages.
    #if bufr.msg_counter == 2: break
bufr.close()
print('大王，数据处理完毕*—*')

值得注意的是其中需要的prepbufr.table文件，它相当于你翻译英文时用到的字典，将对应的二进制数据翻译成text在github的ncepbufr库里prepbufr.table可以下载到，prebufr官网也能下载到。

站点可视化（经纬度分布，类型区分）

在TYP中，以不同的数字对数据类型进行分类一般大于200是风资料，小于的是质量点信息即温压湿等标量资料prepbufr type table
代码如下

import codecs
import numpy as np
import matplotlib
matplotlib.use