在 Numpy 或 Xarray 中对 3D 数组实现 1D 插值

Question

我一直在尝试找到一种利用 np.interp() 来插值 3D 数组最内层维度的方法：

假设我有一个形状为：(lat,lon,深度) = (3,3,4) 的 3D 数组作为玩具示例：

soil_temp = np.asarray([[[ 7.720276, 7.1702576, 4.821167, -1.7201233 ],
              [ 8.950958, 7.4095764, 1.5319519, -0.59176636],
              [ 4.1655273, 3.452118, 0.5239258, -3.0393372 ]],
             [[ 9.18335, 8.561096, 6.048462, -1.2128601 ],
              [ 8.009003, 7.6471863, 6.0042725, -0.17993164],
              [ 5.075775, 4.0643005, 0.55792236, -2.9087524 ]],
             [[13.139526, 12.393707, 9.50769, 0.5014343 ],
             [10.899994, 8.8767395, 1.443573, -0.96343994],
             [9.239685, 7.3506165, 0.4899292, -0.33395386]]])

我的目标是估计活性层厚度（温度达到 0 度时的深度）。我知道如何对单个数组执行此操作：

depths = np.asarray([3.5,17.5,64,194.5])
stemps = np.asarray([7.720276,7.1702576,4.821167,-1.7201233])

thaw_point=[0]

order = stemps.argsort() #sort into monotonically increasing values

y_data = stemps[order] #sort into monotonically increasing values
x_data = depths[order] #sort into monotonically increasing values

thaw_point=0

ALT = np.interp(thaw_point,y_data,x_data) 

这给了我 ALT = 160.183cm

但是，我一直很难弄清楚如何以计算有效的方式将 np.interp() 应用于 3D 数组。

我可以在循环结构中完成它，这对于这个小例子来说很好，但对于更大的数组来说就不行了：

def ALT_interp(data):
    thaw_point = 0
    array_shape = data.shape
    print(array_shape)
    deps = depths
    ALT_array=np.empty([array_shape[0],array_shape[1],1])
    for idx in range(0,array_shape[0]):
        for idx2 in range(0, array_shape[1]):
            stemps = data[idx,idx2,:]
            order = stemps.argsort()   
            y_data = stemps[order]
            x_data = depths[order]
            ALT = np.interp(thaw_point,y_data,x_data)
            ALT_array[idx,idx2,0] = ALT

    return ALT_array
    
depths = np.asarray([3.5,17.5,64,194.5])

soil_temp = np.asarray([[[ 7.720276, 7.1702576, 4.821167, -1.7201233 ],
              [ 8.950958, 7.4095764, 1.5319519, -0.59176636],
              [ 4.1655273, 3.452118, 0.5239258, -3.0393372 ]],
             [[ 9.18335, 8.561096, 6.048462, -1.2128601 ],
              [ 8.009003, 7.6471863, 6.0042725, -0.17993164],
              [ 5.075775, 4.0643005, 0.55792236, -2.9087524 ]],
             [[13.139526, 12.393707, 9.50769, 0.5014343 ],
             [10.899994, 8.8767395, 1.443573, -0.96343994],
             [9.239685, 7.3506165, 0.4899292, -0.33395386]]])
   
ALT_values = ALT_interp(soil_temp)

但是有没有一种计算效率更高的方法来做到这一点（因为我希望将其扩展到 1801x3600x4 数组）？

Answer 1

NumPy 或 SciPy 1D 插值器似乎都不是为此设计的。他们似乎期望 x 坐标是一维数组，但您希望 x 坐标是多维的。

一种方法是使用 SciPy 插值器之一（例如 CubicSpline），其中

depths

作为 x 数据，

soil_temp

作为 y 数据，然后求根。

但是，由于您似乎对线性插值感到满意，所以我会做更简单的事情 - 手动线性插值。看起来要在之间插值的点是每行的最后两个点，所以：

import numpy as np
from scipy import interpolate

depths = np.asarray([3.5,17.5,64,194.5])

soil_temp = np.asarray([[[ 7.720276, 7.1702576, 4.821167, -1.7201233 ],
              [ 8.950958, 7.4095764, 1.5319519, -0.59176636],
              [ 4.1655273, 3.452118, 0.5239258, -3.0393372 ]],
             [[ 9.18335, 8.561096, 6.048462, -1.2128601 ],
              [ 8.009003, 7.6471863, 6.0042725, -0.17993164],
              [ 5.075775, 4.0643005, 0.55792236, -2.9087524 ]],
             [[13.139526, 12.393707, 9.50769, 0.5014343 ],
             [10.899994, 8.8767395, 1.443573, -0.96343994],
             [9.239685, 7.3506165, 0.4899292, -0.33395386]]])

x1 = depths[-2]
x2 = depths[-1]
y1 = soil_temp[:, :, -2]
y2 = soil_temp[:, :, -1]

roots = x1 + (x2 - x1)/(y2 - y1) * (0 - y1)


roots2 = np.empty((3, 3))
for i in range(3):
    for j in range(3):
        stemps = soil_temp[i, j]
        order = stemps.argsort() #sort into monotonically increasing values
        y_data = stemps[order] #sort into monotonically increasing values
        x_data = depths[order] #sort into monotonically increasing values
        thaw_point=0
        roots2[i, j] = np.interp(thaw_point, y_data, x_data) 

np.testing.assert_allclose(roots, roots2)

唯一的失败是数据中没有零交叉的行。您的

np.interp

解决方案没有外推，因此它返回的根是观察到的最大深度 (194.5)，但真正的“活动层厚度”会高于该值。