3

Numpy 数组data有....数据。Numpy 数组z有距离。data 和 z 的形状相同,z 的每个点都是对应的 data 点被测量的距离。更复杂的是,用户将提供具有 3、4 或 5 维的 data/z 数组。

我想将数据内插到 1D numpy array 中的一组距离dists。由于数据结构的原因,插值轴总是从末尾开始的两个轴,即如果数组有3维,插值轴为0;如果数组有 4 个维度,插值轴是 1,等等。因为,AFAICT,所有 numpy/scipy 插值例程都希望在一维数组中给出原始距离,插值数据和 z 在 dists 上似乎是一个有点复杂的任务. 这就是我所拥有的:

def dist_interp(data, z, dists):
    # construct array to hold interpolation results
    num_dims = len(data.shape)
    interp_axis = num_dims-3
    interp_shape = list(data.shape)
    interp_shape[interp_axis] = dists.shape[0]
    interp_res = np.zeros(shape=interp_shape)
    # depending on usage, data could have between 3 and five dimensions.
    # add dims to generalize.  I hate doing it this way.  Must be
    # some other way.
    for n in range(num_dims, 5) :
        data = np.expand_dims(data, axis=0)
        z = np.expand_dims(z, axis=0)
        interp_res = np.expand_dims(interp_res, axis=0)
    for m in range(data.shape[0]):
        for l in range(data.shape[1]):
            for j in range(data.shape[3]):
                for i in range(data.shape[4]):
                    interp_res[m,l,:,j,i]=(
                                  np.interp(dists,z[m,l,:,j,i],
                                            data[m,l,:,j,i]))
    # now remove extra "wrapping" dimensions
    for n in range(0,5-num_dims):
        interp_res = interp_res[0]
    return(interp_res)

我认为这会起作用,但是添加和删除额外的“包装”虚拟尺寸非常不雅,并且会导致代码一点也不紧凑。有更好的想法吗?谢谢。

4

1 回答 1

4

作为这种情况的一般经验法则:

  1. 将你想要做的事情的轴移动到形状元组的末尾
  2. 将结果数组重塑为 2D
  3. 从这个 2D 数组创建新数组
  4. 撤消新阵列上的步骤 2 和 1

对于您的情况,它可能类似于:

# Create some random test data
axis = -2
shape = np.random.randint(10, size=(5,))
data = np.random.rand(*shape)
data = np.sort(data, axis=axis)
z = np.random.rand(*shape)
dists = np.linspace(0,1, num=100)

data = np.rollaxis(data, axis, data.ndim)
new_shape = data.shape
data = data.reshape(-1, data.shape[-1])
z = np.rollaxis(z, axis, z.ndim)
z = z.reshape(-1, z.shape[-1])

out = np.empty(z.shape[:1]+dists.shape)
for o, x, f in zip(out, data, z):
    o[:] = np.interp(dists, x, f)

out = out.reshape(new_shape[:-1]+dists.shape)
out = np.rollaxis(out, -1, axis)
于 2013-10-16T19:37:51.173 回答