python实现栅格滑动窗口统计值

任务目标：在25km的栅格内统计90m栅格的均值

任务转化：构建一个滑动窗口统计均值

主要代码如下：

def zonal_statisticDEM(zonal_tif, data_tif, out_tif):
    zonal = gdal.Open(zonal_tif)
    if zonal is None:
        print('打开数据' + zonal_tif + '失败！')
        sys.exit(1)
    # 读取栅格数据集的x方向像素数，y方向像素数
    zonal_cols = zonal.RasterXSize
    zonal_rows = zonal.RasterYSize

    data = gdal.Open(data_tif)

    if data is None:
        print('打开数据' + data_tif + '失败！')
        sys.exit(1)
    # 读取栅格数据集的x方向像素数，y方向像素数
    data_cols = data.RasterXSize
    data_rows = data.RasterYSize
    data_band = data.GetRasterBand(1)
    data_img = data_band.ReadAsArray(0, 0, data_cols, data_rows)

    d_col = data_cols / zonal_cols
    d_row = data_rows / zonal_rows

    outimg = np.zeros((zonal_rows, zonal_cols))
    for i in range(zonal_rows):
        for j in range(zonal_cols):
            counts = data_img[int(i * d_row):int((1 + i) * d_row), int(j * d_col):int((1 + j) * d_col)]
            # outimg[i, j] = np.std(counts)
            outimg[i, j] = np.average(counts)
    outimg = outimg.astype(np.float)
    # 创建并输出图像
    driver = zonal.GetDriver()
    outimage = driver.Create(out_tif, zonal_cols, zonal_rows, 1, gdal.GDT_Float32)
    geotransform = zonal.GetGeoTransform()
    im_proj = zonal.GetProjection()
    outimage.SetGeoTransform(geotransform)
    outimage.SetProjection(im_proj)
    # 在写入之前，还需要先引入波段对象
    outBand = outimage.GetRasterBand(1)
    # 波段对象支持直接写入矩阵，两个参数分别为x向偏移和y向偏移
    outBand.WriteArray(outimg, 0, 0)
    # 设定NoData值
    outBand.SetNoDataValue(255)
    outBand.FlushCache()
    outBand.GetStatistics(0, 1)
    outimage = None
    print('数据处理成功')