Python导出的GeoTIFF在QGIS中读取时的投影异常问题

问题描述 #

尝试用Python导出Landsat的RGB GeoTIFF文件（已重投影至WGS84），使用的numpy数组包括rgb(H,W,3)、lons(H,W)、lats(H,W)。代码运行正常，用gdalinfo查看输出文件时坐标系统、范围等信息均正确（约东经15°/北纬4°），但导入QGIS后空间范围完全错误，显示为数据的整数维度而非经纬度范围。

用户提供的原始代码：

import numpy as np
from osgeo import gdal, osr

driver = gdal.GetDriverByName('GTiff')
options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF']

# Define the data extent
extent = [np.min(lons), np.min(lats), np.max(lons), np.max(lats)]

# Get dimensions and set geotransform
nx = rgb.shape[0]
ny = rgb.shape[1]
data_type = gdal.GDT_Byte
xres = (extent[2] - extent[0]) / float(nx)
yres = (extent[3] - extent[1]) / float(ny)
geotransform = (extent[0], xres, 0, extent[3], 0, -yres)

# Create a temp grid
grid_data = driver.Create('export_rgb.tif', ny, nx, 3, data_type, options=options)

# Setup projection and geo-transform
# Lat/Lon WSG84 Spatial Reference System
grid_data.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
grid_data.SetProjection(srs.ExportToWkt())

# Write data for each band
grid_data.GetRasterBand(1).WriteArray(rgb[:, :, 0])
grid_data.GetRasterBand(2).WriteArray(rgb[:, :, 1])
grid_data.GetRasterBand(3).WriteArray(rgb[:, :, 2])

# Save the file
grid_data.FlushCache()

# Close the file
driver = None
grid_data = None

问题根源 #

核心问题是GDAL Create方法的宽高参数顺序错误：
GDAL的Create方法要求参数顺序为(filename, xsize, ysize, bands, dtype, options)，其中xsize是图像的列数（对应numpy数组的W，即rgb.shape[1]），ysize是图像的行数（对应numpy数组的H，即rgb.shape[0]）。你当前代码中把ny（列数）和nx（行数）的顺序搞反了，导致GDAL对图像维度的识别与实际数据不匹配，最终QGIS解析空间范围时出现错误。

同时，分辨率计算时也需要对应正确的宽高值，否则会导致空间分辨率与实际不符。

修正后的代码 #

import numpy as np
from osgeo import gdal, osr

driver = gdal.GetDriverByName('GTiff')
options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF']

# Define the data extent
extent = [np.min(lons), np.min(lats), np.max(lons), np.max(lats)]

# 正确获取维度：ysize是行数(H)，xsize是列数(W)
ysize, xsize = rgb.shape[:2]
data_type = gdal.GDT_Byte
# 用xsize计算x方向分辨率，ysize计算y方向分辨率
xres = (extent[2] - extent[0]) / float(xsize)
yres = (extent[3] - extent[1]) / float(ysize)
geotransform = (extent[0], xres, 0, extent[3], 0, -yres)

# 修正Create的xsize和ysize顺序
grid_data = driver.Create('export_rgb.tif', xsize, ysize, 3, data_type, options=options)

# 设置投影和地理变换
grid_data.SetGeoTransform(geotransform)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
grid_data.SetProjection(srs.ExportToWkt())

# 写入各波段数据
grid_data.GetRasterBand(1).WriteArray(rgb[:, :, 0])
grid_data.GetRasterBand(2).WriteArray(rgb[:, :, 1])
grid_data.GetRasterBand(3).WriteArray(rgb[:, :, 2])

# 保存并关闭文件
grid_data.FlushCache()
grid_data = None
driver = None

额外说明 #

如果你的lons和lats是不规则网格（即每个像素的经纬度不是均匀间隔的），上述规则网格的地理变换方式仍会存在偏差，此时需要先将数据重采样到规则格网，或者使用GDAL的地面控制点（GCP）来定义空间映射。但从你描述gdalinfo显示正确的情况来看，只需修正宽高参数顺序即可解决QGIS中的空间范围问题。

内容的提问来源于stack exchange，提问作者AlexT60