【Python&RS】基于GDAL遥感影像分幅裁剪（固定尺寸）

        之前分享过一篇分幅裁剪的文章：【Python&RS】基于GDAL遥感影像分幅裁剪，只不过这篇文章当时编写的逻辑是自己输入需要裁剪多少行多少列，由于大家可能并没有直观地希望自己裁剪多少行列，所以非常局限。今天跟大家分享一下使用固定尺寸对遥感影像进行分幅裁剪，即每张裁剪的影像都是一样大的。

原创作者：RS迷途小书童

博客地址：https://blog.csdn.net/m0_56729804?type=blog

1. 代码逻辑

        逻辑其实很简单，就是利用输入的尺寸计算出需要裁剪多少行多少列。然后按照之前文章的方法进行裁剪即可。当然这里要注意一下边缘处的处理。

2. 代码主函数

        由于和之前发的文章原理差不多，所以我这里就不做过多的解释了，代码中都有详细的注释，如果大家有什么问题可以给我留言。

ds = gdal.Open(filepath)  # 打开数据集dataset
ds_width = ds.RasterXSize  # 获取数据宽度
ds_height = ds.RasterYSize  # 获取数据高度
ds_bands = ds.RasterCount  # 获取波段数
raw = int(ds_height / size) + 1
col = int(ds_width / size) + 1
print("分割后影像行数：", raw)
print("分割后影像列数：", col)
print("分割后影像数为：", raw * col)
for j in range(0, raw):
    print("正在裁剪第%s行......" % (j + 1))
    for k in range(0, col):
        raw_frame = size
        # 计算每幅图像的高度
        col_frame = size
        # 计算每幅图像的宽度
        left_x = j * raw_frame
        left_y = k * col_frame
        # 计算当前影像的左上角像素坐标
        raw_frame = min(ds_height-left_x, raw_frame)
        col_frame = min(ds_width-left_y, col_frame)
        # 防止幅宽超过整幅影像
        driver = gdal.GetDriverByName('GTiff')  # 载入数据驱动，用于存储内存中的数组
        ds_result = driver.Create(out_path+"%s_%s.tif" % (j+1, k+1), col_frame, raw_frame,
                                  bands=ds_bands, eType=gdal.GDT_Byte)  # 创建空tif
        ds_geo = ds.GetGeoTransform()  # 获取仿射地理变换参数
        top_left_x = ds_geo[0]  # 原始影像左上角x的投影坐标
        top_left_y = ds_geo[3]  # 原始影像左上角y投影坐标
        top_left_x = top_left_x + left_y * ds_geo[1]
        top_left_y = top_left_y + left_x * ds_geo[5]
        # 计算得到当前影像的左上角投影坐标
        ds_geo = (top_left_x, ds_geo[1], ds_geo[2], top_left_y, ds_geo[4], ds_geo[5])
        # 新影像的仿射地理变换参数
        ds_result.SetGeoTransform(ds_geo)  # 导入仿射地理变换参数
        ds_result.SetProjection(ds.GetProjection())  # 导入投影信息
        array_band = []
        for i in range(1, ds_bands+1):
            array_band = ds.GetRasterBand(i).ReadAsArray(left_y, left_x, col_frame, raw_frame).astype(np.float64)
            # 根据左上角的像素坐标和幅宽读取指定区域内的数据
            ds_result.GetRasterBand(i).SetNoDataValue(0)  # 将无效值设为0
            ds_result.GetRasterBand(i).WriteArray(array_band)  # 将每个波段写入新的文件中
        ds_result = None
        del ds_result

3. 总结

        今天主要分享的是遥感影像的分幅裁剪，大家可以用这段代码减少数据量，也可以用它制作样本集。同时结合之前发的那篇分幅裁剪的文章，基本上概括了所有的分幅裁剪的情况，友友们也可以根据自己的需求修改一下代码。如果大家在学习Python或者RS时有什么问题，可以随时留言交流！如果大家对批量处理有兴趣同样可以留言给博主，博主会分享相关代码以供学习！