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

最新推荐文章于 2024-10-02 10:53:34 发布
最新推荐文章于 2024-10-02 10:53:34 发布
阅读量943 收藏 27
点赞数 10
文章标签： python 数据库 开发语言
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/shuaibiXiangzai/article/details/141714827
版权

当我们使用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('Agg')
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.lines import Line2D

# 文件路径
filepath = "/bufr_work/prebufr_data/prebufr_output.txt"
save_path = "/bufr_work/2019080112global_prebufr_sit.png"#图片保存地址

# 读取数据
lon, lat, prebufr_type, input_type = [], [], [], []
with codecs.open(filepath, mode='r', encoding='utf-8') as f:
    for line in f:
        a = line.split()
        if len(a) >= 6:
            try:
                b = np.float32(a[5])  # 使用 float32 提高精度
                c = np.float32(a[6])
                d = np.float32(a[8])  # 转换为浮点数
                lon.append(b)
                lat.append(c)
                prebufr_type.append(d)
                #input_type.append(e)
            except ValueError as e:
                print(f"数据转换错误: {e}，行内容: {line}")
        else:
            print(f"行数据不完整，已跳过: {line}")

# 计算 prebufr_type 的唯一值数量
unique_prebufr_types = set(prebufr_type)
num_unique_prebufr_types = len(unique_prebufr_types)

# 打印唯一值和总数量
print(f"Unique prebufr types: {unique_prebufr_types}")
print(f"Number of unique prebufr types: {num_unique_prebufr_types}")
print(f"Total number of prebufr types: {len(prebufr_type)}")

# 定义颜色和标签字典
colors = {
    'SATWND & ASCATW': 'blue',
    'ADPSFC': 'red',
    'SFCSHP': 'green',
    'ADPUPA': 'orange',
    'AIRCFT': 'purple',
    'VADWND & RASSDA': 'brown'
}

categories = {
    'SATWND & ASCATW': [242, 243, 250, 252, 253, 254, 290],
    'ADPSFC': [181, 187, 281, 287, 284, 183],
    'SFCSHP': [180, 280, 282],
    'ADPUPA': [120, 220, 221],
    'AIRCFT': [130, 230],
    'VADWND & RASSDA': [126, 224]
}

# 计算每个类别的点数量
category_counts = {}
for category, prebufr_types in categories.items():
    count = sum(1 for t in prebufr_type if t in prebufr_types)
    category_counts[category] = count

# 根据点数量对类别进行排序
sorted_categories = sorted(category_counts.items(), key=lambda item: item[1], reverse=True)

fig = plt.figure(figsize=(16, 12), dpi=800)
ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

# 添加地图特征
ax.add_feature(cfeature.LAND, facecolor='lightgray')
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=':')

# 绘制实际散点图
for category, count in sorted_categories:
    prebufr_types = categories[category]
    lons = [lon[i] for i in range(len(lon)) if prebufr_type[i] in prebufr_types]
    lats = [lat[i] for i in range(len(lat)) if prebufr_type[i] in prebufr_types]
    size = 2 if category in ['SATWND & ASCATW', 'ADPSFC'] else 10  # 根据类别设置实际大小
    color = colors[category]
    facecolor = 'none' if category == 'SATWND & ASCATW' else color  # 空心圆点
    edgecolor = 'blue' if category == 'SATWND & ASCATW' else color  # 设置边缘颜色
    ax.scatter(lons, lats, marker='o', facecolors=facecolor, edgecolors=edgecolor, linewidths=0.5, s=size, transform=ccrs.PlateCarree())

# 创建图例句柄和标签
legend_handles = []
legend_labels = []
for category, color in colors.items():
    if category == 'SATWND & ASCATW':
        handle = Line2D([0], [0], marker='o', color='w', markerfacecolor='none', markeredgecolor='blue', markeredgewidth=1, markersize=12, linestyle='None')
    else:
        handle = Line2D([0], [0], marker='o', color='w', markerfacecolor=color, markersize=12, linestyle='None')
    legend_handles.append(handle)
    label = f"{category} ({category_counts[category]})"
    legend_labels.append(label)

# 创建图例
ax.legend(handles=legend_handles, labels=legend_labels, title='PreBUFR Type Categories', loc='lower left')

# 设置地图范围
ax.set_extent([0, 359, -89, 89], crs=ccrs.PlateCarree())

# 添加网格线
ax.gridlines(draw_labels=True)

# 添加标题和标签
plt.title('Observation Sites(2024_06_28_12)')
plt.xlabel('Longitude')
plt.ylabel('Latitude')

# 调整布局
plt.tight_layout()

# 保存图像
plt.savefig(save_path, format='png', dpi=800)

print('图已经画好了')

需要注意的地方是，我先打印出我所使用的prepbufr数据所有类型，才根据打印出的类型来设立的分类category，实际的分类当然不止这些，想作全面的可以参考prebufr type的详细分类prepbufr type table，然后由于卫星的风数据站点太多，我将对应的格点做了相应的调整来美化分布图
观测站点全球分布图
你也可只画你同化范围的站点分布
nearby
不颜色的类型分类在官网中都有解释有地面站、卫星风、雷达、飞机等。
附带各种风类分开画的程序

