numpy 移动窗口百分比覆盖

Question

我有一个分类栅格，我正在读入一个 numpy 数组。 (n 类)

我想使用 2d 移动窗口（例如 3 x 3）创建一个 n 维向量，用于存储窗口内每个类的 %cover。因为栅格很大，所以存储这些信息以免每次都重新计算它会很有用......因此我认为最好的解决方案是创建一个 3d 数组作为矢量。将根据这些 %/count 值创建一个新栅格。

我的想法是：

1) 创建一个 3d 数组 n+1 'bands'

2) 波段 1 = 原始分类栅格。彼此的“波段”值 = 计算窗口内某个值的单元格数（即每类一个波段）....例如：

[[2 0 1 2 1]
 [2 0 2 0 0]
 [0 1 1 2 1]
 [0 2 2 1 1]
 [0 1 2 1 1]]
[[2 2 3 2 2]
 [3 3 3 2 2]
 [3 3 2 2 2]
 [3 3 0 0 0]
 [2 2 0 0 0]]
[[0 1 1 2 1]
 [1 3 3 4 2]
 [1 2 3 4 3]
 [2 3 5 6 5]
 [1 1 3 4 4]]
[[2 3 2 2 1]
 [2 3 3 3 2]
 [2 4 4 3 1]
 [1 3 5 3 1]
 [1 3 3 2 0]]

4) 将这些波段读取到 vrt 中，因此只需创建一次 ...并且可以读取更多模块。

问题：在窗口内“计数”最有效的“移动窗口”方法是什么？

目前 - 我正在尝试，但使用以下代码失败：

def lcc_binary_vrt(raster, dim, bands):
    footprint = np.zeros(shape = (dim,dim), dtype = int)+1
    g = gdal.Open(raster)
    data = gdal_array.DatasetReadAsArray(g)

    #loop through the band values
    for i in bands:   
        print i
        # create a duplicate '0' array of the raster
        a_band = data*0
        # we create the binary dataset for the band        
        a_band = np.where(data == i, 1, a_band)
        count_a_band_fname = raster[:-4] + '_' + str(i) + '.tif'        
        # run the moving window (footprint) accross the band to create a 'count'
        count_a_band = ndimage.generic_filter(a_band, np.count_nonzero(x), footprint=footprint, mode = 'constant')
        geoTiff.create(count_a_band_fname, g, data, count_a_band, gdal.GDT_Byte, np.nan)

非常感谢任何建议。

贝基

Answer 1

我对空间科学的东西一无所知，所以我只关注主要问题:)

在窗口内“计数”最有效的“移动窗口”方法是什么？

使用 Numpy 进行移动窗口统计的常用方法是使用 numpy.lib.stride_tricks.as_strided，例如参见 this answer。基本上，这个想法是创建一个包含所有窗口的数组，而不会增加内存使用量：

from numpy.lib.stride_tricks import as_strided

...

m, n = a_band.shape
newshape = (m-dim+1, n-dim+1, dim, dim)
newstrides = a_band.strides * 2  # strides is a tuple
count_a_band = as_strided(ar, newshape, newstrides).sum(axis=(2,3))

但是，对于您的用例，此方法效率低下，因为您一遍又一遍地对相同的数字求和，尤其是在窗口大小增加的情况下。更好的方法是使用cumsum 技巧，例如this answer：

def windowed_sum_1d(ar, ws, axis=None):

    if axis is None:
        ar = ar.ravel()
    else:
        ar = np.swapaxes(ar, axis, 0)

    ans = np.cumsum(ar, axis=0)
    ans[ws:] = ans[ws:] - ans[:-ws]

    ans = ans[ws-1:]

    if axis:
        ans = np.swapaxes(ans, 0, axis)

    return ans


def windowed_sum(ar, ws):
    for axis in range(ar.ndim):
        ar = windowed_sum_1d(ar, ws, axis)
    return ar

...

count_a_band = windowed_sum(a_band, dim)

请注意，在上面的两个代码中，处理边缘情况会很乏味。幸运的是，有一种简单的方法可以包含这些代码并获得与第二个代码相同的效率：

count_a_band = ndimage.uniform_filter(a_band, size=dim, mode='constant') * dim**2

虽然与您已经拥有的非常相似，但这会更快！缺点是您可能需要舍入为整数以消除浮点舍入错误。

最后，您的代码

# create a duplicate '0' array of the raster
a_band = data*0
# we create the binary dataset for the band        
a_band = np.where(data == i, 1, a_band)

有点多余：你可以使用a_band = (data == i)。