用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_post4.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数据的处理技巧，更建立起了一套可复用的山火影响评估方法。在实际项目中，可以进一步结合气象数据、地形因素等进行多维分析，为生态恢复决策提供更全面的数据支持。
