利用python脚本，裁剪TIFF图（卫星地形数据）

import os
from osgeo import gdal

def crop_tif(input_tif, output_dir):
    src_ds = gdal.Open(input_tif)
    if src_ds is None:
        print("无法打开输入的TIF文件")
        return

    # 获取输入TIF的宽度和高度
    width = src_ds.RasterXSize
    height = src_ds.RasterYSize

    # 计算块的宽度和高度
    block_width = int(width / 4)
    block_height = int(height / 4)

    # 创建输出目录
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)

    # 裁剪为十六个小块
    for i in range(4):
        for j in range(4):
            output_tif = os.path.join(output_dir, f"block_{i}_{j}.tif")

            # 计算裁剪区域的范围
            xoff = i * block_width
            yoff = j * block_height

            # 创建输出dataset
            driver = gdal.GetDriverByName("GTiff")
            dst_ds = driver.Create(output_tif, block_width, block_height, src_ds.RasterCount, src_ds.GetRasterBand(1).DataType)

            # 设置裁剪区域
            dst_ds.SetGeoTransform((src_ds.GetGeoTransform()[0] + xoff * src_ds.GetGeoTransform()[1],
                                    src_ds.GetGeoTransform()[1],
                                    0,
                                    src_ds.GetGeoTransform()[3] + yoff * src_ds.GetGeoTransform()[5],
                                    0,
                                    src_ds.GetGeoTransform()[5]))

            # 获取颜色映射表
            color_table = src_ds.GetRasterBand(1).GetColorTable()

            # 设置颜色映射表
            for k in range(1, src_ds.RasterCount + 1):
                dst_band = dst_ds.GetRasterBand(k)
                dst_band.SetColorTable(color_table)

            # 执行裁剪
            for k in range(1, src_ds.RasterCount + 1):
                dst_band = dst_ds.GetRasterBand(k)
                dst_band.WriteArray(src_ds.GetRasterBand(k).ReadAsArray(xoff, yoff, block_width, block_height))

            # 设置输出dataset的地理信息
            dst_ds.SetProjection(src_ds.GetProjection())

            # 关闭dataset
            dst_ds = None

    # 关闭输入dataset
    src_ds = None

# 调用函数进行裁剪示例
input_tif = "E:\\QingDaoCityData\\ownData\\owmPhoto\\xxxxxxx10181717_L18\\chengyang18.tif"  # 输入TIF文件路径
output_dir = "E:\\QingDaoCityData\\ownData\\owmPhoto\\xxxxxxx10181717_L18\\cropping3_3"   # 输出目录
crop_tif(input_tif, output_dir)
print("裁剪完成")