如何用Python绘制HDF4格式MODIS MCD12C1土地覆盖数据

你已经成功读取了MCD12C1的HDF4子数据集，接下来可以按照以下步骤完成土地覆盖数据的地图可视化：

1. 关联子数据集名称与数据（优化数据存储） #

当前代码将所有子数据集存在列表里，但MCD12C1包含多个土地覆盖分类体系（如IGBP全球植被分类、UMD分类等），建议用字典存储，方便后续指定要可视化的分类：

import numpy as np
import pandas as pd
import rioxarray
import rasterio as rio
from osgeo import gdal
import earthpy as et
import earthpy.plot as ep 
import matplotlib.pyplot as plt

title='c:/MCD12C1.A2022001.061.2023244164746'
# 改用字典存储，键为子数据集名称，值为对应数据
lc_data = {}
with rio.open(f'{title}.HDF') as dataset:
    crs = dataset.read_crs()
    for layer_name in dataset.subdatasets:
        print(layer_name)
        with rio.open(layer_name) as subdataset:
            # 提取子数据集的分类体系名称，比如从路径中截取
            lc_type = layer_name.split(':')[-1].split('_')[-1]
            lc_data[lc_type] = subdataset.read(1)
            # 同时保存子数据集的地理变换信息
            lc_data[f"{lc_type}_transform"] = subdataset.transform

2. 处理投影与地理范围（可选但推荐） #

MCD12C1采用正弦投影（Sinusoidal），如果需要和其他常用坐标系（如WGS84）的地图叠加，可通过rioxarray转换投影：

# 以IGBP分类为例，先将数组转为xarray.DataArray
igbp_arr = rioxarray.open_rasterio([name for name in dataset.subdatasets if 'IGBP' in name][0])
# 转换为WGS84坐标系（EPSG:4326）
igbp_wgs84 = igbp_arr.rio.reproject("EPSG:4326")

3. 定义土地覆盖分类的颜色映射与标签 #

MCD12C1的每个分类体系都有官方指定的颜色和类别名称，以IGBP分类为例，你可以定义对应的颜色映射和标签：

# IGBP土地覆盖分类的颜色（RGB格式）和标签
igbp_colors = [
    (0,0,0), (0,0,0), (0,0,0), (0,0,0), (0,0,0), (102,102,102),
    (0,204,0), (0,204,0), (0,204,0), (0,204,0), (102,204,102),
    (102,204,102), (102,204,102), (102,204,102), (204,204,0),
    (204,204,0), (204,204,0), (255,153,0), (204,102,0), (255,255,204),
    (204,204,204), (102,102,204), (0,0,204), (204,204,255), (0,102,204),
    (0,0,102), (204,204,255), (255,255,255), (255,255,255)
]
# 转换为matplotlib可用的颜色格式（0-1范围）
igbp_cmap = plt.matplotlib.colors.ListedColormap([(r/255, g/255, b/255) for r,g,b in igbp_colors])
# IGBP分类标签
igbp_labels = [
    "未分类", "水体", "常绿针叶林", "常绿阔叶林", "落叶针叶林",
    "落叶阔叶林", "混交林", "密灌丛", "疏灌丛", "木质稀树草原",
    "稀树草原", "草地", "永久湿地", "农田", "城市和建成区",
    "农田/自然植被交错带", "永久冰雪", "裸地", "稀树草原", "灌木",
    "草地", "农田", "城市", "水体", "湿地", "裸地", "冰雪", "云"
]
# 设置颜色条的边界（分类值从0到28）
igbp_bounds = np.arange(0, 29)
igbp_norm = plt.matplotlib.colors.BoundaryNorm(igbp_bounds, igbp_cmap.N)

4. 绘制可视化地图 #

使用earthpy.plot或matplotlib绘制土地覆盖图，添加颜色条、标题和地理信息：

# 用转换后的WGS84数据绘制
fig, ax = plt.subplots(figsize=(12, 8))
im = ax.imshow(igbp_wgs84[0], cmap=igbp_cmap, norm=igbp_norm)

# 添加颜色条，设置标签
cbar = ep.draw_legend(im, titles=igbp_labels, classes=igbp_bounds)
cbar.set_label('IGBP土地覆盖类型')

# 设置标题和坐标轴标签
ax.set_title('2022年MCD12C1 IGBP土地覆盖分布', fontsize=14)
ax.set_xlabel('经度')
ax.set_ylabel('纬度')

# 移除坐标轴刻度（可选）
ax.set_xticks([])
ax.set_yticks([])

plt.tight_layout()
plt.show()

额外提示 #

• 如果你不需要转换投影，直接用原始正弦投影的数据绘制即可，只需将igbp_wgs84[0]替换为lc_data['IGBP']，并注意坐标轴是像素坐标而非经纬度。
• 其他分类体系（如UMD、LAI）的颜色映射和标签可以参考MODIS官方文档自行定义。

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