用SciPy进行规则网格空间插值

气象上，分析不同空间分辨率的网格数据时，常需要将它们统一插值到一种网格上。本文介绍用Python的科学分析库SciPy完成这一任务。

数据1：CERES的气溶胶光学厚度数据（time x lat x lon），空间分辨率1 x 1度。

数据2：ERA5气象场数据，空间分辨率0.25 x 0.25度。

目标：将CERES气溶胶光学厚度插值到ERA5的0.25 x 0.25度网格上。

动手干吧：

首先，调入所需的函数库，读入源数据（CERES光学厚度）及其经纬度坐标信息。

from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt
import netCDF4 as nc

CERES_name = "./CERES_SYN1deg-Month_Terra-Aqua-MODIS_Ed4.1_Subset_200003-202308.nc"
ERA5_name = "./ERA5_surface_hrly_20190701.nc"
F_CERES_in = nc.Dataset(CERES_name, "r")
F_ERA5_in = nc.Dataset(ERA5_name,"r")

AOD_550 = np.array(F_CERES_in.variables["ini_aod55_mon"][0,:,:])
lat_CERES = np.array(F_CERES_in.variables["lat"])
lon_CERES = np.array(F_CERES_in.variables["lon"])

接下来，让我们从ERA5数据获取目标网格的信息：

lat_ERA5 = np.array(F_ERA5_in.variables["latitude"])
lon_ERA5 = np.array(F_ERA5_in.variables["longitude"])

到了最关键的步骤了——进行插值。我们使用的是SciPy的interpolate子库的RegularGridInterpolator函数。这一函数是专门用来进行规则网格插值的，支持线性（linear）、临近（nearist）、三次样条插值（cubic）等。

这一插值要分两部走。第一步：根据源数据空间曲面建立插值函数；第二步：利用该函数和新网格信息，获取新网格上的值（这里需要先将1维经纬度转为与新网格对应的经度、纬度二维矩阵）。注意，np.meshgrid的indexing有两个选项：‘ij’和‘xy’，代表两种不同数据排列顺序，例如 lat x lon 大小是11 x 22，则indexing=‘ij’的结果是 11x22 大小的矩阵，而 indexing = ‘xy’的结果是 22x11大小的矩阵。

# Step 1:
interp = interpolate.RegularGridInterpolator((lat_CERES, lon_CERES), AOD_550, method='linear')

# Step 2: 
lat_ERA5_g, lon_ERA5_g  = np.meshgrid(lat_ERA5, lon_ERA5)
AOD_550_new = interp((LAT_ERA5_g, lon_ERA5_g))

至此，插值完成，AOD_550_new即为插值出来的新网格上的数据。

接着进行可视化，代码和插值前后的图如下：

（细心的朋友也许会发现，两组数据的纬度排列顺序是不一样的，CERES从南向北，而ERA5从北向南。但是这不影响插值过程，只是在可视化时需至北向南排列）

（细心的朋友还会发现，插值前后数据的区域范围不同，本例子用全球数据插值到中国区域，多出的数据并不影响插值效果。）

fig, ax = plt.subplots(2,1,layout='constrained')

extent = (min(lon_CERES), max(lon_CERES), min(lat_CERES), max(lat_CERES))
plot1 = ax[0].imshow(AOD_550[::-1,:],cmap='RdBu_r', vmin=0, vmax=1.2, extent=extent)
ax[0].set_xlabel('longitude')
ax[0].set_ylabel('latitude')
ax[0].set_title('AOD550 (Before Interpolation)')
fig.colorbar(plot1, shrink=0.7)

extent = (min(lon_ERA5), max(lon_ERA5), min(lat_ERA5), max(lat_ERA5))
plot2 = ax[1].imshow(AOD_550_new,cmap='RdBu_r', vmin=0, vmax=1.2, extent=extent)
ax[1].set_xlabel('longitude')
ax[1].set_ylabel('latitude')
ax[1].set_title("AOD550 (After Interpolation)")
fig.colorbar(plot2, shrink=0.7)
plt.show()