用Python和rioxarray搞定MODIS数据：从下载到可视化，手把手教你分析科罗拉多州山火前后变化

当科罗拉多州的冷泉山火在2016年肆虐时，卫星遥感数据成为科学家评估灾害影响的重要工具。MODIS作为全球覆盖最频繁的中分辨率卫星传感器，其地表反射率产品（MOD09GA）为植被变化分析提供了关键数据支撑。本文将带您用Python生态中的 rioxarray 和 earthpy 工具包，完成从数据获取到灾害评估的全流程实战，特别适合GIS分析师和环境科学研究者需要快速掌握遥感处理核心技能的场景。

1. 环境配置与数据准备

1.1 搭建Python遥感分析环境

推荐使用conda创建专属的遥感分析环境，避免库版本冲突：

conda create -n modis_analysis python=3.8
conda activate modis_analysis
conda install -c conda-forge rioxarray earthpy matplotlib geopandas

关键库功能说明：

库名称	核心功能	版本要求
rioxarray	地理栅格数据处理	≥0.11.0
earthpy	遥感专用工具集	≥0.9.4
geopandas	矢量数据处理	≥0.10.0

1.2 获取冷泉山火数据集

EarthPy库内置了本次分析的示范数据集，通过以下代码自动下载：

import earthpy as et

# 下载科罗拉多州冷泉山火数据集
data_path = et.data.get_data('cold-springs-fire')
print(f"数据下载到：{data_path}")

数据集包含：

• MODIS地表反射率数据（火灾前后各一期）
• 火灾边界矢量文件（GeoMAC格式）
• 地形辅助数据

提示：首次运行会自动从Figshare下载约850MB数据，建议保持网络畅通

2. MODIS数据预处理实战

2.1 理解MOD09GA数据特性

MODIS地表反射率产品采用HDF-EOS格式存储，但示范数据已转换为GeoTIFF。关键参数需要特别注意：

import numpy as np

# MOD09GA产品关键元数据
modis_metadata = {
    "scale_factor": 0.0001,  # 缩放因子
    "fill_value": -28672,    # 无效值标记
    "valid_range": (-100, 16000),  # 有效值范围
    "bands": {
        1: "Red (620-670nm)",
        2: "NIR1 (841-876nm)",
        4: "Green (545-565nm)",
        6: "SWIR1 (1628-1652nm)"
    }
}

2.2 数据加载与质量控制

使用rioxarray进行智能加载，自动处理无效值和坐标系统：

import rioxarray as rxr
from glob import glob
import os

# 加载火灾前MODIS数据
modis_pre_path = os.path.join(data_path, "modis/reflectance/07_july_2016/crop/*.tif")
modis_pre = rxr.open_rasterio(glob(modis_pre_path)[0], masked=True)

# 查看数据属性
print(f"空间分辨率: {modis_pre.rio.resolution()}米")
print(f"坐标系: {modis_pre.rio.crs}")
print(f"无效值: {modis_pre.rio.nodata}")

2.3 空间裁剪与值域转换

针对研究区域进行空间裁剪，并应用缩放因子：

import geopandas as gpd
from shapely.geometry import box

# 加载火灾边界
fire_boundary = gpd.read_file(
    os.path.join(data_path, "vector_layers/fire-boundary-geomac/co_cold_springs_20160711_2200_dd83.shp")
)

# 坐标系统统一
fire_boundary = fire_boundary.to_crs(modis_pre.rio.crs)

# 创建裁剪范围
clip_box = [box(*fire_boundary.total_bounds)]

# 执行裁剪并缩放
modis_clipped = modis_pre.rio.clip(clip_box) * modis_metadata["scale_factor"]

3. 植被指数计算与分析

3.1 归一化燃烧指数(NBR)原理

NBR通过近红外(NIR)和短波红外(SWIR)波段计算，公式为：

$$ NBR = \frac{NIR - SWIR}{NIR + SWIR} $$

在MODIS中对应波段：

