Python使用GDAL写入新TIF文件显示全黑、大小异常问题

问题原因 #

• 文件大小差异：你生成的新文件是无附加信息的裸栅格，4633KB刚好是对应尺寸float32格式无压缩的纯像元体积。原始TIF额外包含了概览（金字塔）、像元统计值、内置元数据，部分场景还会使用压缩存储，因此体积和裸栅格存在差异，属于正常现象，不是数据写入错误。
• 显示全黑：你没有复制原始波段的NoData值，也没有为新文件计算像元统计值。ArcMap加载栅格默认按波段统计的最大最小值做对比度拉伸，你修改的像元值100远大于原始数据的正常数值范围，会导致拉伸阈值偏移，绝大多数正常像元被映射为纯黑色，看起来全图发黑。

修复方案 #

写入时补全NoData值配置、像元统计值计算，可按需配置压缩参数匹配原始文件的存储规则，修正后代码如下：

import numpy as np
from osgeo import gdal
gdal.AllRegister()

# 只读方式打开原始TIF
tif = gdal.Open('data/pre_heilj_mean90_15.tif', gdal.GA_ReadOnly)
imwidth = tif.RasterXSize
imheight = tif.RasterYSize
in_band = tif.GetRasterBand(1)
# 读取原始NoData值
nodata_val = in_band.GetNoDataValue()
data = in_band.ReadAsArray()

# 修改目标像元值
data[100][100] = 100
data = data.astype(np.float32)

# 创建输出文件，可按需添加压缩创建选项减少体积
driver = gdal.GetDriverByName("GTiff")
create_ops = ["COMPRESS=LZW", "TILED=YES"] # 若原文件用其他压缩，替换为对应参数即可
dataset = driver.Create(
    'data/res.tif',
    imwidth,
    imheight,
    1,
    gdal.GDT_Float32,
    options=create_ops
)

# 复制空间参考信息
dataset.SetSpatialRef(tif.GetSpatialRef())
dataset.SetGeoTransform(tif.GetGeoTransform())
dataset.SetProjection(tif.GetProjection())

out_band = dataset.GetRasterBand(1)
# 配置和原文件一致的NoData值
if nodata_val is not None:
    out_band.SetNoDataValue(nodata_val)
# 写入像元数组
out_band.WriteArray(data)
# 计算波段统计值，供GIS软件做正确的显示拉伸
out_band.ComputeStatistics(False)

# 刷写缓存释放资源
dataset.FlushCache()
del dataset, in_band, tif, data

补充说明 #

• 如果需要新文件和原文件体积完全一致，给新文件构建和原文件层级相同的概览（金字塔）即可，两者的体积差就是概览、元数据占用的空间。
• 如果加载后仍显示异常，可在ArcMap图层属性的符号系统选项中，手动设置拉伸的最大最小值和原文件保持一致，即可正常渲染。

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