如何将shapely.Polygon叠加到带GCP的倾斜TIFF影像上？

1. 确认GCP是否被正确应用 #

倾斜影像的地理配准依赖GCP，但部分TIFF文件仅嵌入GCP未生成带旋转的有效变换矩阵。用rasterio检查并手动生成正确变换：

import rasterio
from rasterio.transform import from_gcps

with rasterio.open("your_tilted_image.tif") as src:
    # 查看文件中存储的GCP
    print("Embedded GCPs:", src.gcps[0])
    # 对比默认变换与GCP生成的变换
    gcp_transform = from_gcps(src.gcps[0]) if src.gcps[0] else None
    print("Default transform:", src.transform)
    print("GCP-derived transform:", gcp_transform)

若两个变换差异显著，后续处理必须使用gcp_transform而非默认的src.transform——这是旋转参数未被正确处理的核心原因。

2. 精准转换多边形坐标系统 #

将EPSG:4326的Shapely多边形转换至影像的EPSG:32617坐标系，必须使用专业投影工具避免顺序错误：

from shapely.ops import transform
from pyproj import Transformer

# 定义坐标转换器，强制xy顺序（符合GIS与Shapely规范）
proj_transformer = Transformer.from_crs("EPSG:4326", "EPSG:32617", always_xy=True)
# 执行多边形投影转换
projected_polygon = transform(proj_transformer.transform, your_shapely_polygon)

3. 基于带旋转的变换矩阵可视化叠加 #

倾斜影像的extent需结合旋转参数计算，不能直接使用轴对齐的bounds：

import matplotlib.pyplot as plt

with rasterio.open("your_tilted_image.tif") as src:
    img_data = src.read(1)
    # 优先使用GCP生成的变换矩阵
    current_transform = from_gcps(src.gcps[0]) if src.gcps[0] else src.transform
    
    # 计算包含旋转的影像地理范围
    left = current_transform.c
    right = current_transform.c + current_transform.a * src.width + current_transform.b * src.height
    bottom = current_transform.f + current_transform.d * src.width + current_transform.e * src.height
    top = current_transform.f

    # 绘制影像与多边形
    plt.figure(figsize=(12,12))
    plt.imshow(img_data, extent=(left, right, bottom, top), cmap="gray")
    poly_x, poly_y = projected_polygon.exterior.xy
    plt.plot(poly_x, poly_y, color="#ff3333", linewidth=2)
    plt.axis("equal")
    plt.show()

4. 排除数据质量问题 #

若上述操作后仍存在大偏差：

• 验证GCP的坐标是否确实为EPSG:32617的UTM坐标，排查是否误将经纬度作为UTM输入
• 手动选取影像中清晰特征点（如墙角），对比其GCP坐标与多边形对应点的投影坐标，确认偏差来源是GCP精度还是处理流程

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