pygrib#
本节介绍如何使用 pygrib 加载 GRIB2 要素场,并使用 cedarkit-maps 绘图。
安装#
使用 conda 安装 pygrib
conda install -c conda-forge pygrib
准备#
导入需要的包
import xarray as xr
import pandas as pd
import pygrib
设置绘图的数据参数,使用 CMA-MESO 2024 年 4 月 1 日 00 时次 024 时效数据。
system_name = "CMA-MESO"
data_type = "cma_meso_3km/grib2/orig"
start_time = pd.to_datetime("2024-04-01 00:00:00")
forecast_time = pd.to_timedelta("24h")
加载数据#
设置 GRIB2 数据文件路径
注:可以使用 reki 库查找本地文件路径
file_path = '/g3/COMMONDATA/OPER/CEMC/MESO_3KM/Prod-grib/2024040100/ORIG/rmf.hgra.2024040100024.grb2'
file_path
'/g3/COMMONDATA/OPER/CEMC/MESO_3KM/Prod-grib/2024040100/ORIG/rmf.hgra.2024040100024.grb2'
注:可以使用 reki 库查找本地文件路径
Show code cell source
from reki.data_finder import find_local_file
file_path_use_reki = find_local_file(
data_type,
start_time=start_time,
forecast_time=forecast_time,
)
file_path_use_reki
Show code cell output
PosixPath('/g3/COMMONDATA/OPER/CEMC/MESO_3KM/Prod-grib/2024040100/ORIG/rmf.hgra.2024040100024.grb2')
使用 pygrib 加载 2 米温度要素场。
grbs = pygrib.open(file_path)
selected_grbs = grbs.select(shortName="2t")
selected_grb = selected_grbs[0]
selected_grb
10:2 metre temperature:K (instant):regular_ll:heightAboveGround:level 2 m:fcst time 24 mins:from 202404010000
获取经纬度坐标
lats, lons = selected_grb.latlons()
lats[:, 0], lons[0, :]
(array([60.099998, 60.069998, 60.039998, ..., 10.059998, 10.029998,
10. ]),
array([ 70. , 70.03, 70.06, ..., 144.94, 144.97, 145. ]))
获取要素场二维数组
values = selected_grb.values
values.shape
(1671, 2501)
将要素场包装为 xarray.DataArray
对象
t_2m_field = xr.DataArray(
values,
coords={
"latitude": lats[:, 0],
"longitude": lons[0, :],
},
dims=["latitude", "longitude"]
) - 273.15
t_2m_field
<xarray.DataArray (latitude: 1671, longitude: 2501)> Size: 33MB array([[-0.71931641, -0.75931641, -0.80931641, ..., -4.16931641, -2.87931641, -2.95931641], [-3.65931641, -3.66931641, -2.27931641, ..., -6.91931641, -6.31931641, -3.72931641], [-3.54931641, -3.58931641, -0.85931641, ..., -7.60931641, -7.47931641, -3.43931641], ..., [26.60068359, 26.60068359, 26.59068359, ..., 25.80068359, 25.79068359, 25.97068359], [26.59068359, 26.60068359, 26.59068359, ..., 25.78068359, 25.78068359, 26.00068359], [26.59068359, 26.59068359, 26.60068359, ..., 25.78068359, 25.79068359, 25.97068359]]) Coordinates: * latitude (latitude) float64 13kB 60.1 60.07 60.04 ... 10.06 10.03 10.0 * longitude (longitude) float64 20kB 70.0 70.03 70.06 ... 144.9 145.0 145.0
绘图#
使用 cedarkit-maps 绘制 2 米温度填充图
Show code cell source
from cedarkit.maps.style import ContourStyle
from cedarkit.maps.chart import Panel
from cedarkit.maps.domains import EastAsiaMapDomain
from cedarkit.maps.colormap import get_ncl_colormap
t_2m_level = [-24, -20, -16, -12, -8, -4, 0, 4, 8, 12, 16, 20, 24, 28, 32]
color_index = [2, 12, 22, 32, 42, 52, 62, 72, 82, 92, 102, 112, 122, 132, 142, 152]
t_2m_color_map = get_ncl_colormap("BlAqGrYeOrReVi200", index=color_index)
t_2m_style = ContourStyle(
colors=t_2m_color_map,
levels=t_2m_level,
fill=True,
)
domain = EastAsiaMapDomain()
panel = Panel(domain=domain)
panel.plot(t_2m_field, style=t_2m_style)
domain.set_title(
panel=panel,
graph_name="2m Temperature (C)",
system_name=system_name,
start_time=start_time,
forecast_time=forecast_time,
)
domain.add_colorbar(panel=panel, style=t_2m_style)
panel.show()