利用控制点实现图像重投影的专业指南

def simple_correction(n, m, set_points):
    """
    一个简单的线性缩放和位移函数，仅基于两个控制点。
    适用于没有旋转、倾斜或复杂非线性畸变的情况。
    """
    if len(set_points) != 2:
        raise ValueError("Simple correction requires exactly two control points.")

    p1_old, p1_new = (set_points[0]["old_x"], set_points[0]["old_y"]), (set_points[0]["new_x"], set_points[0]["new_y"])
    p2_old, p2_new = (set_points[1]["old_x"], set_points[1]["old_y"]), (set_points[1]["new_x"], set_points[1]["new_y"])

    # 计算X和Y方向的缩放因子
    scale_x = (p2_new[0] - p1_new[0]) / (p2_old[0] - p1_old[0]) if (p2_old[0] - p1_old[0]) != 0 else 1
    scale_y = (p2_new[1] - p1_new[1]) / (p2_old[1] - p1_old[1]) if (p2_old[1] - p1_old[1]) != 0 else 1

    # 计算X和Y方向的偏移量
    offset_x = p1_new[0] - p1_old[0] * scale_x
    offset_y = p1_new[1] - p1_old[1] * scale_y

    ans = []
    for j in range(m):
        row = []
        for i in range(n):
            new_x = offset_x + i * scale_x
            new_y = offset_y + j * scale_y
            row.append([new_x, new_y])
        ans.append(row)
    return ans

# 示例调用
# correction(5,5,[{"old_x":1,"old_y":1,"new_x":100,"new_y":150},{"old_x":3,"old_y":3,"new_x":200,"new_y":200}])

from osgeo import gdal, osr
import numpy as np
import os

def reproject_image_with_gcps(input_image_path, output_image_path, gcps_data, target_srs_wkt=None,
                              output_resolution_x=None, output_resolution_y=None,
                              resampling_algorithm=gdal.GRIORA_NearestNeighbour,
                              output_type=gdal.GDT_Int16,
                              use_tps=True):
    """
    使用GDAL通过控制点对图像进行重投影。

    Args:
        input_image_path (str): 输入图像文件的路径。
        output_image_path (str): 输出重投影图像文件的路径。
        gcps_data (list): 包含GCP数据的列表。每个元素应为字典，包含
                          'target_x', 'target_y', 'source_pixel', 'source_line'。
                          例如：[{'target_x': -111.931, 'target_y': 41.745, 'source_pixel': 1078, 'source_line': 648}, ...]
        target_srs_wkt (str, optional): 目标空间参考系统的WKT字符串。
                                        例如：'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,...]]'
                                        或通过 osr.SpatialReference().SetWellKnownGeogCS('WGS84').ExportToWkt() 生成。
                                        如果为None，则不设置目标SRS。
        output_resolution_x (float, optional): 输出图像的X方向分辨率。
        output_resolution_y (float, optional): 输出图像的Y方向分辨率。
        resampling_algorithm (gdal.GRIORA_*): 重采样算法，如 gdal.GRIORA_NearestNeighbour, gdal.GRIORA_Bilinear。
        output_type (gdal.GDT_*): 输出图像的数据类型，如 gdal.GDT_Int16, gdal.GDT_Byte。
        use_tps (bool): 是否使用薄板样条（TPS）变换。如果为False，则使用默认的GDAL变换（通常是多项式）。
    """
    # 1. 打开输入数据集
    dataset = gdal.Open(input_image_path, gdal.GA_ReadOnly)
    if dataset is None:
        print(f"无法打开图像文件: {input_image_path}")
        return

    # 2. 定义GCPs
    gcps = []
    for gcp_info in gcps_data:
        # Z坐标通常为0，除非进行三维变换
        gcps.append(gdal.GCP(gcp_info['target_x'], gcp_info['target_y'], 0,
                             gcp_info['source_pixel'], gcp_info['source_line']))

    # 3. 设置目标空间参考系统
    srs = None
    if target_srs_wkt:
        srs = osr.SpatialReference()
        srs.ImportFromWkt(target_srs_wkt)

    # 4. 应用GCPs并执行Warp
    # 创建一个内存中的数据集副本，用于设置GCPs，避免直接修改原始文件
    # 也可以直接对原始dataset设置GCPs，但如果原始文件是只读的，或者不想修改原始文件，则需要副本
    mem_ds = gdal.GetDriverByName('MEM').CreateCopy('', dataset)
    if srs:
        mem_ds.SetGCPs(gcps, srs.ExportToWkt())
    else:
        mem_ds.SetGCPs(gcps, '') # 如果没有目标SRS，则只设置GCPs，Warp会尝试推断或使用默认

    warp_options = {
        'format': 'GTiff',
        'resampleAlg': resampling_algorithm,
        'outputType': output_type,
        'dstNodata': 65535,  # 示例，根据实际数据类型调整
        'srcNodata': 65535,  # 示例，根据实际数据类型调整
    }

    if use_tps:
        warp_options['tps'] = True
    else:
        # 如果不使用TPS，可以指定多项式阶数，例如 order=1 (仿射), order=2, order=3
        # warp_options['order'] = 1
        pass # 默认行为通常是基于GCPs数量选择最佳多项式或通用变换

    if output_resolution_x and output_resolution_y:
        warp_options['xRes'] = output_resolution_x
        warp_options['yRes'] = output_resolution_y

    print(f"正在将图像重投影到: {output_image_path}...")
    dst_ds = gdal.Warp(output_image_path, mem_ds, **warp_options)

    if dst_ds is not None:
        print("图像重投影完成。")
        dst_ds = None # 释放数据集
    else:
        print("图像重投影失败。")

    dataset = None # 释放原始数据集
    mem_ds = None # 释放内存数据集

