如何为栅格数据中的特定地类分配固定对应颜色

问题背景 #

我需要展示某区域的土地利用类型，现有多个包含地类信息的栅格文件，每个像素值对应特定地类。已定义地类值、对应颜色及文本标签：

from matplotlib.colors import ListedColormap
import rasterio
import numpy as np
import earthpy.plot as ep
import matplotlib.pyplot as plt

cmap_values = [0, 11, 22, 33, 40, 55, 66, 77]
cmap_colors = ['white', ## 0: 无数据
               'black', ## 11: 城镇用地
               'darkorange', ## 22: 耕地
               'brown', ## 33: 牧场
               'darkgreen', ## 40: 林地
               'purple', ## 55: 草地/灌丛
               'gray', ## 66: 其他土地
               'blue'  ## 77: 水域
              ] 
cmap = ListedColormap(cmap_colors)
value_text = ['无数据', '城镇用地', '耕地', '牧场', 
             '林地', '草地/灌丛', 
             '其他土地', '水域']

但部分栅格文件缺少部分地类，例如某栅格的唯一像素值为：

src = rasterio.open('myFile.tif')
data = src.read(1)
print(np.unique(data))
# 输出：array([11, 22, 33, 40, 55, 77], dtype=uint8)

使用以下代码可视化时，颜色与地类对应错误（如城镇用地被显示为白色）：

f, ax = plt.subplots()
im = ax.imshow(data, cmap=cmap)
ep.draw_legend(im, titles=value_text, classes=cmap_values)

解决方案 #

问题根源是ListedColormap按索引匹配颜色，而imshow会自动将数据的最小值映射到色图第一个颜色、最大值映射到最后一个，导致缺失地类时颜色错位。需明确指定数值与颜色的对应关系，以下是两种可行方法：

方法1：用BoundaryNorm绑定数值与颜色 #

通过BoundaryNorm定义每个地类值的区间，确保特定数值对应固定颜色，不受数据缺失影响：

from matplotlib.colors import BoundaryNorm
from matplotlib.patches import Patch

# 为每个地类值创建边界区间（避免数值落在区间外）
bounds = [val - 0.5 for val in cmap_values] + [cmap_values[-1] + 0.5]
# 创建Norm对象，绑定边界与色图
norm = BoundaryNorm(bounds, cmap.N)

# 可视化时传入norm参数
f, ax = plt.subplots()
im = ax.imshow(data, cmap=cmap, norm=norm)

# 筛选当前数据存在的地类，生成对应图例
present_values = np.unique(data)
present_indices = [cmap_values.index(val) for val in present_values]
legend_elements = [
    Patch(facecolor=cmap_colors[i], label=value_text[i]) 
    for i in present_indices
]
ax.legend(handles=legend_elements, loc='lower right', bbox_to_anchor=(1.2, 0))

plt.show()

这种方法适合需要统一所有栅格颜色规则，或保留完整地类体系的场景。

方法2：过滤色图与标签，仅保留存在的地类 #

如果只需要展示当前栅格存在的地类，可先筛选对应颜色和标签，再创建新色图：

# 获取当前数据包含的地类值
present_values = np.unique(data)
# 匹配对应的颜色和文本标签
filtered_colors = [cmap_colors[cmap_values.index(val)] for val in present_values]
filtered_texts = [value_text[cmap_values.index(val)] for val in present_values]
# 创建仅包含当前地类的色图
filtered_cmap = ListedColormap(filtered_colors)

# 可视化并绘制对应图例
f, ax = plt.subplots()
im = ax.imshow(data, cmap=filtered_cmap)
ep.draw_legend(im, titles=filtered_texts, classes=present_values)

plt.show()

该方法更简洁，适合聚焦当前栅格地类的场景。

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