GRIB笔记：使用 xarray 插值
目录
上一篇文章《xarray指南：插值》中介绍使用 xarray 实现插值，本文介绍如何将该插值应用到 GRIB 2 文件中的要素场。
本文使用 nwpc-oper/nwpc-data 库从本地 GRIB 2 文件中加载要素场。 nwpc-data 库介绍请浏览《nwpc-data 库简介》
准备
加载需要的库
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr
from nwpc_data.data_finder import find_local_file
from 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,
)
t850<xarray.DataArray 't' (latitude: 720, longitude: 1440)>
array([[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.8
Attributes:
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.5
Attributes:
GRIB_edition: 2
...skip...
long_name: Temperature
units: K
插值
t850.interp(
latitude=39.8,
longitude=116.4667,
)<xarray.DataArray 't' ()>
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.5
Attributes:
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]
构建站点坐标对象，并创建新的维度 location
lat = 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
)<xarray.DataArray 't' (location: 3)>
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) <U2 'bj' 'sh' 'gz'
Attributes:
GRIB_edition: 2
...skip...
long_name: Temperature
units: K
降尺度
比较使用 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_t850<xarray.DataArray 't' (latitude: 851, longitude: 1001)>
array([[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.0
Attributes:
GRIB_edition: 2
...skip...
long_name: Temperature
units: K
使用 interp_like 按 down_t850 的网格对 t850 进行插值
interpolated = t850.interp_like(down_t850)
interpolated<xarray.DataArray 't' (latitude: 851, longitude: 1001)>
array([[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.0
Attributes:
GRIB_edition: 2
...skip...
long_name: Temperature
units: K
计算两者的差值
d = interpolated - down_t850
d<xarray.DataArray 't' (latitude: 851, longitude: 1001)>
array([[ 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()<xarray.DataArray 't' ()>
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.0003937500001711669
d.max().item()0.00040625000019645086
可以看到，两者的差值在 0.001 级别，考虑到 JPEG 2000 压缩的精度，误差在可接受范围。
