GDAL创建栅格写入坐标系无效的解决办法

这个问题其实是CGCS2000坐标系数据处理的通病了，是PROJ库的坐标系数据库版本太旧不包含CSCG2000坐标系造成的。

---

今天写脚本时遇到一个bug一直解决不了，我只是想实现一个简单的图片重采样功能，将原栅格缩放至指定大小并保留坐标信息（CGCS2000的投影坐标系），但遇到的问题是生成的图像用arcgis打开后显示未知坐标系。脚本代码如下（GPT生成的，我一开始还以为是代码有问题）

from osgeo import gdal,osr
def resample_image(input_path, output_path, new_width, new_height):
    # 打开输入图像
    input_ds = gdal.Open(input_path)
    if input_ds is None:
        print("无法打开输入图像。")
        return

    # 获取输入图像的原始大小和投影信息
    input_width = input_ds.RasterXSize
    input_height = input_ds.RasterYSize
    input_projection = input_ds.GetProjection()
    input_geo_transform = input_ds.GetGeoTransform()

    # 创建输出图像
    driver = gdal.GetDriverByName('GTiff')
    output_ds = driver.Create(output_path, new_width, new_height, input_ds.RasterCount, input_ds.GetRasterBand(1).DataType)
    output_ds.SetProjection(input_projection)

    # 计算新的地理变换信息
    x_res = input_geo_transform[1] * (input_width / new_width)
    y_res = input_geo_transform[5] * (input_height / new_height)
    output_geo_transform = (input_geo_transform[0], x_res, 0, input_geo_transform[3], 0, y_res)

    output_ds.SetGeoTransform(output_geo_transform)

    # 执行重采样
    resampling_method = gdal.GRA_Bilinear  # 使用双线性插值
    gdal.ReprojectImage(input_ds, output_ds, input_projection, input_projection, resampling_method)

    # 关闭数据集
    input_ds = None
    output_ds = None

    print("重采样完成。")

坐标信息，无非是仿射矩阵和投影参数，代码逻辑很简单，所以不是代码问题。我尝试许久之后用c++写了一遍发现问题解决了，我想起了C++的代码中重设置了PROJ库的环境变量，而Python没有，那么真相就明确了，由于环境变量中有POSTGIS数据库的PROJ库，于是脚本使用了那个版本相对较老的PROJ.DB坐标系数据库，而那个库里没有新的坐标系信息。

加上环境变量设置后，问题解决

import os
os.environ['PROJ_LIB'] = r'D:\coding\gdal_compile\GDAL_353_DEBUG\bin'

要不是我之前写c++和GDAL时遇到了proj库冲突的warning，我估计想到这个bug的来源问题会花更长的时间。

另附个人博客：www.gislxz.com 内有GDAL学习笔记