这个问题其实是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学习笔记