屏蔽插值返回常量值

Question

我想沿第一维插入 3D 数组。

就数据而言，这意味着我想在地理值中插入缺失时间，换句话说，稍微平滑一下这个动画：

我通过调用执行此操作：

new = ma.apply_along_axis(func1d=masked_interpolation, axis=0, arr=dst_data, x=missing_bands, xp=known_bands)

其中插值函数如下：

def masked_interpolation(data, x, xp, propagate_mask=True):
    import math
    import numpy as np
    import numpy.ma as ma

    # The x-coordinates (missing times) at which to evaluate the interpolated values.
    assert len(x) >= 1

    # The x-coordinates (existing times) of the data points (where returns a tuple because each element of the tuple refers to a dimension.)
    assert len(xp) >= 2

    # The y-coordinates (value at existing times) of the data points, that is the valid entries
    fp = np.take(data, xp)
    assert len(fp) >= 2

    # Returns the one-dimensional piecewise linear interpolant to a function with given discrete data points (xp, fp), evaluated at x.
    new_y = np.interp(x, xp, fp.filled(np.nan))

    # interpolate mask & apply to interpolated data
    if propagate_mask:
        new_mask = data.mask[:]
        new_mask[new_mask]  = 1
        new_mask[~new_mask] = 0
        # the mask y values at existing times
        new_fp = np.take(new_mask, xp)
        new_mask = np.interp(x, xp, new_fp)
        new_y = np.ma.masked_array(new_y, new_mask > 0.5)

    print(new_y) # ----> that seems legit
    data[x] = new_y # ----> from here it goes wrong
    return data

打印时new_y，插值看起来一致（分布在 [0,1] 区间，这是我想要的）。但是，当我打印最终输出（new 数组）时，它肯定更平滑（更多条带），但所有非屏蔽值都更改为 -0.1（没有任何意义）：

< a href="https://i.sstatic.net/TNiib.gif" rel="nofollow noreferrer">

将其写入光栅文件的代码是：

# Writing the new raster
meta = source.meta
meta.update({'count' : dst_shape[0] })
meta.update({'nodata' : source.nodata})
meta.update(fill_value = source.nodata)
assert new.shape == (meta['count'],meta['height'],meta['width'])
with rasterio.open(outputFile, "w", **meta) as dst:
    dst.write(new.filled(fill_value=source.nodata))

Answer 1

这很难弄清楚。所发生的情况是，插值函数必须用 nan 填充，以便插值起作用，但然后用有限值替换剩余的 nan（例如，来自整个 fp 向量为 nan 时）。然后应用插值掩码无论如何都会隐藏这些值。事情是这样的：

def masked_interpolation(data, x, xp, propagate_mask=True):
    import math
    import numpy as np
    import numpy.ma as ma

    # The x-coordinates (missing times) at which to evaluate the interpolated values.
    assert len(x) >= 1
    # The x-coordinates (existing times) of the data points (where returns a tuple because each element of the tuple refers to a dimension.)
    assert len(xp) >= 2
    # The y-coordinates (value at existing times) of the data points, that is the valid entries
    fp = np.take(data, xp)
    assert len(fp) >= 2

    # Returns the one-dimensional piecewise linear interpolant to a function with given discrete data points (xp, fp), evaluated at x.
    new_y = np.interp(x, xp, fp.filled(np.nan))
    np.nan_to_num(new_y, copy=False)

    # interpolate mask & apply to interpolated data
    if propagate_mask:
        new_mask = data.mask[:]
        new_mask[new_mask]  = 1
        new_mask[~new_mask] = 0
        # the mask y values at existing times
        new_fp = np.take(new_mask, xp)
        new_mask = np.interp(x, xp, new_fp)
        new_y = np.ma.masked_array(new_y, new_mask > 0.5)

    data[x] = new_y
    return data

结果是：