将一副遥感影像裁剪为多幅图像

参考：https://blog.csdn.net/zsc201825/article/details/89359995

以Landsat8某景影像为例：原图（底图）是假彩色的，上面单波段图像为裁剪的9景影像。

代码如下：

from osgeo import gdal
import math

# 读取要裁剪的原始遥感影像
in_ds = gdal.Open("D:\ProfessionalProfile\DataRelevant\L134036_20170808.tif")
print("open tif file succeed")
width = in_ds.RasterXSize  # 获取数据宽度
height = in_ds.RasterYSize  # 获取数据高度
outbandsize = in_ds.RasterCount  # 获取数据波段数
# outbandsize = in_ds.GetRasterBand
im_geotrans = in_ds.GetGeoTransform()  # 获取仿射矩阵信息
im_proj = in_ds.GetProjection()  # 获取投影信息
datatype = in_ds.GetRasterBand(1).DataType
im_data = in_ds.ReadAsArray()  # 获取数据

# 读取原图中的每个波段
for i in range(1,outbandsize+1):
    locals()['in_band' + str(i)] =in_ds.GetRasterBand(i)

# 定义切图的起始点坐标，自己选定即可。
offset_x = 0
offset_y = 0

# 定义切图的大小（矩形框）
block_xsize = 2000  # 行
block_ysize = 2000  # 列

k = 0
# ‘//’: 除法结果向下取整；‘/’: 除法结果为浮点数。
for i in range( width // block_xsize):
    for j in range( height // block_xsize):

        # # 从每个波段中切需要的矩形框内的数据(注意读取的矩形框不能超过原图大小)
        for x in range(1, outbandsize + 1):
            locals()['out_band' + str(x)] = locals()['in_band'+str(x)].ReadAsArray(i * block_xsize, j * block_xsize, block_xsize, block_ysize)

        # print(out_band7)
        k += 1

        # 获取Tif的驱动，为创建切出来的图文件做准备
        gtif_driver = gdal.GetDriverByName("GTiff")

        # 创建切出来的要存的文件（3代表3个不都按，最后一个参数为数据类型，跟原文件一致）
        out_ds = gtif_driver.Create(str(k) + 'Clip.tif', block_xsize, block_ysize, outbandsize, datatype)
        print("create new tif file succeed")

        # 获取原图的原点坐标信息
        ori_transform = in_ds.GetGeoTransform()
        if ori_transform:
            print(ori_transform)
            print("Origin = ({}, {})".format(ori_transform[0], ori_transform[3]))
            print("Pixel Size = ({}, {})".format(ori_transform[1], ori_transform[5]))

        # 读取原图仿射变换参数值
        top_left_x = ori_transform[0]  # 左上角x坐标
        w_e_pixel_resolution = ori_transform[1]  # 东西方向像素分辨率
        top_left_y = ori_transform[3]  # 左上角y坐标
        n_s_pixel_resolution = ori_transform[5]  # 南北方向像素分辨率

        # 根据反射变换参数计算新图的原点坐标
        top_left_x = top_left_x + i * block_xsize * w_e_pixel_resolution
        top_left_y = top_left_y + j * block_xsize * n_s_pixel_resolution

        # 将计算后的值组装为一个元组，以方便设置
        dst_transform = (top_left_x, ori_transform[1], ori_transform[2], top_left_y, ori_transform[4], ori_transform[5])

        # 设置裁剪出来图的原点坐标
        out_ds.SetGeoTransform(dst_transform)

        # 设置SRS属性（投影信息）
        out_ds.SetProjection(in_ds.GetProjection())

        # # 写入目标文件
        for y in range(1, outbandsize + 1):
            out_ds.GetRasterBand(y).WriteArray( locals()['out_band' + str(y)] )

        # 将缓存写入磁盘
        out_ds.FlushCache()
        print("FlushCache succeed")

        # del out_ds

        print("End!")