GRIB笔记：使用xarray实现二维要素场插值

meteorology/data, programming/python
目录

介绍如何使用 xarray 实现二维矩阵插值。

准备

加载 Python 常用科学库

import numpy as np
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

使用 nwpc-data 库查找数据文件，并加载要素场

from nwpc_data.data_finder import find_local_file
from nwpc_data.grib.eccodes import (
    load_field_from_file,
)

数据

data_path = find_local_file(
    "grapes_gfs_gmf/grib2/orig",
    start_time="2021031100",
    forecast_time="24h"
)
data_path

PosixPath('/g1/COMMONDATA/OPER/NWPC/GRAPES_GFS_GMF/Prod-grib/2021031100/ORIG/gmf.gra.2021031100024.grb2')

使用 load_field_from_file() 加载要素场，返回 xarray.DataArray 对象

t2m = load_field_from_file(
    data_path,
    parameter="2t",
) - 273.15
t2m


t2m.values.shape

(720, 1440)

使用 xarray 内置的 plot() 函数进行绘图

def plot(field):
    fig = plt.figure(figsize=(10, 5))
    ax = plt.axes(projection=ccrs.PlateCarree())
    field.plot(
        ax=ax
    )
    ax.coastlines()
    plt.show()

plot(t2m)

插值

interp() 函数实现多维插值

1 维数据：scipy.interpolate.interp1d()

插值方法：

  • linear
  • nearest
  • zero
  • slinear
  • quadratic
  • cubic

多维数据：scipy.interpolate.interpn()

插值方法：

  • linear：线性插值
  • nearest：最近邻插值

单点插值

默认使用线性插值

t2m.interp(
    latitude=39.8062,
    longitude=116.4693
)

最近邻插值

t2m.interp(
    latitude=39.8062,
    longitude=116.4693,
    method="nearest"
)

sel() 方法也可以实现最近邻插值

t2m.sel(
    latitude=39.8062,
    longitude=116.4693,
    method="nearest"
)

多点插值

构造新的维度 station

station_lats = xr.DataArray(
    [39.8062, 41.80], 
    coords=[["Beijing", "Shenyang"]], 
    dims="station"
)
station_lons = xr.DataArray(
    [116.4693, 123.40], 
    coords=[["Beijing", "Shenyang"]], 
    dims="station"
)
station_lats, station_lons

(<xarray.DataArray (station: 2)>
 array([39.8062, 41.8   ])
 Coordinates:
   * station  (station) <U8 'Beijing' 'Shenyang',
 <xarray.DataArray (station: 2)>
 array([116.4693, 123.4   ])
 Coordinates:
   * station  (station) <U8 'Beijing' 'Shenyang')

interp() 返回数据包含 station 维度

t2m.interp(
    latitude=station_lats,
    longitude=station_lons
)

网格插值

grid_field = t2m.interp(
    latitude=np.arange(89.95, -90, -0.1),
    longitude=np.arange(0, 360, 0.1)
)
grid_field

设置 fille_value 参数填充 nan 值

grid_field = t2m.interp(
    latitude=np.arange(89.95, -90, -0.1),
    longitude=np.arange(0, 360, 0.1),
    kwargs={
        "fill_value": "extrapolate"
    }
)
grid_field

绘图

plot(grid_field)

interp_like() 函数按另一个 DataArray 对象的维度对数据进行插值

grid = xr.DataArray(
    coords=[
        ("latitude", np.arange(89.95, -90, -0.1)),
        ("longitude", np.arange(0, 360, 0.1))
    ]
)
grid


t2m.interp_like(
    grid,
    kwargs={
        "fill_value": "extrapolate"
    }
)

参考

xarray 数据插值指南

http://xarray.pydata.org/en/stable/interpolation.html

《xarray指南：插值》

GRIB 2 数据插值

使用scipy实现二维要素场插值

使用xarray实现二维要素场插值

使用pyinterp实现二维要素场插值