import codecs
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import os

# 文件路径
filepath = "/obs/prebufr_output.txt"
save_folder = "/bufr_work/prebufr site sort"

# 确保保存文件夹存在
os.makedirs(save_folder, exist_ok=True)

# 读取数据
lon, lat, prebufr_type = [], [], []
with codecs.open(filepath, mode='r', encoding='utf-8') as f:
    for line in f:
        a = line.split()
        if len(a) >= 6:
            try:
                b = np.float32(a[2])
                c = np.float32(a[3])
                d = np.float32(a[4])
                lon.append(b)
                lat.append(c)
                prebufr_type.append(d)
            except ValueError as e:
                print(f"数据转换错误: {e}，行内容: {line}")
        else:
            print(f"行数据不完整，已跳过: {line}")

# 定义分类及其颜色
categories = {
    'SATWND & ASCATW': [242, 243, 250, 252, 253, 254, 290],
    'ADPSFC': [181, 187, 281, 287, 284, 183],
    'SFCSHP': [180, 280, 282],
    'ADPUPA': [120, 220, 221],
    'AIRCFT': [130, 230],
    'VADWND & RASSDA': [126, 224]
}

colors = {
    'SATWND & ASCATW': 'blue',
    'ADPSFC': 'red',
    'SFCSHP': 'green',
    'ADPUPA': 'orange',
    'AIRCFT': 'purple',
    'VADWND & RASSDA': 'brown'
}

# 为每个分类绘制图像
for category, prebufr_types in categories.items():
    lons = [lon[i] for i in range(len(lon)) if prebufr_type[i] in prebufr_types]
    lats = [lat[i] for i in range(len(lat)) if prebufr_type[i] in prebufr_types]

    fig = plt.figure(figsize=(10, 6), dpi=800)
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())

    # 添加地图特征
    ax.add_feature(cfeature.LAND, facecolor='lightgray')
    ax.add_feature(cfeature.COASTLINE)
    ax.add_feature(cfeature.BORDERS, linestyle=':')

    # 设置散点大小
    size = 1 if category in ['SATWND & ASCATW', 'ADPSFC'] else 3

    # 绘制散点，点大小设置为根据分类调整
    ax.scatter(lons, lats, color=colors[category], marker='o', s=size, transform=ccrs.PlateCarree())

    # 计算当前类别的点数量并构建图例标签
    count = len(lons)
    legend_label = f'{category} ({count})'

    # 创建图例句柄
    ax.legend([plt.Line2D([0], [0], marker='o', color='w', markerfacecolor=colors[category], markersize=10)], 
              [legend_label], loc='lower left')

    # 设置地图范围
    ax.set_extent([0, 359, -89, 89], crs=ccrs.PlateCarree())

    # 添加网格线
    ax.gridlines(draw_labels=True)

    # 添加标题和标签
    plt.title(f'Observation Sites - {category}')
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')

    # 保存图像
    save_path = os.path.join(save_folder, f'{category.replace(" ", "_")}.png')
    plt.savefig(save_path, format='png', dpi=800)
    plt.close()  # 关闭图形，以释放内存

    print(f'图像已保存: {save_path}')

放两张效果图
地面站
船只走航

