上一篇文章《xarray指南：插值》中介绍使用 xarray 实现插值，本文介绍如何将该插值应用到 GRIB 2 文件中的要素场。

本文使用 nwpc-oper/nwpc-data 库从本地 GRIB 2 文件中加载要素场。nwpc-data 库介绍请浏览《nwpc-data 库简介》

准备

加载需要的库import numpy as npimport pandas as pdimport matplotlib.pyplot as pltimport xarray as xrfrom nwpc_data.data_finder import find_local_filefrom nwpc_data.grib.eccodes import load_field_from_file

站点插值

获取 850hPa 温度场data_file = find_local_file( "grapes_gfs_gmf/grib2/orig", start_time=pd.to_datetime("2020-04-26"), forecast_time="24h",)data_filePosixPath("/g1/COMMONDATA/OPER/NWPC/GRAPES_GFS_GMF/Prod-grib/2020042521/ORIG/gmf.gra.2020042600024.grb2")t850 = load_field_from_file( data_file, parameter="t", level_type="pl", level=500,)t850array([[262.26335937, 262.28335938, 262.29335938, ..., 262.26335937, 262.27335937, 262.30335938], ..., [243.63335937, 243.48335937, 243.64335938, ..., 243.56335938, 243.53335938, 243.73335937]])Coordinates: time datetime64[ns] 2020-04-26 step timedelta64[ns] 1 days valid_time datetime64[ns] 2020-04-27 pl float64 850.0 * latitude (latitude) float64 89.88 89.62 89.38 ... -89.38 -89.62 -89.88 * longitude (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8Attributes: GRIB_edition: 2 ...skip... long_name: Temperature units: K

获取北京站温度，[39.8, 116.4667]

索引，最近邻查找，返回距离该点最新的点。注意观察返回结果的坐标值。t850.sel( latitude=39.8, longitude=116.4667, method="nearest")array(279.25335938)Coordinates: time datetime64[ns] 2020-04-26 step timedelta64[ns] 1 days valid_time datetime64[ns] 2020-04-27 pl float64 850.0 latitude float64 39.88 longitude float64 116.5Attributes: GRIB_edition: 2 ...skip... long_name: Temperature units: K

插值t850.interp( latitude=39.8, longitude=116.4667,)array(279.34686817)Coordinates: time datetime64[ns] 2020-04-26 step timedelta64[ns] 1 days valid_time datetime64[ns] 2020-04-27 pl float64 850.0 latitude float64 39.8 longitude float64 116.5Attributes: GRIB_edition: 2 ...skip long_name: Temperature units: K

实际上，两种方法返回的数据不一样。(t850.sel( latitude=39.8, longitude=116.4667, method="nearest") - t850.interp( latitude=39.8, longitude=116.4667,)).item()-0.09350879999999506

同时插值多个站点

获取多个站点的温度：北京：[39.8, 116.4667]

上海：[31.3908, 121.4447]

广州：[23.21, 113.4822]

构建站点坐标对象，并创建新的维度 locationlat = xr.DataArray( [39.8, 31.3908, 23.21], dims="location", coords={ "location": ["bj", "sh", "gz"], })lon = xr.DataArray( [116.4667, 121.4447, 113.4822], dims="location", coords={ "location": ["bj", "sh", "gz"], })

使用新维度插值t850.interp( latitude=lat, longitude=lon)array([279.34686817, 282.03347057, 285.58804322])Coordinates: time datetime64[ns] 2020-04-26 step timedelta64[ns] 1 days valid_time datetime64[ns] 2020-04-27 pl float64 850.0 latitude (location) float64 39.8 31.39 23.21 longitude (location) float64 116.5 121.4 113.5 * location (location)

降尺度

比较使用 xarray 降尺度和目前业务系统中使用 gribpost.exe 生成的降尺度产品的差异。

获取降尺度产品中 850hPa 的温度场down_data_file = ( "/g2/nwp_pd/NWP_PST_DATA/GMF_GRAPES_GFS_POST" "/togrib2/output_togrib2_downscaling/2020042600/" "grapes_downscaling.2020042600024.grb2")down_t850 = load_field_from_file( down_data_file, parameter="t", level_type="pl", level=850,)down_t850array([[263.58795313, 263.59595312, 263.60395313, ..., 268.64095313, 268.70695312, 268.77295313], ..., [291.28795313, 291.38595313, 291.48395313, ..., 290.36095312, 290.11695313, 289.87295312]])Coordinates: time datetime64[ns] 2020-04-26 step timedelta64[ns] 1 days valid_time datetime64[ns] 2020-04-27 pl float64 850.0 * latitude (latitude) float64 70.0 69.9 69.8 69.7 ... -14.8 -14.9 -15.0 * longitude (longitude) float64 45.0 45.1 45.2 45.3 ... 144.8 144.9 145.0Attributes: GRIB_edition: 2 ...skip... long_name: Temperature units: K

使用 interp_like 按 down_t850 的网格对 t850 进行插值interpolated = t850.interp_like(down_t850)interpolatedarray([[263.58835937, 263.59635937, 263.60435938, ..., 268.64135938, 268.70735938, 268.77335938], ..., [291.28835938, 291.38635938, 291.48435938, ..., 290.36135937, 290.11735937, 289.87335938]])Coordinates: time datetime64[ns] 2020-04-26 step timedelta64[ns] 1 days valid_time datetime64[ns] 2020-04-27 pl float64 850.0 * latitude (latitude) float64 70.0 69.9 69.8 69.7 ... -14.8 -14.9 -15.0 * longitude (longitude) float64 45.0 45.1 45.2 45.3 ... 144.8 144.9 145.0Attributes: GRIB_edition: 2 ...skip... long_name: Temperature units: K

计算两者的差值d = interpolated - down_t850darray([[ 4.06250000e-04, 4.06250000e-04, 4.06250000e-04, ..., 4.06250000e-04, 4.06250000e-04, 4.06250000e-04], [ 4.06250000e-04, -1.93750000e-04, 2.06250000e-04, ..., 6.25000001e-06, 2.06250000e-04, 4.06250000e-04], [ 4.06250000e-04, -1.93750000e-04, 2.06250000e-04, ..., -1.93750000e-04, -3.93750000e-04, 4.06250000e-04], ..., [ 4.06250000e-04, 4.06250000e-04, 4.06250000e-04, ..., 6.24999996e-06, 2.06250000e-04, 4.06250000e-04], [ 4.06250000e-04, 2.06250000e-04, 6.25000001e-06, ..., -3.93750000e-04, 6.25000001e-06, 4.06250000e-04], [ 4.06250000e-04, 4.06250000e-04, 4.06250000e-04, ..., 4.06250000e-04, 4.06250000e-04, 4.06250000e-04]])Coordinates: time datetime64[ns] 2020-04-26 step timedelta64[ns] 1 days valid_time datetime64[ns] 2020-04-27 pl float64 850.0 * latitude (latitude) float64 70.0 69.9 69.8 69.7 ... -14.8 -14.9 -15.0 * longitude (longitude) float64 45.0 45.1 45.2 45.3 ... 144.8 144.9 145.0

计算所有差值绝对值之和np.abs(d).sum()array(256.68321875)Coordinates: time datetime64[ns] 2020-04-26 step timedelta64[ns] 1 days valid_time datetime64[ns] 2020-04-27 pl float64 850.0

查看差值的最大值和最小值d.min().item()-0.0003937500001711669d.max().item()0.00040625000019645086

可以看到，两者的差值在 0.001 级别，考虑到 JPEG 2000 压缩的精度，误差在可接受范围。

参考

https://github.com/nwpc-oper/nwpc-data

题图来自 Pixabay。