• NIR: 波段2 (841-876nm)
• SWIR: 波段6 (1628-1652nm)

3.2 Python实现NBR计算

def calculate_nbr(dataset, nir_band=1, swir_band=5):
    """计算归一化燃烧指数"""
    nir = dataset[nir_band]
    swir = dataset[swir_band]
    return (nir - swir) / (nir + swir)

# 计算火灾前后NBR
nbr_pre = calculate_nbr(modis_clipped)

3.3 结果可视化对比

使用matplotlib创建专业级专题图：

import matplotlib.pyplot as plt
import earthpy.plot as ep

# 创建对比子图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

# 火灾前NBR
ep.plot_bands(
    nbr_pre,
    ax=ax1,
    cmap="RdYlGn",
    title="火灾前NBR (2016-07)",
    vmin=-1, vmax=1
)

# 添加图例和比例尺
fire_boundary.boundary.plot(ax=ax1, color='black', linewidth=1)
ax1.annotate('植被健康区', xy=(0.1, 0.9), xycoords='axes fraction', color='darkgreen')

4. 变化检测与灾害评估

4.1 差异NBR(dNBR)计算

# 加载火灾后数据并计算NBR
modis_post_path = os.path.join(data_path, "modis/reflectance/09_sept_2016/crop/*.tif")
modis_post = rxr.open_rasterio(glob(modis_post_path)[0], masked=True)
nbr_post = calculate_nbr(modis_post.rio.clip(clip_box) * 0.0001)

# 计算dNBR
dnbr = nbr_pre - nbr_post

4.2 烧伤严重度分级

根据USGS标准划分烧伤等级：

dNBR范围	烧伤程度	RGB颜色值
>0.66	高重度烧伤	#FF0000
0.44-0.66	中度烧伤	#FFA500
0.27-0.44	低度烧伤	#FFFF00
<0.27	未烧伤/恢复区	#00FF00

实现分级可视化：

from matplotlib.colors import ListedColormap

# 自定义颜色映射
burn_colors = ['#00FF00', '#FFFF00', '#FFA500', '#FF0000']
burn_cmap = ListedColormap(burn_colors)

# 创建分类边界
bounds = [0, 0.27, 0.44, 0.66, 1]
ep.plot_bands(
    dnbr,
    cmap=burn_cmap,
    title="冷泉山火烧伤严重度分级",
    cbar=False
)

# 添加自定义图例
legend_labels = {
    '未烧伤/恢复区': '#00FF00',
    '低度烧伤': '#FFFF00',
    '中度烧伤': '#FFA500',
    '高重度烧伤': '#FF0000'
}
patches = [plt.plot([],[], marker="s", ms=10, ls="", color=color, 
                    label=label)[0] for label, color in legend_labels.items()]
plt.legend(handles=patches, bbox_to_anchor=(1.05, 1))

4.3 统计烧伤面积

# 创建掩膜分类
burn_classes = np.digitize(dnbr.values, bounds)

# 计算各类像素数量
class_counts = np.bincount(burn_classes.flatten())

# 转换为面积(假设500m分辨率)
pixel_area = 500 * 500  # 平方米
burn_areas = class_counts * pixel_area / 10000  # 转换为公顷

print(f"高重度烧伤面积: {burn_areas[3]:.1f} 公顷")
print(f"总受影响面积: {sum(burn_areas[1:]):.1f} 公顷")

在完成整个分析流程后，建议将关键结果保存为GeoTIFF以便后续使用：

# 保存dNBR结果
dnbr.rio.to_raster("cold_springs_dnbr.tif", dtype="float32")

# 保存分类结果
rxr.DataArray(burn_classes).rio.to_raster("burn_severity_classes.tif")

通过这个完整案例，我们不仅掌握了MODIS数据的处理技巧，更建立起了一套可复用的山火影响评估方法。在实际项目中，可以进一步结合气象数据、地形因素等进行多维分析，为生态恢复决策提供更全面的数据支持。