使用QGIS Python API实现地理配准变换

嗨，这个问题我刚好折腾过好几次，用QGIS Python API直接调用这些地理配准变换真的可行，核心是用QGIS自带的QgsGeorefTransform类来处理不同的变换类型，再结合控制点和栅格操作就能搞定。我给你捋清楚具体步骤，再附上完整示例：

一、先搞懂核心类和变换类型对应关系 #

首先得知道你在Georeferencer插件里看到的那些变换，在API里是QgsGeorefTransform类的枚举值，对应关系如下：

• Linear → QgsGeorefTransform.Linear
• Polynomial 1 → QgsGeorefTransform.Polynomial1
• Polynomial 2 → QgsGeorefTransform.Polynomial2
• Helmert → QgsGeorefTransform.Helmert
• Projective → QgsGeorefTransform.Projective

二、具体实现步骤 #

1. 导入必要模块 #

先把需要的QGIS核心模块导入：

from qgis.core import (
    QgsRasterLayer, QgsGeorefTransform, QgsGeorefDataPoint,
    QgsPointXY, QgsRasterFileWriter, QgsCoordinateReferenceSystem, QgsGeorefUtils
)

2. 准备输入栅格和控制点 #

• 加载你要配准的栅格：

input_raster_path = "/path/to/your/input_raster.tif"
input_raster = QgsRasterLayer(input_raster_path, "Input Raster")
if not input_raster.isValid():
    print("输入栅格加载失败！请检查文件路径")

• 定义控制点：每个控制点需要原图的像素坐标和对应的目标地理坐标，注意不同变换需要的最少控制点数量：
  • Linear/Helmert：至少2个
  • Polynomial 1：至少3个
  • Polynomial 2：至少6个
  • Projective：至少4个
    示例控制点定义：

# 格式：(像素X, 像素Y) → (目标地理X, 目标地理Y)
control_points = [
    ((100, 200), (123456.789, 987654.321)),
    ((300, 400), (123556.789, 987554.321)),
    ((500, 600), (123656.789, 987454.321))
]

# 转换成QGIS识别的控制点对象列表
georef_points = []
for (pixel_xy, geo_xy) in control_points:
    pixel_point = QgsPointXY(pixel_xy[0], pixel_xy[1])
    geo_point = QgsPointXY(geo_xy[0], geo_xy[1])
    georef_points.append(QgsGeorefDataPoint(pixel_point, geo_point))

3. 初始化变换并计算参数 #

选择你需要的变换类型，实例化变换对象后传入控制点计算参数：

# 选择变换类型，这里以Polynomial 1为例，可按需替换
transform_type = QgsGeorefTransform.Polynomial1
transform = QgsGeorefTransform(transform_type)

# 传入控制点并计算变换参数
transform.setPoints(georef_points)
transform.updateParameters()

# 检查变换是否有效
if not transform.isValid():
    print("变换无效！请检查控制点数量是否达标或坐标是否正确")
else:
    print("变换参数计算成功！")
    # 可查看每个控制点的残差，验证配准精度
    for i, point in enumerate(georef_points):
        residual = transform.residual(point)
        print(f"控制点{i+1}的残差: {residual:.4f}")

4. 应用变换并输出配准后的栅格 #

用QGIS内置工具完成栅格的配准输出：

output_path = "/path/to/your/registered_raster.tif"
# 目标坐标系，按需修改为你需要的EPSG编码
target_crs = QgsCoordinateReferenceSystem("EPSG:4326")

# 执行栅格配准（重投影）
success = QgsGeorefUtils.warpRaster(
    input_raster,
    output_path,
    target_crs,
    transform,
    input_raster.width(),
    input_raster.height(),
    1,  # 重采样方法：0=最近邻，1=双线性，2=三次卷积，按需修改
    QgsRasterFileWriter.driverForExtension("tif")
)

if success:
    print(f"配准完成！输出文件路径：{output_path}")
else:
    print("配准失败，请检查参数设置")

三、一些实用小提醒 #

• 控制点尽量选均匀分布、特征明显的位置（比如道路交叉口、建筑角点），能大幅提升配准精度
• 重采样方法按需选择：遥感影像常用双线性/三次卷积，分类图或矢量转栅格用最近邻能保留原始值
• 如果需要批量处理，可以把上述代码封装成函数，循环读取控制点文件和栅格文件即可

备注：内容来源于stack exchange，提问作者Shubham_geo