GDAL ReadAsArray读取Landsat-8像素值异常的原因及转换方法咨询

嘿，刚接触Landsat-8和Python的话，遇到这个问题太正常啦！我来给你掰扯清楚：

为什么你读到的是数万的数值？ #

你用GDAL读到的其实是DN值（Digital Number，数字量化值），这是Landsat-8 L1级产品（就是你手里的L1TP）存储的原始定标后数据——传感器捕获的光信号经过量化后变成了0-65535之间的整数，所以出现数万的数值完全是正常的。

而《Landsat8DataUsersHandbook》里提到的0.636-0.673μm是B4波段的光谱波长范围，对应的是这个波段探测的光的波长区间；你想要的应该是地表反射率（反映地表对该波段光的反射程度，一般在0-1之间），这需要从DN值转换过来才行。

怎么用GDAL把DN值转成反射率？ #

Landsat-8的元数据文件（和你的TIFF同目录，名字是LC08_L1TP_172039_20150509_20170411_01_T1_MTL.txt）里存着所有定标需要的参数，我们可以用它来完成转换：

步骤1：读取DN值和元数据

先读取TIFF的像素数组，再解析MTL文件里的定标参数：

from osgeo import gdal
import numpy as np

# 打开B4波段的TIFF文件
tiff_path = "LC08_L1TP_172039_20150509_20170411_01_T1_B4.tiff"
ds = gdal.Open(tiff_path)
b4_dn = ds.GetRasterBand(1).ReadAsArray()

# 解析MTL元数据文件
mtl_path = tiff_path.replace("_B4.tiff", "_MTL.txt")
mtl_params = {}
with open(mtl_path, "r") as f:
    for line in f:
        if "=" in line:
            key, val = line.strip().split("=", 1)
            mtl_params[key.strip()] = val.strip().strip('"')

步骤2：用定标公式计算反射率

Landsat官方提供的反射率转换公式是：

地表反射率 = (DN × 反射率增益 + 反射率偏移) / sin(太阳高度角)

我们从元数据里提取对应的参数，代入计算：

# 获取B4波段的反射率定标系数和太阳高度角
reflect_mult = float(mtl_params["REFLECTANCE_MULT_BAND_4"])
reflect_add = float(mtl_params["REFLECTANCE_ADD_BAND_4"])
sun_elev = float(mtl_params["SUN_ELEVATION"])

# 计算反射率（注意把角度转成弧度）
b4_reflectance = (b4_dn * reflect_mult + reflect_add) / np.sin(np.radians(sun_elev))

# 把反射率限制在合理范围（0-1，超过的部分截断）
b4_reflectance = np.clip(b4_reflectance, 0, 1)

# 输出转换后的范围看看
print(f"转换后的B4波段反射率范围：{np.min(b4_reflectance):.4f} - {np.max(b4_reflectance):.4f}")

额外说明 #

• 如果需要计算辐亮度（Radiance），可以用MTL里的RADIANCE_MULT_BAND_4和RADIANCE_ADD_BAND_4，公式更简单：辐亮度 = DN × 辐亮度增益 + 辐亮度偏移
• 转换后的反射率一般在0-1之间，如果你习惯用百分比，乘以100就可以得到0-100的数值
• 一定要确保MTL文件和TIFF文件在同一个目录，路径替换要准确哦

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