如何用Rasterio生成可在系统查看器正常显示的NDVI伪彩色GeoTIFF

问题描述 #

计算NDVI数据集并将其从[-1,1]缩放至[0,255]，希望保存为可在普通图像查看器（如Gwenview）中显示伪彩色的GeoTIFF，但生成的文件始终显示为灰度图。因目标应用不支持带地理参考的PNG，必须使用TIFF格式。

现有代码 #

sds = rasterio.open(ds.GetSubDatasets()[c.DS_10m][0])
bands = sds.read([c.BAND_NIR, c.BAND_RED])
ndvi = np.zeros(bands[0].shape, dtype=rasterio.float32)

ndvi = scaleOnes(
    (bands[0].astype(float) - bands[1].astype(float)) / (bands[0] + bands[1])
)

kwargs = sds.meta

kwargs.update(driver="GTIFF", dtype=rasterio.uint8, count=1, compress="lzw")

with rasterio.open("ndvitest2.tif", "w", **kwargs) as dst:
    dst.write(ndvi, indexes=1)
    dst.write_colormap(1, c.COLORMAP_NDVI)

颜色映射表 #

COLORMAP_NDVI = {
    102: (191, 191, 191, 255),
    114: (173, 173, 173, 255),
    127: (255, 255, 224, 255),
    130: (255, 249, 204, 255),
    133: (237, 232, 181, 255),
    137: (222, 217, 156, 255),
    140: (204, 199, 130, 255),
    143: (189, 184, 107, 255),
    146: (176, 194, 97, 255),
    149: (163, 204, 89, 255),
    153: (145, 191, 82, 255),
    159: (128, 179, 71, 255),
    165: (112, 163, 64, 255),
    172: (97, 150, 54, 255),
    178: (79, 138, 46, 255),
    184: (64, 125, 36, 255),
    191: (48, 110, 28, 255),
    197: (33, 97, 18, 255),
    204: (15, 84, 10, 255),
}

问题根源与解决步骤 #

1. 数据类型与缩放范围不匹配 #

现有代码中scaleOnes生成的NDVI为float32类型，缩放范围未对齐颜色映射表的键区间（102-204），且未显式转换为uint8，导致查看器无法正确关联调色板。

修正代码：

# 替换原NDVI计算与缩放逻辑
ndvi_raw = (bands[0].astype(float) - bands[1].astype(float)) / (bands[0] + bands[1])
ndvi_raw = np.nan_to_num(ndvi_raw, nan=-1)  # 处理除以0产生的NaN
ndvi_raw = np.clip(ndvi_raw, -1, 1)  # 限制值在标准NDVI范围内
# 将[-1,1]映射到颜色表对应的102-204区间，转为uint8
ndvi_scaled = ((ndvi_raw + 1) * (204 - 102) / 2 + 102).astype(np.uint8)

2. 颜色映射表不完整 #

当前颜色表仅包含离散值，普通图像查看器大多不支持自动插值缺失颜色，会默认显示为灰度。需补全0-255全范围的颜色映射。

补全代码：

from scipy.interpolate import interp1d

# 提取现有颜色表的键和颜色值
keys = sorted(c.COLORMAP_NDVI.keys())
colors = np.array([c.COLORMAP_NDVI[k] for k in keys])

# 创建各颜色通道的线性插值函数
r_interp = interp1d(keys, colors[:,0], kind='linear')
g_interp = interp1d(keys, colors[:,1], kind='linear')
b_interp = interp1d(keys, colors[:,2], kind='linear')
a_interp = interp1d(keys, colors[:,3], kind='linear')

# 生成完整的0-255颜色映射
full_colormap = {}
for i in range(256):
    if i in c.COLORMAP_NDVI:
        full_colormap[i] = c.COLORMAP_NDVI[i]
    else:
        # 插值后确保颜色值在0-255范围内
        r = int(np.clip(r_interp(i), 0, 255))
        g = int(np.clip(g_interp(i), 0, 255))
        b = int(np.clip(b_interp(i), 0, 255))
        a = int(np.clip(a_interp(i), 0, 255))
        full_colormap[i] = (r, g, b, a)

3. 元数据缺失 photometric 标记 #

GeoTIFF需明确设置photometric='PALETTE'，告知查看器这是调色板图像，而非灰度图。

更新元数据：

kwargs.update(driver="GTIFF", dtype=rasterio.uint8, count=1, compress="lzw", photometric='PALETTE')

修正后的完整代码 #

import numpy as np
import rasterio
from scipy.interpolate import interp1d

# 读取遥感数据
sds = rasterio.open(ds.GetSubDatasets()[c.DS_10m][0])
bands = sds.read([c.BAND_NIR, c.BAND_RED])

# 计算原始NDVI并处理异常值
ndvi_raw = (bands[0].astype(float) - bands[1].astype(float)) / (bands[0] + bands[1])
ndvi_raw = np.nan_to_num(ndvi_raw, nan=-1)
ndvi_raw = np.clip(ndvi_raw, -1, 1)

# 映射到颜色表对应的数值区间并转换类型
ndvi_scaled = ((ndvi_raw + 1) * (204 - 102) / 2 + 102).astype(np.uint8)

# 补全颜色映射表
keys = sorted(c.COLORMAP_NDVI.keys())
colors = np.array([c.COLORMAP_NDVI[k] for k in keys])
r_interp = interp1d(keys, colors[:,0], kind='linear')
g_interp = interp1d(keys, colors[:,1], kind='linear')
b_interp = interp1d(keys, colors[:,2], kind='linear')
a_interp = interp1d(keys, colors[:,3], kind='linear')

full_colormap = {}
for i in range(256):
    if i in c.COLORMAP_NDVI:
        full_colormap[i] = c.COLORMAP_NDVI[i]
    else:
        r = int(np.clip(r_interp(i), 0, 255))
        g = int(np.clip(g_interp(i), 0, 255))
        b = int(np.clip(b_interp(i), 0, 255))
        a = int(np.clip(a_interp(i), 0, 255))
        full_colormap[i] = (r, g, b, a)

# 更新输出元数据
kwargs = sds.meta
kwargs.update(driver="GTIFF", dtype=rasterio.uint8, count=1, compress="lzw", photometric='PALETTE')

# 写入伪彩色GeoTIFF
with rasterio.open("ndvitest2.tif", "w", **kwargs) as dst:
    dst.write(ndvi_scaled, indexes=1)
    dst.write_colormap(1, full_colormap)

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