如何在python中使用地面控制点对未引用的航空影像进行地理配准

import gdal
import osr

image_path = '未引用的航空影像路径'
control_points_path = '地面控制点文件路径'

image_dataset = gdal.Open(image_path, gdal.GA_Update)
control_points_dataset = gdal.Open(control_points_path)

control_points_layer = control_points_dataset.GetLayer()
control_points = []

for feature in control_points_layer:
    geometry = feature.GetGeometryRef()
    x = geometry.GetX()
    y = geometry.GetY()
    image_x = feature.GetField('影像X坐标')
    image_y = feature.GetField('影像Y坐标')
    control_points.append((x, y, image_x, image_y))

image_projection = image_dataset.GetProjection()
image_geotransform = image_dataset.GetGeoTransform()

target_projection = osr.SpatialReference()
target_projection.ImportFromEPSG(4326)  # 设置目标投影坐标系，这里以WGS84为例

transform = osr.CoordinateTransformation(target_projection, image_projection)

for point in control_points:
    x, y, image_x, image_y = point
    image_x_geo, image_y_geo, _ = transform.TransformPoint(x, y)

    # 更新影像的地理坐标
    image_geotransform[0] = image_x_geo
    image_geotransform[3] = image_y_geo
    image_dataset.SetGeoTransform(image_geotransform)

    # 更新影像的像素坐标
    image_x_pixel = int((image_x_geo - image_geotransform[0]) / image_geotransform[1])
    image_y_pixel = int((image_y_geo - image_geotransform[3]) / image_geotransform[5])
    image_dataset.GetRasterBand(1).WriteArray([image_x_pixel], image_x, image_y)
    image_dataset.GetRasterBand(2).WriteArray([image_y_pixel], image_x, image_y)

image_dataset.FlushCache()

image_dataset = None
control_points_dataset = None