Python GDAL 实现栅格叠加运算 - 支持不同行列数

要实现栅格的叠加运算，可以使用 gdal 库中的 RasterIO 函数和 numpy 库来处理栅格数据。以下是一个示例代码：

import gdal
import numpy as np

def raster_overlay(raster1, raster2, output_raster):
    # 打开栅格文件
    ds1 = gdal.Open(raster1, gdal.GA_ReadOnly)
    ds2 = gdal.Open(raster2, gdal.GA_ReadOnly)

    # 获取栅格文件的行列数和数据类型
    rows1 = ds1.RasterYSize
    cols1 = ds1.RasterXSize
    dtype1 = ds1.GetRasterBand(1).DataType

    rows2 = ds2.RasterYSize
    cols2 = ds2.RasterXSize
    dtype2 = ds2.GetRasterBand(1).DataType

    # 读取栅格数据
    data1 = ds1.GetRasterBand(1).ReadAsArray().astype(np.float32)
    data2 = ds2.GetRasterBand(1).ReadAsArray().astype(np.float32)

    # 创建输出栅格文件
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(output_raster, cols1, rows1, 1, dtype1)

    # 将栅格数据写入输出栅格文件
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(data1 + data2)

    # 设置输出栅格文件的地理参考信息
    out_ds.SetProjection(ds1.GetProjection())
    out_ds.SetGeoTransform(ds1.GetGeoTransform())

    # 关闭栅格文件
    ds1 = None
    ds2 = None
    out_ds = None

# 示例用法
raster1 = 'input1.tif'
raster2 = 'input2.tif'
output_raster = 'output.tif'

raster_overlay(raster1, raster2, output_raster)

此代码中，'raster1' 和 'raster2' 是输入栅格文件的路径，'output_raster' 是输出栅格文件的路径。代码首先打开输入栅格文件，然后获取其行列数和数据类型。接下来，使用 'ReadAsArray' 函数读取栅格数据，并将其转换为 'np.float32' 类型。然后，使用 'GetDriverByName' 函数创建输出栅格文件，并将栅格数据写入输出栅格文件。最后，设置输出栅格文件的地理参考信息，并关闭栅格文件。

需要注意的是，如果输入的栅格文件行列数不同，本例中的代码会按照第一个栅格文件的行列数进行处理。如果需要根据行列数较大的栅格文件进行处理，可以修改代码中的行列数变量。另外，栅格文件的投影信息和地理变换信息也会按照第一个栅格文件的信息设置，如果需要根据另一个栅格文件的信息设置，可以相应地修改代码。