Meteva笔记:加载GRIB 2要素场
Meteva 是由 nmc 开源的全流程检验程序库,提供了常用的各种气象预报检验评估的算法函数,气象检验分析的图片和表格型产品的制作函数,以及检验评估系统示例。
本文介绍如何通过 nwpc-data 库将本地 GRIB 2 文件接入到 Meteva 工具中。
准备
本文代码均在 CMA-PI 高性能计算机上运行。
载入需要的库
import pandas as pd
import numpy as np
import xarray as xr
import meteva.base as meb
import meteva.method as mem
from nwpc_data.grib.eccodes import load_field_from_file
from nwpc_data.data_finder import find_local_file
meb
提供 IO 和绘图相关函数,mem
提供计算检验指标的函数。
从 GDS 加载数据
在加载本地数据文件前,首先使用 Meteva 内置的函数从 GDS 服务中获取要素场,用于后续对比验证。
GDS 服务的相关配置方法请访问同样由 nmcdev 开源的 nmcdev/nmc_met_io 项目。
Meteva 支持 nmc_met_io 项目的配置文件。本文使用已配置好的文件:
gds_config_file = "/g1/u/wangdp/.config/.nmcdev/config.ini"
使用 meb.io.read_gds_ip_port
从配置文件中读取 IP 地址和端口号
ip, port = meb.io.read_gds_ip_port(gds_config_file)
ip, port
('10.32.8.164', 8080)
格点数据
加载 GRAPES GFS 的格点数据。
构建数据路径
start_time = pd.Timestamp("2020-09-16 08:00:00")
start_time
Timestamp('2020-09-16 08:00:00')
path = meb.get_path("GRAPES_GFS/TMP/850/YYMMDDHH.TTT", start_time, 24)
path
'GRAPES_GFS/TMP/850/20091608.024'
使用 meb.read_griddata_from_gds()
函数从 GDS 服务中获取数据,返回 xr.DataArray
对象
t850_grid = meb.read_griddata_from_gds(ip, port, path)
t850_grid
使用 meb.set_griddata_coords()
函数修改坐标
meb.set_griddata_coords(
t850_grid,
name="t",
member_list = ["GRAPES_GFS"],
gtime = [start_time],
dtime_list = [24]
)
t850_grid
使用 Meteva 内置的绘图函数 meb.tool.plot_tools.contourf_2d_grid()
绘制等值线图
meb.tool.plot_tools.contourf_2d_grid(
t850_grid,
)
读取本地 GRIB 2 数据
载入
文件路径
file_path = find_local_file(
"grapes_gfs_gmf/grib2/orig",
start_time="2020091600",
forecast_time="24h",
)
file_path
PosixPath('/g1/COMMONDATA/OPER/NWPC/GRAPES_GFS_GMF/Prod-grib/2020091600/ORIG/gmf.gra.2020091600024.grb2')
读取 850hPa 温度场
单个要素场默认只有两个维度:latitude 和 longitude
field = load_field_from_file(
file_path,
parameter="t",
level_type="pl",
level=850
)
field
转换
将 time,step 和 pl 都扩展为维度,并将单位转为摄氏度
field = field.expand_dims(["time", "step", "pl"])
field = field - 273.15
field
使用 meb.xarray_to_griddata()
函数将要素场对象转为 meb.grid_data()
函数生成的 xr.DataArray
对象
可以看到,对于单个要素场,该函数自动生成了 memeber
维度,坐标值为 0
grid_data = meb.xarray_to_griddata(
field,
level_dim="pl",
time_dim="time",
dtime_dim="step",
lat_dim="latitude",
lon_dim="longitude",
)
grid_data
裁剪成与 t850_grid
相同的网格
cropped_grid_data = grid_data.sel(
lat=slice(-9.875, 79.88),
lon=slice(0.0, 180.0)
)
cropped_grid_data
验证
对比本地提取的要素场和从 GDS 中获取的要素场是否相同。
从下图中可以看到,member
的名称不同,同时 time
也不同。GDS 使用的是北京时间,而本地文件使用世界时。
修改坐标值,与 t850_grid
保持一致:
- 添加 member 名称
- 将时间从世界时改为北京时
meb.set_griddata_coords(
cropped_grid_data,
member_list=["GRAPES_GFS"],
gtime=cropped_grid_data.time + pd.Timedelta(hours=8)
)
cropped_grid_data
求两个场的偏差
diff_t850 = cropped_grid_data - t850_grid
diff_t850
求偏差场中最大偏差
abs(diff_t850).max()
0.01001473
差值可能是因为压缩精度的问题,在可以接受的范围内。说明本地读取的 GRIB 2 文件可以代替 GDS 中的数据。
计算
计算 024 时效与该时刻分析场的均方根误差
载入数据
整合函数,实现如下功能:
- 使用 nwpc-data 从 GRIB 2 文件中加载要素场
- 将返回的要素场转换为
xr.DataArray
对象
def load_grid(
file_path,
parameter,
level_type=None,
level=None,
) -> xr.DataArray:
field = load_field_from_file(
file_path,
parameter=parameter,
level_type=level_type,
level=level,
)
field = field.expand_dims(["time", "step", "pl"])
grid = meb.xarray_to_griddata(
field,
level_dim="pl",
time_dim="time",
dtime_dim="step",
lat_dim="latitude",
lon_dim="longitude",
)
return grid
分析场
file_path = find_local_file(
"grapes_gfs_gmf/grib2/orig",
start_time="2020091700",
forecast_time="0h",
)
file_path
PosixPath('/g1/COMMONDATA/OPER/NWPC/GRAPES_GFS_GMF/Prod-grib/2020091700/ORIG/gmf.gra.2020091700000.grb2')
anal_grid = load_grid(
file_path,
parameter="t",
level_type="pl",
level=850
)
anal_grid
预报场
file_path = find_local_file(
"grapes_gfs_gmf/grib2/orig",
start_time="2020091600",
forecast_time="24h",
)
file_path
PosixPath('/g1/COMMONDATA/OPER/NWPC/GRAPES_GFS_GMF/Prod-grib/2020091600/ORIG/gmf.gra.2020091600024.grb2')
fcst_grid = load_grid(
file_path,
parameter="t",
level_type="pl",
level=850
)
fcst_grid
计算指标
计算均方根误差 RMSE
使用 squeeze
方法删掉长度为 1 的维度,将数据变为二维矩阵
mem.rmse(
anal_grid.squeeze(),
fcst_grid.squeeze(),
)
计算多个预报数据的指标
加载另一个数据:48 小时预报
file_path = find_local_file(
"grapes_gfs_gmf/grib2/orig",
start_time="2020091500",
forecast_time="48h",
)
fcst_grid_48 = load_grid(
file_path,
parameter="t",
level_type="pl",
level=850
)
两个数据仅预报时效不同,沿 dtime
维度合并两个要素场
fcst_grids = xr.concat(
[
grid.squeeze().reset_coords(["member", "level", "time"], drop=True)
for grid in [fcst_grid, fcst_grid_48]
],
dim="dtime"
)
fcst_grids
计算均方根误差
mem.rmse(
anal_grid.squeeze(),
fcst_grids.values,
)
array([1.13974815, 1.85073222])
参考
Meteva 项目:
https://github.com/nmcdev/meteva
nmc_met_io 项目:
https://github.com/nmcdev/nmc_met_io
nwpc-data 项目:
https://github.com/nwpc-oper/nwpc-data
- 关于字符串匹配查找的总结(43天)
- 一条sql语句导致的数据库宕机问题及分析(42天)
- 外部表的导入导出问题 (41天)
- 当我们和计算机交互时,它看到的是什么?
- 一条sql语句“导致”的数据库宕机问题及分析 (38天)
- rman数据备份恢复学习笔记(49天)
- 虚拟专用数据库VPD应用 (48天)
- 关于创建视图的问题(48天)
- 性能调优之redo切换频率(47天)
- 关于oracle中session跟踪的总结(56天)
- oracle中关于小数中0的格式化(55天)
- 关于trigger过滤最大值的问题(54天)
- oracle共享服务器配置汇总(53天)
- 关于drop user的cascade选项解惑(52天)
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法