屏蔽插值 returns 常量值

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（没有任何意义）：

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

# 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

弄清楚这件事非常棘手。发生的情况是插值函数必须用 nans 填充，这样插值才有效，但随后用有限值替换剩余的 nans（例如来自整个 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

导致：

屏蔽插值 returns 常量值

Masked interpolation returns constant values

python

numpy

interpolation

rasterio