如何基于压力数组对不含压力信息的数组进行插值

Question

假设我有一个变量x，形状为(Time: 127, bottom_top: 58, south_north: 76, west_east: 96)，其中bottom_top没有明确包含原始压力的值。也就是说，原来的压力水平来自一个新的变量p，其形状与(Time: 127, bottom_top: 58, south_north: 76, west_east: 96)相同

我想做的是在bottom_top方向x进行插值，从原始压力p到新的压力水平p_new。例如

x[0,:,1,1].values =

array([404.13263, 404.13217, 404.13174, 404.11353, 404.11615, 404.1033 ,
       404.0999 , 404.0954 , 404.0902 , 404.07092, 404.05087, 404.04208,
       404.0345 , 404.01535, 403.98822, 403.95865, 403.92953, 403.90146,
       403.88913, 403.88214, 403.87677, 403.8698 , 403.86185, 403.8513 ,
       403.83362, 403.8124 , 403.791  , 403.76935, 403.7525 , 403.72525,
       403.6679 , 403.61584, 403.52277, 403.44308, 403.36752, 403.2787 ,
       403.17352, 403.0464 , 402.89468, 402.69656, 402.52823, 402.40897,
       402.3129 , 402.2052 , 402.11267, 402.06073, 402.03674, 401.99417,
       401.81857, 401.57474, 401.34845, 401.14847, 400.97156, 400.95645,
       400.92035, 400.7942 , 400.16495, 398.33502], dtype=float32)

p[0,:,1,1].values =

array([86983.83  , 86819.19  , 86656.23  , 86487.69  , 86326.62  , 86120.52  ,
       85831.95  , 85461.48  , 85050.03  , 84598.    , 84103.71  , 83613.43  ,
       83075.71  , 82500.    , 81923.46  , 81346.94  , 80730.56  , 80030.49  ,
       79289.44  , 78507.66  , 77642.914 , 76696.05  , 75707.52  , 74719.38  ,
       73690.03  , 72619.51  , 71549.4   , 70438.445 , 69286.41  , 68093.85  ,
       66819.19  , 65503.797 , 64148.29  , 62711.4   , 61234.523 , 59716.75  ,
       58158.11  , 56557.582 , 54874.125 , 53109.348 , 51303.953 , 49458.043 ,
       47530.21  , 45519.61  , 43423.426 , 41243.23  , 38992.438 , 36712.22  ,
       34389.49  , 31879.646 , 29135.926 , 26400.299 , 23733.72  , 20932.85  ,
       18026.195 , 15077.396 , 11595.89  ,  7294.6274], dtype=float32)

p_new =

array([1.00758e+05, 9.00880e+04, 8.03100e+04, 7.13940e+04, 6.32770e+04,
       5.59420e+04, 4.93250e+04, 4.33770e+04, 3.80270e+04, 3.32140e+04,
       2.89180e+04, 2.50580e+04, 2.16130e+04, 1.85420e+04, 1.58370e+04,
       1.34660e+04, 1.14090e+04, 9.64300e+03, 8.15100e+03, 6.90200e+03,
       5.86000e+03, 4.98700e+03, 4.25500e+03, 3.63600e+03, 3.11200e+03,
       2.66800e+03, 2.29100e+03, 1.97200e+03, 1.69900e+03, 1.46600e+03,
       1.26700e+03, 1.09600e+03, 9.49000e+02, 8.23000e+02, 7.15000e+02,
       6.21000e+02, 5.40000e+02, 4.70000e+02, 4.10000e+02, 3.58000e+02,
       3.12000e+02, 2.73000e+02, 2.40000e+02, 2.10000e+02, 1.85000e+02,
       1.63000e+02, 1.43000e+02, 1.26000e+02, 1.11000e+02, 9.80000e+01,
       8.70000e+01, 7.70000e+01, 6.70000e+01, 5.90000e+01, 5.20000e+01,
       4.60000e+01, 4.00000e+01, 3.50000e+01, 3.10000e+01, 2.70000e+01,
       2.40000e+01, 2.10000e+01, 1.80000e+01, 1.60000e+01, 1.40000e+01,
       1.20000e+01, 1.00000e+01, 9.00000e+00, 8.00000e+00, 7.00000e+00,
       6.00000e+00])

有什么方法可以直接插值吗

x_new = interpolation(x,p,p_new)

提前谢谢

Answer 1

我需要沿着不同的轴插值相当频繁。我使用以下方法作为一般解决方案。有更简单的具体解决方案。e、 g使用嵌套for循环沿每个轴迭代。我在numpy实用程序文件中有iter\u1d函数。我不确定这是不是最好的方法，但对我来说效果很好

# A general solution

import numpy as np
from scipy.interpolate import interp1d

# I've found this useful but never know whether there's a better way
def iter_1d( shape, axis=-1 ):
    """ Iterator for all 1d slices along 'axis' in an array of 'shape'.
        It returns a generic selection, not a selection of a specific array 
        Usage:  for s in iterate_1d(arr.shape, 2):
                    print(arr[s], b[s], c[s]) 
    """

    while axis<0:
        axis+=len(shape)
    limits = list(shape)        # [ 127, 58, 76, 96 ]
    current = [0]*len(shape)    # [   0,  0,  0,  0 ]
    current[axis] = slice(None) # [   0,  :,  0,  0 ]
    target = list(range(len(shape)-1, -1, -1))  # Reversed array of axis pointers
    # e.g 4d case:  target = [ 3, 2, 1, 0 ]
    target.remove(axis)  # Remove the 'axis' index from this list
    # e,g, case: axis = 1 from above target = [ 3, 2, 0 ]

    def incr_current():
        """ Increments the 'right hand' index in current.  
            If this overflows, increnent the one before etc...
            Returns True until axis 0 overflows, then returns False
        """
        nonlocal current
        for ix in target:
            current[ix] += 1
            if current[ix] < limits[ix]: return True
            current[ix] = 0
        return False

    cont=True
    while cont:
        yield tuple(current)
        cont = incr_current()

# Thos can iterate 1d slices along any axis.  The question is for axis 1
def do_interpolate( x_new, x, y, axis=-1 ):
    """ Will interpolate all 1D slices of the axis specified."""
    result = np.empty_like(x_new)

    for sel in iter_1d( x.shape, axis ):
        do_interp = interp1d(x[sel], y[sel], fill_value='extrapolate')
        result[sel] = do_interp(x_new[sel])
    return result

x_new = do_interpolate( p_new, p, x, axis = 1 )

插值函数都是从一个一维自变量开始工作的，因此每一组p与一组不同的x相关联，必须为每个切片生成插值对象

对于问题中的特定需求，简单地使用循环进行迭代可能更容易

# Specific solution
def interpolate( x, p, p_new):
    result = np.empty_like(p_new)

    for t in range(127):
        for ns in range(76):
            for ew in range(96):
                sel = np.s_[t, :, ns, ew]
                do_interp = interp1d(p[sel], x[sel], fill_value='extrapolate')
                result[sel] = do_interp(p_new[sel])
    return result

x_new = interpolate( x, p, p_new )