阿翔_ocean
关注
常见气象数据获取方式及批量下载代码汇总
10-20 2万+
26个来源的气象数据获取代码
MATLAB读取gdas1气象数据,HYSPLIT模式简介及单机版使用图文教程
weixin_32258183的博客
03-19 1754
模式介绍HYSPLIT - Hybrid Single Particle Lagrangian IntegratedTrajectory ModelHYSPLIT - 4 模型是由美国国家海洋和大气管理局 (NOAA)的空气资源实验室和澳大利亚气象局在过去 20年间联合研发的一种用于计算和分析大气污染物输送、扩散轨迹的专业模型。该模型具有处理多种气象要素输入场、多种物理过程和不同类型污染物排放源功...
气象数据bufr格式解析
qq_33381052的博客
06-24 3841
气象数据、BUFR数据解析
利用seawind数据对WRF模式和GSI同化的风结果进行评估（python）
shuaibiXiangzai的博客
08-31 864
评估是模式预报和资料同化的重要一步，这里主要介绍一下用一般的bias和RMSE来检验模拟的可靠性，虽然结果不是很乐观，但不影响评估方法。接下来简单介绍一下实验所用的数据和方法GFS数据是由全球预报系统（GFS）模式gfs官网生成的全球数值天气预报数据。该数据覆盖全球，从地表到高层大气，提供多种气象变量的信息，如温度、风速、湿度、降水量、气压等。GFS数据以较高的空间分辨率（通常为0.25度或0.5度）和时间分辨率（1小时或3小时）提供，适用于天气预报、气候分析、环境监测等领域。
管道操作——为GIS准备GDAS气象数据
丁劲犇技术发布空间
06-11 6530
很多出国深造的同学，都对国外高校中的计算机教学、使用记忆犹新。国内一般院校的老师很多都是从微软的DOS起步开始捣鼓微型计算机的，基本上对unix系统用的不多。对命令行操作，也停留在dos命令的概念上。最近，一位同学毕业设计遇到了读取天气预报数据的问题，在linux下摸索了很久，还是蛮有意思的。征得他的同意，把他的问题作为unix管道操作一例分享一下。 1 问题背景 该学生毕业设计是做一个天......
关于运行GDAS源码的坑及解决方案
huqiaolong的博客
05-08 367
前提：请明确源码一般是没什么问题的，大多数运行不成功都是配置问题 坑：运行scripts-cnn/train-cifar.sh内存溢出的问题 原因：这是pytorch源码，指定GPU不能在命令行中指定，需要在入口py文件中指定 措施：在exps-cnn/train_base.py中加两行 import os os.environ["CUDA_VISIBLE_DEVICES"] = "2"...
GDAS下载方式
Fairy_Nan的博客
12-20 5256
https://zhuanlan.zhihu.com/p/136989134
python安装用于解析GDAS资料同化数据ARLreader环境
BelievIamthebest的博客
06-19 220
【代码】python安装用于解析GDAS资料同化数据ARLreader环境。
python——GDAS1数据处理并转为netCDF格式
ljt_217的博客
02-13 1581
GDAS1数据处理，并转为nc格式文件
Python库 | gdas-0.2.4.tar.gz
04-09
资源分类：Python库 所属语言：Python 资源全名：gdas-0.2.4.tar.gz 资源来源：官方 安装方法：https://lanzao.blog.csdn.net/article/details/101784059
NCEPLIBS-bufr-develop工具；气象数据解析插件，bufr功能
06-24
这种格式可以有效地压缩数据，减少存储空间，并且能够容纳各种类型和来源的气象信息，包括气象站观测、雷达回波、卫星图像等。 NCEPLIBS-bufr-develop工具是基于Python编写的，这使得它具有良好的跨平台性和易用性...
【训练营数据】气象业务数据处理GP2.rar
08-23
预处理的BUFR文件（prepbufr.gdas）是气象数据的一种标准编码格式，由世界气象组织（WMO）制定。BUFR（Binary Universal Form for the Representation of Meteorological Data）是一种二进制编码，能够高效地存储...
hysplit4气象数据数据下载地址
05-05
hysplit是NOAA研发的一款进行气流反演的一款软件，通过气象数据进行分析大气输送通道。数据下载地址
（笔记）第三期书生·浦语大模型实战营（十一卷王场）--书生入门岛通关第2关Python 基础知识
最新发布
haidizym的博客
10-02 477
学员闯关手册：https://aicarrier.feishu.cn/wiki/ZcgkwqteZi9s4ZkYr0Gcayg1n1g?课程视频：https://www.bilibili.com/video/BV1mS421X7h4/课程文档：https://github.com/InternLM/Tutorial/tree/camp3/docs/L0/Python。
遥感影像-实例分割数据集：iSAID 从切图到YOLO格式数据集制作详细介绍
GIS潮流
09-30 1236
开源数据集isaid标注包含实例分割，但是原始影像太大，很吃显存，一般显卡无法用原始影像直接训练，所以需要对影像进行裁剪，并生成对应的标签，因为想用yolo系列跑模型，所以将标签需要转为txt格式。
神经网络在多分类问题中的应用
小高~的博客
09-29 895
神经网络在分类任务中的应用越来越广泛，尤其是在图像识别、自然语言处理等领域。本文将介绍如何使用PyTorch构建一个简单的神经网络来处理多分类问题。我们将通过一个实战案例，展示数据构造、模型训练、模型保存与加载、以及评估结果的整个过程。
酷炫音乐盒: 使用Python和Tkinter配合Pygame打造自己的音乐播放器
宇宙第一小趴菜的博客
09-29 1238
Tkinter是一个功能强大且简单易用的库，适合构建轻量级的GUI应用程序。它在不同的操作系统上具有相同的外观和行为，因此可以实现跨平台的应用程序。当然，Tkinter还有更多的功能和细节，通过阅读官方文档和教程，你可以更深入地了解和学习如何运用Tkinter来创建各种GUI应用程序。
gdas气象数据在哪下载
12-17
gdas气象数据可以在美国国家环境信息中心（National Centers for Environmental Information, NCEI）的官方网站上进行下载。在该网站上，用户可以找到GDAS气象数据的各种格式和版本，包括GRIB和NetCDF格式等。用户可以根据自己的需求选择合适的数据格式进行下载，而且下载过程非常简单和便捷。 另外，用户们也可以通过一些气象数据共享平台或者气象数据交换组织来获取GDAS气象数据，这些组织和平台通常会提供GDAS气象数据的下载链接或者共享渠道。通过这些途径，用户可以充分利用GDAS气象数据进行气象预测、研究和分析等工作。 除此之外，一些科研机构或者大学也可能提供GDAS气象数据的下载服务，用户可以通过这些机构的官方网站或者相关科研项目的网站来获取GDAS气象数据。总之，用户可以通过多种途径来下载GDAS气象数据，可以选择最适合自己需求的渠道来获取所需的气象数据。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

阿翔_ocean

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值