获取两个EPSG:32637坐标系TIFF文件的交集并裁剪保存方法

核心思路 #

先获取两个TIFF的地理边界，计算出它们的交集范围，再使用GDAL工具（命令行或Python API）将两个TIFF统一裁剪到该交集区域，确保输出文件的坐标和尺寸完全一致。

方案一：GDAL命令行实现（快速便捷） #

步骤1：获取两个TIFF的地理边界 #

使用gdalinfo命令提取每个文件的左上角（Upper Left）和右下角（Lower Right）坐标：

# 获取第一个TIFF的边界
gdalinfo tiff1.tif | grep "Upper Left\|Lower Right"
# 获取第二个TIFF的边界
gdalinfo tiff2.tif | grep "Upper Left\|Lower Right"

记录结果中的坐标值，例如：

• tiff1: Upper Left (x1_ul, y1_ul)，Lower Right (x1_lr, y1_lr)
• tiff2: Upper Left (x2_ul, y2_ul)，Lower Right (x2_lr, y2_lr)

步骤2：计算交集区域的边界 #

EPSG:32637属于UTM北带，y值越大越靠北，因此交集范围的计算规则为：

• 交集左上角x = max(x1_ul, x2_ul)
• 交集左上角y = min(y1_ul, y2_ul)
• 交集右下角x = min(x1_lr, x2_lr)
• 交集右下角y = max(y1_lr, y2_lr)

注意：如果计算后出现交集左上角x >= 右下角x或左上角y <= 右下角y，说明两个TIFF无交集，无法完成裁剪。

步骤3：统一裁剪两个TIFF #

使用gdalwarp命令，指定交集范围（-te参数顺序为xmin ymin xmax ymax），并保证分辨率一致：

# 裁剪第一个TIFF，若原分辨率不同，添加-tr 参数指定统一分辨率（例如-tr 0.5 0.5）
gdalwarp -te [交集xmin] [交集ymin] [交集xmax] [交集ymax] \
         -t_srs EPSG:32637 -co COMPRESS=LZW \
         tiff1.tif tiff1_cropped.tif

# 用完全相同的参数裁剪第二个TIFF
gdalwarp -te [交集xmin] [交集ymin] [交集xmax] [交集ymax] \
         -t_srs EPSG:32637 -co COMPRESS=LZW \
         tiff2.tif tiff2_cropped.tif

方案二：Python GDAL API实现（自动化批量处理） #

适合需要批量处理或集成到工作流的场景，脚本如下：

from osgeo import gdal

def get_raster_extent(raster_path):
    """获取栅格的地理范围（xmin, ymin, xmax, ymax）"""
    ds = gdal.Open(raster_path)
    if not ds:
        raise FileNotFoundError(f"无法打开文件: {raster_path}")
    geotrans = ds.GetGeoTransform()
    x_min = geotrans[0]
    y_max = geotrans[3]
    x_max = x_min + geotrans[1] * ds.RasterXSize
    y_min = y_max + geotrans[5] * ds.RasterYSize
    ds = None
    return (x_min, y_min, x_max, y_max)

def crop_raster(input_path, output_path, bounds, res=None):
    """裁剪栅格到指定范围，统一分辨率"""
    warp_opts = gdal.WarpOptions(
        outputBounds=bounds,
        dstSRS="EPSG:32637",
        creationOptions=["COMPRESS=LZW"],
        xRes=res[0] if res else None,
        yRes=res[1] if res else None
    )
    gdal.Warp(output_path, input_path, options=warp_opts)

# 主流程
if __name__ == "__main__":
    # 输入文件路径
    tiff1_path = "tiff1.tif"
    tiff2_path = "tiff2.tif"
    
    # 获取两个栅格的范围
    ext1 = get_raster_extent(tiff1_path)
    ext2 = get_raster_extent(tiff2_path)
    
    # 计算交集范围
    intersect_bounds = (
        max(ext1[0], ext2[0]),
        max(ext1[1], ext2[1]),
        min(ext1[2], ext2[2]),
        min(ext1[3], ext2[3])
    )
    
    # 检查交集有效性
    x_min, y_min, x_max, y_max = intersect_bounds
    if x_min >= x_max or y_min >= y_max:
        raise ValueError("两个TIFF文件不存在交集区域")
    
    # 统一分辨率（采用第一个栅格的分辨率）
    ds = gdal.Open(tiff1_path)
    geotrans = ds.GetGeoTransform()
    target_res = (geotrans[1], abs(geotrans[5]))
    ds = None
    
    # 执行裁剪
    crop_raster(tiff1_path, "tiff1_cropped.tif", intersect_bounds, target_res)
    crop_raster(tiff2_path, "tiff2_cropped.tif", intersect_bounds, target_res)

内容的提问来源于stack exchange，提问作者Muhammad Muneeb Ur Rahman