Python如何实现二维数组插值？Matlab代码迁移numpy.interp报错怎么解决

错误根本原因 #

• numpy.interp仅支持1维线性插值，无法替代Matlab的interp2（2维插值函数），传入2D网格和2D数组时会触发维度不匹配错误
• 你调用np.interp时传入了5个参数，而该函数仅支持numpy.interp(x, xp, fp, left=None, right=None, period=None)的参数规则，API使用完全错误，是触发报错的直接原因
• 你的Python代码还存在几处语法/逻辑错误：
  • 循环变量range(1,b+1)中的b未定义，应该是bands
  • 计算ny数组时错误使用了old_size[0]，应该替换为old_size[1]
  • 未提前初始化new_imagen数组，直接赋值会触发索引错误

解决方案 #

使用scipy.interpolate.RegularGridInterpolator实现和Matlabinterp2一致的2D插值效果，该方法支持nearest、linear等插值方法，匹配你原代码的需求。

修正后的完整代码 #

import numpy as np
from osgeo import gdal
from scipy.interpolate import RegularGridInterpolator

source = "multi_rgb.tif"
new_multi = gdal.Open(source)

old_pixel = 4
new_pixel = 1
old_size = np.array([old_pixel,old_pixel]).astype(np.float32)
new_size = np.array([new_pixel,new_pixel]).astype(np.float32)

# 插值方法和原Matlab对应
metod = 'nearest'

f = new_multi.RasterYSize
c = new_multi.RasterXSize
bands = new_multi.RasterCount
o_s = old_size / 2
n_s = new_size / 2

# 修正ny计算错误
ox = np.arange(o_s[0], c*old_size[0] + o_s[0], old_size[0])
oy = np.arange(o_s[1], f*old_size[1] + o_s[1], old_size[1])
nx = np.arange(n_s[0], c*old_size[0] + n_s[0], new_size[0])
ny = np.arange(n_s[1], f*old_size[1] + n_s[1], new_size[1])

# 读取所有波段
band_list = []
for i in range(1, bands+1):
    band = new_multi.GetRasterBand(i).ReadAsArray().astype(np.uint8)
    band_list.append(band)
old_imagen = np.stack(band_list, axis=2)

# 提前初始化新数组
new_imagen = np.zeros((len(ny), len(nx), bands), dtype=np.uint8)

# 逐波段插值
for i in range(bands):
    print(f"处理第 {i+1} 个波段")
    # 构建2D插值器
    interp_func = RegularGridInterpolator((oy, ox), old_imagen[:,:,i], method=metod, bounds_error=False, fill_value=0)
    # 生成新网格的坐标点
    Nx, Ny = np.meshgrid(nx, ny)
    points = np.stack((Ny.flatten(), Nx.flatten()), axis=1)
    # 插值后还原维度
    new_imagen[:,:,i] = interp_func(points).reshape(len(ny), len(nx)).astype(np.uint8)

优化方案 #

如果你的场景只是做影像重采样，更简单的方案是直接用GDAL自带的重采样接口，不需要手动实现插值，性能和稳定性更高：

# GDAL重采样示例
output_path = "resampled_multi_rgb.tif"
gdal.Warp(output_path, new_multi, xRes=new_pixel, yRes=new_pixel, resampleAlg=gdal.GRA_NearestNeighbour)

内容的提问来源于stack exchange，提问作者Daniel Ariza