GDAL重采样与裁剪图像示例

        GDAL重采样，可以通过写文件时改变图像尺寸和geo_transformes的分辨率信息实现。核心代码示例如下：

in_ds = gdal.Open(fi, gdal.GA_ReadOnly)
geotrans = in_ds.GetGeoTransform()
geotrans = list(geotrans)
geoproj = in_ds.GetProjection()
# 更新为采样后分辨率
geotrans[1] = tar_resolution  #
geotrans[5] = -tar_resolution  # 
in_band = in_ds.GetRasterBand(2)
xsize = in_band.XSize
ysize = in_band.YSize
# 重新计算尺寸
ysize_t = int(ysize * img_resolution / tar_resolution)
xsize_t = int(xsize * img_resolution / tar_resolution)
# 按新的尺寸写
out_ds = in_ds.GetDriver().Create(fo_, xsize_t, ysize_t, 1,
                                          in_band.DataType)  # 创建一个构建重采样影像的句柄
out_ds.SetProjection(in_ds.GetProjection())  # 设置投影信息
# 更新geotrans
geotrans = tuple(geotrans)

        GDAL裁切先计算索引，通过索引获取裁切后的图像，同时重新计算geo_transformes里的左上角点坐标。核心代码示例如下：

xsize_t = in_band.XSize
ysize_t = in_band.YSize
# 计算偏移
y_offset = ysize_t - height
x_offset = xsize_t - width
data = in_band.ReadAsArray()

if pro_mode == "r_b": # 假设裁剪右下部分
    # 更新左上角点坐标 
    geotrans[3] = geotrans[3] + geotrans[5] * y_offset
    geotrans[0] = geotrans[0] + geotrans[2] * x_offset
    data = data[y_offset:, x_offset:]

完整代码示例：

from os import makedirs, remove
from os.path import exists, join
from osgeo import gdal
import numpy as np


def get_params(img_info, img_config):
    return img_info["indir"], img_info["outdir"], img_info["img_list"], img_info["img_resolution"], \
           img_info["pro_mode"], img_config["height"], img_config["width"], img_config["channel"], img_config["tar_resolution"]


def large2small(img_info, img_config):
    indir, outdir, img_list, img_resolutions, pro_modes, height, width, channel, tar_resolution = get_params(img_info, img_config)

    for f, img_resolution, pro_mode in zip(img_list, img_resolutions, pro_modes):
        fi = join(indir, f)
        fo_ = join(outdir, "_" + f)
        fo = join(outdir, f)

        # resize image
        in_ds = gdal.Open(fi, gdal.GA_ReadOnly)
        geotrans = in_ds.GetGeoTransform()
        geotrans = list(geotrans)
        geoproj = in_ds.GetProjection()
        geotrans[1] = tar_resolution  # 
        geotrans[5] = -tar_resolution  # 
        in_band = in_ds.GetRasterBand(1)
        xsize = in_band.XSize
        ysize = in_band.YSize
        ysize_t = int(ysize * img_resolution / tar_resolution)
        xsize_t = int(xsize * img_resolution / tar_resolution)
        if exists(fo_):
            remove(fo_)
        out_ds = in_ds.GetDriver().Create(fo_, xsize_t, ysize_t, 1,
                                          in_band.DataType)  # 创建一个构建重采样影像的句柄
        out_ds.SetProjection(in_ds.GetProjection())  # 设置投影信息
        geotrans = tuple(geotrans)
        out_ds.SetGeoTransform(geotrans)  # 设置地理变换信息
        data = np.empty((ysize_t, xsize_t), np.uint8)  # 设置一个与重采样影像行列号相等的矩阵去接受读取所得的像元值
        in_band.ReadAsArray(buf_obj=data)
        out_band = out_ds.GetRasterBand(1)
        out_band.WriteArray(data)
        out_band.FlushCache()
        out_band.ComputeStatistics(False)  # 计算统计信息
        # out_ds.BuildOverviews('average', [1, 2, 4, 8, 16, 32])  # 构建金字塔
        del in_ds  # 删除句柄
        del out_ds

        # reload and clip
        in_ds = gdal.Open(fo_, gdal.GA_ReadOnly)
        geotrans = in_ds.GetGeoTransform()
        geotrans = list(geotrans)
        geoproj = in_ds.GetProjection()
        in_band = in_ds.GetRasterBand(1)
        xsize_t = in_band.XSize
        ysize_t = in_band.YSize
        y_offset = ysize_t - height
        x_offset = xsize_t - width
        data = in_band.ReadAsArray()
        #
        # pro_mode = "l_u"
        if pro_mode == "l_b":
            geotrans[3] = geotrans[3] + geotrans[5] * y_offset
            data = data[y_offset:, :width]
        elif pro_mode == "l_u":
            data = data[:height, :width]
        elif pro_mode == "r_u":
            geotrans[0] = geotrans[0] + geotrans[2] * x_offset
            data = data[:height, x_offset:]
        elif pro_mode == "r_b":
            geotrans[3] = geotrans[3] + geotrans[5] * y_offset
            geotrans[0] = geotrans[0] + geotrans[2] * x_offset
            data = data[y_offset:, x_offset:]

        if exists(fo):
            remove(fo)
        options = ['COMPRESS=LZW']
        out_ds = in_ds.GetDriver().Create(fo, width, height, 1,
                                          in_band.DataType, options=options)  # 创建一个构建重采样影像的句柄
        out_ds.SetProjection(in_ds.GetProjection())  # 设置投影信息
        geotrans = tuple(geotrans)
        out_ds.SetGeoTransform(geotrans)  # 设置地理变换信息
        # data = np.empty((height, width), np.uint8)  # 设置一个与重采样影像行列号相等的矩阵去接受读取所得的像元值
        # in_band.ReadAsArray(buf_obj=data)

        out_band = out_ds.GetRasterBand(1)
        out_band.WriteArray(data)
        out_band.FlushCache()
        out_band.ComputeStatistics(False)  # 计算统计信息
        out_ds.BuildOverviews('average', [1, 2, 4, 8, 16, 32])  # 构建金字塔
        del in_ds  # 删除句柄
        del out_ds

        if exists(fo_):
            remove(fo_)