# --- 示例用法 ---
if __name__ == "__main__":
    # 创建一个虚拟的TIFF文件用于测试
    # 实际应用中，请替换为您的图像文件路径
    dummy_image_path = 'test_input.tiff'
    output_image_path = 'test_reprojected.tiff'

    # 生成一个简单的5x5的虚拟图像
    driver = gdal.GetDriverByName('GTiff')
    rows, cols = 50, 50
    dummy_ds = driver.Create(dummy_image_path, cols, rows, 1, gdal.GDT_Byte)
    dummy_band = dummy_ds.GetRasterBand(1)
    dummy_band.WriteArray(np.random.randint(0, 256, (rows, cols), dtype=np.uint8))
    dummy_ds = None # 关闭并保存虚拟图像

    # 定义控制点数据
    # 这里的GCPs数据是示例，实际中应根据图像和目标坐标系的实际情况来提供。
    # 目标坐标(-111.931075, 41.745836) 对应 图像像素(1078, 648)
    # 目标坐标(-111.901655, 41.749269) 对应 图像像素(531, 295)
    # 目标坐标(-111.899180, 41.739882) 对应 图像像素(722, 334)
    # 目标坐标(-111.930510, 41.728719) 对应 图像像素(102, 548)
    # 注意：示例中的GCPs是针对一个假设的更大图像的，此处为演示目的，将像素坐标调整到50x50范围内
    gcps_example = [
        {'target_x': -111.931075, 'target_y': 41.745836, 'source_pixel': 45, 'source_line': 5}, # 右上角附近
        {'target_x': -111.901655, 'target_y': 41.749269, 'source_pixel': 5, 'source_line': 5},   # 左上角附近
        {'target_x': -111.899180, 'target_y': 41.739882, 'source_pixel': 45, 'source_line': 45}, # 右下角附近
        {'target_x': -111.930510, 'target_y': 41.728719, 'source_pixel': 5, 'source_line': 45}    # 左下角附近
    ]

    # 设置目标空间参考系统为WGS84地理坐标系
    target_srs = osr.SpatialReference()
    target_srs.SetWellKnownGeogCS('WGS84')
    target_srs_wkt = target_srs.ExportToWkt()

    # 执行重投影
    reproject_image_with_gcps(
        input_image_path=dummy_image_path,
        output_image_path=output_image_path,
        gcps_data=gcps_example,
        target_srs_wkt=target_srs_wkt,
        output_resolution_x=0.0001, # 示例分辨率，根据目标坐标系和实际需求调整
        output_resolution_y=0.0001,
        resampling_algorithm=gdal.GRIORA_Bilinear, # 使用双线性插值
        output_type=gdal.GDT_Byte, # 虚拟图像是Byte类型
        use_tps=True # 使用TPS变换
    )

    # 清理虚拟文件
    if os.path.exists(dummy_image_path):
        os.remove(dummy_image_path)
    if os.path.exists(output_image_path):
        print(f"成功生成重投影图像: {output_image_path}")
        # os.remove(output_image_path) # 如果需要，可以取消注释以删除输出文件