使用Rasterio from_gcps重投影仅调边界未移像素问题求助

我的栅格与Google Earth地面真值错位，尝试用rasterio.warp.reproject结合手动获取的GCPs做对齐，简化代码如下：

import rasterio
import rasterio.crs
from rasterio.control import GroundControlPoint
from rasterio.transform import from_gcps
from rasterio.warp import reproject, Resampling
from rasterio.crs import CRS

if __name__ == '__main__':
    # Open the input TIFF file
    with rasterio.open(input_img_path) as src:
        gcps = [
            GroundControlPoint(col=(...),row=(...),x=(...),y=(...)),
            (...)
        ] # x and y are lon and lat extracted from Google Earth
        transform = from_gcps(gcps)
        gcp_crs = CRS.from_epsg(4326)  # This is the crs of the GCPs
        # Warp the image
        kwargs = src.profile.copy()
        kwargs.update(transform=transform, crs=gcp_crs)
        with rasterio.open(output_img_path, "w", **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    dst_transform=dst.transform,
                    dst_crs=dst.crs,
                    resampling=Resampling.nearest
                )

运行后仅栅格边界正确对齐，像素位置完全没变化，QGIS中可见：

• 调整前：栅格与地面真值错位
• 重投影后：栅格边界（黑色方框）位置正确
• 原图与重投影图对比：仅边界变化，像素未移动

试过直接在reproject传GCPs、保留源CRS、更换GCP组合等方法，都没解决，求方案。

解决方案 #

问题核心是直接用GCP生成的transform作为输出transform，但未让重投影过程基于GCP做真正的几何校正。正确步骤如下：

修正后的代码 #

import rasterio
from rasterio.control import GroundControlPoint
from rasterio.warp import reproject, Resampling, calculate_default_transform
from rasterio.crs import CRS

if __name__ == '__main__':
    input_img_path = "你的输入栅格路径.tif"
    output_img_path = "校正后的输出路径.tif"
    
    # 手动获取的GCPs：col/row为原图0-based像素坐标，x/y为WGS84经纬度
    gcps = [
        GroundControlPoint(col=100, row=200, x=116.38, y=39.90),
        GroundControlPoint(col=500, row=200, x=116.39, y=39.90),
        GroundControlPoint(col=100, row=600, x=116.38, y=39.89),
        GroundControlPoint(col=500, row=600, x=116.39, y=39.89),
        # 推荐4个及以上GCP以提升校正精度
    ]
    
    with rasterio.open(input_img_path) as src:
        # 定义GCP对应的CRS（WGS84）
        gcp_crs = CRS.from_epsg(4326)
        
        # 计算校正后输出的transform、宽高
        # 目标CRS可按需替换（如UTM投影），这里用WGS84
        dst_crs = CRS.from_epsg(4326)
        dst_transform, dst_width, dst_height = calculate_default_transform(
            gcp_crs, dst_crs, src.width, src.height,
            # 若需固定输出范围，可手动指定left/bottom/right/top
            # 保留原分辨率可添加参数：dst_resolution=(src.res[0], src.res[1])
        )
        
        # 更新输出文件配置
        kwargs = src.profile.copy()
        kwargs.update({
            'crs': dst_crs,
            'transform': dst_transform,
            'width': dst_width,
            'height': dst_height,
            'nodata': src.nodata
        })
        
        with rasterio.open(output_img_path, 'w', **kwargs) as dst:
            for i in range(1, src.count + 1):
                reproject(
                    source=rasterio.band(src, i),
                    destination=rasterio.band(dst, i),
                    src_transform=src.transform,
                    src_crs=src.crs,
                    # 关键：传入GCP列表，让rasterio基于GCP做几何校正
                    src_gcps=[(gcps, gcp_crs)],
                    dst_transform=dst_transform,
                    dst_crs=dst_crs,
                    resampling=Resampling.nearest
                )

关键要点 #

• 必须传入src_gcps参数：这是让rasterio启用GCP几何校正的核心，而非单纯执行投影转换。
• 避免直接用from_gcps(gcps)作为输出transform：该transform仅为GCP拟合结果，未结合原图像像素网格重采样，会导致输出像素与原图一一对应、仅边界坐标变更。
• 保证GCP质量：
  • 至少3个GCP可做 affine 变换，4个及以上支持 projective 变换，精度更高。
  • 确认col/row为原图0-based像素坐标，x/y为正确的WGS84经纬度。
• 按需调整输出分辨率/范围：通过calculate_default_transform的参数可手动控制输出图像的分辨率和边界，适配不同需求。

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