Python中使用GDAL从GeoTIFF栅格提取旋转矩形区域的方法

你之前使用的ReadAsArray(xoffset, yoffset, xsize, ysize)仅支持读取和栅格坐标轴平行的矩形区域，无法直接读取旋转矩形。可通过内存虚拟数据集+Warp重采样的方式实现旋转区域提取，最终直接输出可直接使用的numpy数组，无需生成中间临时文件，输出数组格式和原方法读取结果完全兼容。

所需依赖 #

确保Python环境已安装以下依赖：

• osgeo.gdal（GDAL的Python绑定）
• numpy

实现代码 #

from osgeo import gdal
import numpy as np

def extract_rotated_rect(
    raster_path: str,
    rect_corners: np.ndarray,
    out_pixel_width: int,
    out_pixel_height: int,
    resample_method: int = gdal.GRA_Bilinear
) -> np.ndarray:
    """
    从GeoTIFF中提取旋转矩形区域，返回轴对齐的numpy数组
    参数说明：
        raster_path: 输入GeoTIFF文件路径
        rect_corners: 旋转矩形4个角点的地理坐标，shape为(4,2)，顺序按[左上,右上,右下,左下]排列
        out_pixel_width: 输出结果的像素宽度（列数）
        out_pixel_height: 输出结果的像素高度（行数）
        resample_method: 重采样方法，默认双线性，分类数据建议替换为gdal.GRA_NearestNeighbour
    """
    # 打开源栅格
    src_ds = gdal.Open(raster_path)
    src_srs = src_ds.GetSpatialRef()
    nodata_val = src_ds.GetRasterBand(1).GetNoDataValue()

    # 计算旋转矩形在地理坐标系下的实际宽高，推导输出分辨率
    geo_width = np.linalg.norm(rect_corners[1] - rect_corners[0])
    geo_height = np.linalg.norm(rect_corners[2] - rect_corners[1])
    pixel_size_x = geo_width / out_pixel_width
    pixel_size_y = geo_height / out_pixel_height

    # 定义输出栅格的仿射变换，让旋转矩形刚好轴对齐铺满整个输出数组
    out_geotransform = (
        rect_corners[0][0],
        pixel_size_x,
        0,
        rect_corners[0][1],
        0,
        -pixel_size_y  # 栅格y轴向下，分辨率取负值
    )

    # 创建内存中的目标栅格，不写入本地磁盘
    mem_driver = gdal.GetDriverByName('MEM')
    dst_ds = mem_driver.Create(
        '',
        out_pixel_width,
        out_pixel_height,
        src_ds.RasterCount,
        src_ds.GetRasterBand(1).DataType
    )
    dst_ds.SetGeoTransform(out_geotransform)
    dst_ds.SetSpatialRef(src_srs)

    # 构造旋转矩形的闭合多边形作为裁剪边界
    closed_corners = np.vstack([rect_corners, rect_corners[0]])
    cutline_wkt = "POLYGON (({}))".format(
        ",".join([f"{x} {y}" for x,y in closed_corners])
    )

    # 执行裁剪、重采样、旋转变换，结果直接写入内存数据集
    gdal.Warp(
        dst_ds,
        src_ds,
        cutlineWKT=cutline_wkt,
        cutlineSR=src_srs,
        resampleAlg=resample_method,
        dstNodata=nodata_val
    )

    # 读取为标准numpy数组
    out_arr = dst_ds.ReadAsArray()

    # 释放数据集资源
    src_ds = None
    dst_ds = None

    return out_arr


# 调用示例
if __name__ == "__main__":
    # 替换为目标旋转矩形的4个角点地理坐标，顺序为左上、右上、右下、左下
    target_corners = np.array([
        [120.10, 30.20],
        [120.15, 30.22],
        [120.13, 30.17],
        [120.08, 30.15]
    ])
    # 提取结果，输出为500*300像素的numpy数组
    data = extract_rotated_rect(
        raster_path="path/to/raster.tiff",
        rect_corners=target_corners,
        out_pixel_width=500,
        out_pixel_height=300,
        resample_method=gdal.GRA_NearestNeighbour # 分类数据（如土地覆盖）用最邻近插值
    )
    # 可直接对data做numpy计算
    print(data.shape, data.dtype)

注意事项 #

• 如果没有提前计算4个角点坐标，仅知道矩形中心点地理坐标、矩形地理宽高、旋转角度，可先通过平面几何三角函数算出4个角点再传入函数。
• 重采样方法按需选择：连续值数据（如DEM、气温、遥感反射率）用双线性/三次卷积，分类数据必须用最邻近插值，避免类别值被篡改。
• 提取出的数组中，旋转矩形范围外的区域会自动填充源数据的NoData值，和原ReadAsArray返回的数组使用逻辑完全一致。

内容的提问来源于stack exchange，提问作者tobycoleman