GDAL给定一套格网重投影坐标进行影像投影转换

要使用GDAL进行影像投影转换，需要进行以下步骤：

1. 打开原始影像和目标影像的数据集，例如：

from osgeo import gdal

src_ds = gdal.Open('input.tif')
dst_ds = gdal.Open('output.tif', gdal.GA_Update)

2. 获取原始影像和目标影像的投影信息和地理转换信息，例如：

src_proj = src_ds.GetProjection()
dst_proj = dst_ds.GetProjection()

src_geotrans = src_ds.GetGeoTransform()
dst_geotrans = dst_ds.GetGeoTransform()

3. 创建一个投影转换器，并设置投影转换的参数，例如：

from osgeo import osr

src_srs = osr.SpatialReference()
src_srs.ImportFromWkt(src_proj)

dst_srs = osr.SpatialReference()
dst_srs.ImportFromWkt(dst_proj)

transform = osr.CoordinateTransformation(src_srs, dst_srs)

4. 读取原始影像的像素数据，并进行投影转换，例如：

import numpy as np

band = src_ds.GetRasterBand(1)
data = band.ReadAsArray()

cols, rows = data.shape

x = np.arange(0, cols) * src_geotrans[1] + src_geotrans[0]
y = np.arange(0, rows) * src_geotrans[5] + src_geotrans[3]

x, y = np.meshgrid(x, y)

x_new, y_new, _ = transform.TransformPoints(np.vstack((x.flatten(), y.flatten())).T).T
x_new = x_new.reshape((rows, cols))
y_new = y_new.reshape((rows, cols))

new_geotrans = (x_new.min(), src_geotrans[1], 0, y_new.max(), 0, -src_geotrans[5])

data_new = np.zeros((rows, cols), dtype=data.dtype)

gdal.ReprojectImage(src_ds, dst_ds, src_proj, dst_proj, gdal.GRA_Bilinear, 0, 0, transform)

5. 保存目标影像的像素数据，并关闭数据集，例如：

band = dst_ds.GetRasterBand(1)
band.WriteArray(data_new)

# 设置目标影像的地理转换信息
dst_ds.SetGeoTransform(new_geotrans)
dst_ds.SetProjection(dst_srs.ExportToWkt())

# 关闭数据集
dst_ds = None
src_ds = None