使用GDAL地理配准影像时输出空白TIFF文件问题咨询

你遇到的空白问题由至少2个明确错误导致，geotransform的经纬度对应逻辑你写对了，但存在计算错误和数据类型匹配问题：

• 经度对应geotransform里的x轴（水平方向），纬度对应y轴（垂直方向），你写的(minLon, resLon, 0, maxLat, 0, -resLat)参数顺序逻辑是正确的
• 第一个硬错误：计算纬度方向分辨率时分母错用了x方向的像素数nLons，导致计算出的y方向分辨率完全错误，栅格地理范围和实际像素尺寸不匹配，GIS软件加载时无法定位到正确的像素位置
• 第二个高频错误：你固定使用gdal.GDT_Int32作为输出数据类型，如果你的array存储的是0-1区间浮点值、或者小数值浮点观测值，写入时强转为32位整型会全部截断为0，视觉上就是全空白影像
• 另外你没有手动释放GDAL数据集对象，部分环境下会导致文件写入不完整，也会出现空白问题

修复后可直接运行的代码 #

#!/usr/bin/env python3
import numpy as np
from osgeo import gdal
from osgeo import osr

gdal.UseExceptions() # 开启异常提示，方便定位报错

# 加载形状为(197, 250, 3)的数组，第三维度依次为像素值、经度、纬度
data = np.load("data.npy")
array = data[:, :, 0]
lon = data[:, :, 1]
lat = data[:, :, 2]

nLons = array.shape[1]  # x方向（经度）列数
nLats = array.shape[0]  # y方向（纬度）行数

# 计算经纬度范围
maxLon, minLon = lon.max(), lon.min()
maxLat, minLat = lat.max(), lat.min()

# 修正分辨率计算：两个方向分辨率分别除以对应方向的像素数
resLon = (maxLon - minLon) / nLons
resLat = (maxLat - minLat) / nLats
geotransform = (minLon, resLon, 0, maxLat, 0, -resLat)

# 自动匹配GDAL数据类型，避免强转丢值
if np.issubdtype(array.dtype, np.floating):
    gdal_dtype = gdal.GDT_Float32
else:
    gdal_dtype = gdal.GDT_Int32

# 创建输出栅格
driver = gdal.GetDriverByName('GTiff')
output_raster = driver.Create('myRaster_fixed.tif', nLons, nLats, 1, gdal_dtype)
output_raster.SetGeoTransform(geotransform)

# 设置WGS84坐标系
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
output_raster.SetProjection(srs.ExportToWkt())

# 写入数组
output_band = output_raster.GetRasterBand(1)
output_band.WriteArray(array)
# 如有需要可设置无效值
# output_band.SetNoDataValue(-9999)
output_raster.FlushCache()

# 手动释放资源，确保文件写入完整
del output_raster, output_band

后续排查提示 #

• 如果运行后加载影像还是显示空白，先查看栅格的像素统计值：如果统计值存在非0值，只是软件显示为白色，是GIS软件的渲染拉伸问题，手动调整像素值拉伸范围即可
• 如果你的原始经纬度不是等间隔规则网格，不能直接用仿射变换参数做配准，需要逐像素生成GCP控制点，再调用gdal.Warp做校正
• 写入前先打印array.min()、array.max()确认数组本身不是全0值，排除数据源本身的问题

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