三维数组的Python二维插值

Question

我有一个全球降水数据，其维度为(time=410,lat=72,lon=144)。我希望将纬度内插到180，将纬度内插到360（即把数据从2.5度降到1度）。

在Matlab中，我曾经这样做。

%LAT,LON,XI,YI are old and new meshgridded lat-lons
for t = 1:size(precip,1)
    newPrecip(t,:,:)=interp2(LON,LAT,squeeze(precip(t,:,:)),XI,YI);
end

在python中我试过 interp2d, map_coordinates和 interp (basemap)函数，但没有成功。

下面是我的代码，其中有 map_coordinates这让我最接近得到一个结果。

new_indicies = np.mgrid[0:410, -89.5:89.5:180J, -179.5:179.5:360J]
newPrecip = ndimage.map_coordinates(precip, new_indicies, order=1)

新的数据确实有我想要的尺寸(410,180,360)，但它实际上并没有进行插值，而只是在新添加的网格(在原始数据周围)中填充了0个值。

我是python的新手，如果能帮助我解决这个问题，我将感激不尽。

Answer 1

我通常会使用 底图.

import mpl_toolkits.basemap as mp
Zg = mp.interp(dataIN,lonIN,latIN,lonOUT,latOUT,
               checkbounds=False, masked=False, order=1)

这里，lonIN和latIN是原始网格的1D经纬度（假设是常规网格），dataIn是原始数据的2D数组，lonOUT和latOUT是要插值到的2D网格，而Zg是新网格上的输出。输出网格可以使用。

x=np.arange(-180,180,1)
y=np.arange(-90,90,1)
lonOUT,latOUT=np.meshgrid(x,y)

希望能帮到你 干杯，Trond

Answer 2

我想看看 scipy.interpolate. 看起来像 网格数据 可能是你所需要的。

Answer 3

感谢Trond，我意识到我做错了什么。这是最后的代码，对我来说是可行的。希望对其他人有参考作用。这是我第一次在Stackoverflow上发帖，很高兴我的疑问得到了快速而正确的回答！我的疑问是什么？

#this code assumes a input data 'precip' of dimensions (410,72,144)
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import interp

LAT=np.arange(-89.5,90.5,1)
LON=np.arange(-179.5,180.5,1)
LON,LAT=np.meshgrid(LON,LAT)
lat=np.arange(-88.75,91.25,2.5)
lon=np.arange(-178.75,181.25,2.5)

newPrecip=np.zeros((410,180,360), dtype='float')
for i in range(410):
    newPrecip[i,:,:]=interp(np.squeeze(precip[i,:,:]),lon,lat,LON,LAT,order=1)

plt.figure(1)
plt.pcolor(lon,lat,precip.mean(axis=0))
plt.figure(2)
plt.pcolor(LON,LAT,newPrecip.mean(axis=0))

Answer 4

如果你想要普通Kriging插值，还有一个选择就是PyKrige库!