Python中10m分辨率浮点型栅格5×5邻域均值平滑报错求助

嗨，我来帮你排查这个报错的问题～

首先看你抛出的错误TypeError: 'float' object cannot be interpreted as an integer，核心原因很明确：你计算新数组形状shape = (h/bh, w/bw, bh, bw)的时候用了普通除法/，得到了浮点数（比如106/5=21.2、144/5=28.8），但np.lib.stride_tricks.as_strided要求形状的每个元素必须是整数——它需要把原数组划分成规整的块，浮点数尺寸根本没法划分。

接下来给你一步步的修复方案，还有更简单的替代写法：

一、直接修复原代码的问题 #

1. 修正数组形状计算 #

原数组的高106、宽144都不能被5整除，你可以选择裁剪原数组到能被5整除的最大尺寸，这样就能得到整数形状：

bh, bw = 5,5
# 计算裁剪后的最大规整尺寸
h_crop = (h // bh) * bh
w_crop = (w // bw) * bw
# 裁剪数组
ds_array_cropped = ds_array[:h_crop, :w_crop]
# 重新计算整数形状
shape = (h_crop // bh, w_crop // bw, bh, bw)

2. 修复Python3兼容性问题 #

你的代码里用了Python2的xrange，但Python3里已经统一成range了；另外time.clock()在Python3.3+已经被废弃，换成time.perf_counter()计时更准确：

# 替换原来的计时代码
start_time_script = time.perf_counter()
# ...
end_time_blocks = time.perf_counter() - start_time_blocks

3. 完善空间参考设置 #

你注释掉了投影设置，其实可以直接复用原栅格的投影，不用手动输入EPSG：

# 直接读取原栅格的投影信息
srs = osr.SpatialReference()
srs.ImportFromWkt(ds.GetProjectionRef())

二、更简单的替代方案：用Scipy做滑动窗口均值 #

手动分块写循环不仅容易出错，效率也不高。scipy.ndimage.uniform_filter可以直接帮你计算指定窗口的均值，还能处理任意尺寸的数组，甚至可以选择边界的处理方式：

import time
import glob
import os
import gdal
import osr
import numpy as np
from scipy.ndimage import uniform_filter  # 导入这个库

start_time_script = time.perf_counter()
path_ras=r'D:\Firm_SM\F1A/'

for rasterfile in glob.glob(os.path.join(path_ras,'*.tif')):
    rasterfile_name = rasterfile[rasterfile.find('IMG'):rasterfile.find('.tif')]
    print(f'Processing: {rasterfile_name}')
    
    ds = gdal.Open(rasterfile,gdal.GA_ReadOnly)
    ds_xform = ds.GetGeoTransform()
    ds_driver = gdal.GetDriverByName('Gtiff')
    # 复用原栅格投影
    srs = osr.SpatialReference()
    srs.ImportFromWkt(ds.GetProjectionRef())
    
    ds_array = ds.ReadAsArray()
    h,w = ds_array.shape
    print(f'Original array size: {h} x {w}')
    
    # 计算5x5窗口均值，mode选nearest避免边界值异常
    smoothed_array = uniform_filter(ds_array, size=(5,5), mode='nearest')
    # 如果需要降采样到50x50米（原10x10，5倍分辨率）
    downsampled_array = smoothed_array[::5, ::5]
    new_h, new_w = downsampled_array.shape
    
    # 创建输出栅格
    resized_array = ds_driver.Create(
        f'{rasterfile_name}_resized_to_50m.tif',
        new_w, new_h, 1, gdal.GDT_Float32
    )
    # 更新地理变换：分辨率乘以5
    new_xform = (
        ds_xform[0], ds_xform[1]*5, ds_xform[2],
        ds_xform[3], ds_xform[4], ds_xform[5]*5
    )
    resized_array.SetGeoTransform(new_xform)
    resized_array.SetProjection(srs.ExportToWkt())
    
    # 写入数组
    band = resized_array.GetRasterBand(1)
    band.WriteArray(downsampled_array)
    # 务必FlushCache确保写入磁盘
    resized_array.FlushCache()
    
    print(f'Finished processing {rasterfile_name}\n')

end_time = time.perf_counter() - start_time_script
print(f'Program ran for: {end_time:.2f} seconds')

这个写法不仅代码更简洁，还避免了手动分块的各种坑，效率也更高。

内容的提问来源于stack exchange，提问作者user286076