如何利用另一影像的GCP对偏移GeoTIFF影像进行地理配准

第一步：修正现有代码的错误 #

你提供的代码存在变量名调用、未定义等问题，先修正如下：

import cv2
from PIL import Image
import numpy as np
from osgeo import gdal
from affine import Affine  # 需提前安装：pip install affine

# 读取影像
reference_image = Image.open(r'reference_image.tif')
target_image = Image.open(r'target_image.tif')
reference = np.array(reference_image)
target = np.array(target_image)

# SIFT特征匹配
sift = cv2.SIFT_create()
keypoints1, descriptors1 = sift.detectAndCompute(reference, None)
keypoints2, descriptors2 = sift.detectAndCompute(target, None)

bf = cv2.BFMatcher()
matches = bf.match(descriptors1, descriptors2)
matches = sorted(matches, key=lambda x: x.distance)

# 获取前7个最优匹配点对
points1 = [(keypoints1[match.queryIdx].pt[0], keypoints1[match.queryIdx].pt[1]) for match in matches[:7]]
points2 = [(keypoints2[match.trainIdx].pt[0], keypoints2[match.trainIdx].pt[1]) for match in matches[:7]]

# 读取地理变换参数
reference_tif = gdal.Open(r'reference_image.tif')
target_tif = gdal.Open(r'target_image.tif')
reference_gt = reference_tif.GetGeoTransform()
target_gt = target_tif.GetGeoTransform()  # 修正：调用GDAL影像对象的地理变换方法

# 将影像坐标转换为地理坐标
reference_coords = []
target_coords = []

for pt in points1:
    affine_transform = Affine.from_gdal(*reference_gt)
    x, y = affine_transform * (pt[0], pt[1])
    reference_coords.append((x, y))

for pt in points2:
    affine_transform = Affine.from_gdal(*target_gt)
    x, y = affine_transform * (pt[0], pt[1])
    target_coords.append((x, y))

第二步：计算空间变换模型并对齐影像 #

下面提供两种主流方法，分别适用于GIS分析场景和快速视觉验证场景：

方法1：GDAL驱动的GCP重采样（保留地理参考）

该方法会完整保留影像的地理坐标信息，适合后续GIS分析流程：

from osgeo import osr

# 1. 创建GCP配对：将目标影像的像素坐标与参考影像对应的地理坐标绑定
gcps = []
for idx in range(len(points2)):
    # GCP参数：目标地理X, 目标地理Y, 高程(设为0), 目标影像像素X, 目标影像像素Y
    gcp = gdal.GCP(reference_coords[idx][0], reference_coords[idx][1], 0, 
                   points2[idx][0], points2[idx][1])
    gcps.append(gcp)

# 2. 为目标影像设置GCP和参考投影
target_tif.SetGCPs(gcps, reference_tif.GetProjection())

# 3. 创建输出影像（匹配参考影像的尺寸、地理范围）
output_path = r'aligned_target.tif'
driver = gdal.GetDriverByName('GTiff')
output_ds = driver.Create(output_path, 
                          reference_tif.RasterXSize, reference_tif.RasterYSize,
                          target_tif.RasterCount, target_tif.GetRasterBand(1).DataType)
output_ds.SetGeoTransform(reference_gt)
output_ds.SetProjection(reference_tif.GetProjection())

# 4. 重采样对齐（双线性插值适合连续型影像）
gdal.ReprojectImage(target_tif, output_ds, 
                    target_tif.GetProjection(), reference_tif.GetProjection(),
                    gdal.GRA_Bilinear)

# 释放资源
output_ds = None
target_tif = None
reference_tif = None

方法2：OpenCV扭曲变换（快速视觉验证）

如果仅需要快速得到对齐后的视觉结果，可直接用OpenCV计算变换矩阵并应用扭曲：

# 将匹配点转换为OpenCV所需的数组格式
src_pts = np.float32(points2).reshape(-1,1,2)
dst_pts = np.float32(points1).reshape(-1,1,2)

# 计算透视变换矩阵（RANSAC过滤异常匹配点）
M, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)

# 应用扭曲变换，输出与参考影像同尺寸的对齐结果
aligned_image = cv2.warpPerspective(target, M, (reference.shape[1], reference.shape[0]),
                                    flags=cv2.INTER_LINEAR)

# 保存预览图（无地理参考，仅用于视觉验证）
cv2.imwrite(r'aligned_preview.png', aligned_image)

关键注意事项 #

• 若影像仅存在平移、旋转、缩放类的仿射偏移，可替换为cv2.getAffineTransform计算变换矩阵，效率更高
• GDAL方法是GIS生产场景的首选，能保证地理信息的完整性；OpenCV方法仅适合快速验证对齐效果
• 匹配点数量建议不少于4个（透视变换的最低要求），偏移越复杂，需要的匹配点数量越多

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