使用python gdal将影像重采样之取和值(sum)

我们都知道gdal为我们提供了多种重采样的方法，这里我们不再叙述，有兴趣可以自己去看重采样方法
下面我们来介绍怎样进行取和值重采样，针对少数特别的需求。

思路

1. 把整个图像看成一张网格，一个个的像元，我们要改变格子的大小，将小的格子合并成大的或者将大的分割成小的，我们这里讲合并成大的，因为我们要取和值；
   
2. 使用gdal将图像的矩阵数组读出，取得这些格子的数据，一个矩阵，然后进行计算；
3. 使用numpy，将矩阵进行分割，分割成我们想要的想要大小，对应矩阵的大小。
   均分：

np.hsplit(b1, 3)  # 参数 3 指的是将 b1 矩阵均分为3份，矩阵第 1 维的列数必须是3的倍数，否则报错。
np.hsplit(b2, 2)  # 参数 2 指的是将 b2 矩阵均分为2份。

非均分：

np.split(b1, 3, axis=1)
np.split(b2, 2, axis=1)

这里我们使用非均分，因为图像有可能不是刚好的矩形，可能会又一点误差，但是不影响我们使用，效果已测试，几乎看不到误差。

4. 代码

# 重采样栅格sum
def resample_tif(in_raster, out_raster,out_size):
    start=time.clock()
    ds = gdal.Open(in_raster)
    in_band = ds.GetRasterBand(1)
    x_size = in_band.XSize
    y_size = in_band.YSize
    in_geotrans = list(ds.GetGeoTransform())
    x_res = list(ds.GetGeoTransform())[1]
    y_res = list(ds.GetGeoTransform())[5]
    width = int(abs(x_size * x_res / out_size)) + 1
    height = int(abs(y_size * y_res / out_size)) + 1
    raster_data = in_band.ReadAsArray()
    print("width:" + str(width))
    print("height:" + str(height))
    # 创建空矩阵来放置数据
    result_data = np.empty((height, width), np.float32)
    # 水平分割矩阵
    arr_split_h = np.array_split(raster_data, width, axis=1)

    for i in range(len(arr_split_h)):
        #     垂直分割矩阵
        arr_split_v = np.array_split(arr_split_h[i], height, axis=0)
        for j in range(len(arr_split_v)):
            # 求和
            result_data[j][i] = np.sum(arr_split_v[j])

    out_geotrans = tuple([in_geotrans[0], out_size, 0.0, in_geotrans[3], 0.0, -out_size])
    out_ds = ds.GetDriver().Create(out_raster, width, height, 1,
                                   in_band.DataType)  # 创建一个构建重采样影像的句柄
    out_ds.SetProjection(ds.GetProjection())  # 设置投影信息
    out_ds.SetGeoTransform(out_geotrans)  # 设置地理变换信息
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(result_data)
    out_band.SetNoDataValue(-999)
    out_band.FlushCache()

    end=time.clock()
    print("耗时："+str(end-start)+"秒")
    return