Scipy 内插网格数据的问题

Question

我有一个netcdf file with a spatial resolution of 0.05º and I want to regrid it to a spatial resolution of 0.01º like this other netcdf。我尝试使用 scipy.interpolate.griddata，但我并没有真正做到这一点，我认为我缺少一些东西。

original_dataset = xr.open_dataset('to_regrid.nc')
target_dataset= xr.open_dataset('SSTA_L4_MED_0_1dg_2022-01-18.nc')

根据scipy.interpolate.griddata documentation，我需要构建我的插值管道如下：

grid = griddata(points, values, (grid_x_new, grid_y_new), method='nearest')

所以在我的例子中，我假设如下：

#Saving in variables the old and new grids
grid_x_new = target_dataset['lon']
grid_y_new = target_dataset['lat']
grid_x_old = original_dataset ['lon']
grid_y_old = original_dataset ['lat']

points = (grid_x_old,grid_y_old)
values = original_dataset['analysed_sst'] #My variable in the netcdf is the sea surface temp.

现在，当我运行 griddata:

from scipy.interpolate import griddata
grid = griddata(points, values, (grid_x_new, grid_y_new),method='nearest')

我收到以下错误：

ValueError: shape mismatch: objects cannot be broadcast to a single shape

我认为它与 lat/lon 数组形状有关。我是 netcdf 领域的新手，真的不知道这里的问题是什么。任何帮助将不胜感激！

Answer 1

我建议使用 xesm 重新网格化 xarray 数据集。下面的代码将重新网格化您的数据集：

import xarray as xr
import xesmf as xe
original_dataset = xr.open_dataset('to_regrid.nc')
target_dataset= xr.open_dataset('SSTA_L4_MED_0_1dg_2022-01-18.nc')
regridder = xe.Regridder(original_dataset, target_dataset, "bilinear")
dr_out = regridder(original_dataset)

Answer 2

在您的原始代码中，grid_x_old 和 grid_y_old 中的索引应对应于数据集中的每个唯一坐标。为了让事情正常工作，像下面这样的东西会起作用：

import xarray as xr
from scipy.interpolate import griddata
original_dataset = xr.open_dataset('to_regrid.nc')
target_dataset= xr.open_dataset('SSTA_L4_MED_0_1dg_2022-01-18.nc')
#Saving in variables the old and new grids
grid_x_old = original_dataset.to_dataframe().reset_index().loc[:,["lat", "lon"]].lon
grid_y_old = original_dataset.to_dataframe().reset_index().loc[:,["lat", "lon"]].lat

grid_x_new = target_dataset.to_dataframe().reset_index().loc[:,["lat", "lon"]].lon
grid_y_new = target_dataset.to_dataframe().reset_index().loc[:,["lat", "lon"]].lat
values = original_dataset.to_dataframe().reset_index().loc[:,["lat", "lon", "analysed_sst"]].analysed_sst
points = (grid_x_old,grid_y_old)
grid = griddata(points, values, (grid_x_new, grid_y_new),method='nearest')

Scipy 内插网格数据的问题

Issues with Scipy interpolate griddata

python

grid

interpolation

scipy

netcdf