用Python和rioxarray搞定MODIS数据从下载到可视化手把手教你分析科罗拉多州山火前后变化当科罗拉多州的冷泉山火在2016年肆虐时卫星遥感数据成为科学家评估灾害影响的重要工具。MODIS作为全球覆盖最频繁的中分辨率卫星传感器其地表反射率产品MOD09GA为植被变化分析提供了关键数据支撑。本文将带您用Python生态中的rioxarray和earthpy工具包完成从数据获取到灾害评估的全流程实战特别适合GIS分析师和环境科学研究者需要快速掌握遥感处理核心技能的场景。1. 环境配置与数据准备1.1 搭建Python遥感分析环境推荐使用conda创建专属的遥感分析环境避免库版本冲突conda create -n modis_analysis python3.8 conda activate modis_analysis conda install -c conda-forge rioxarray earthpy matplotlib geopandas关键库功能说明库名称核心功能版本要求rioxarray地理栅格数据处理≥0.11.0earthpy遥感专用工具集≥0.9.4geopandas矢量数据处理≥0.10.01.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], maskedTrue) # 查看数据属性 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_band1, swir_band5): 计算归一化燃烧指数 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, axax1, cmapRdYlGn, title火灾前NBR (2016-07), vmin-1, vmax1 ) # 添加图例和比例尺 fire_boundary.boundary.plot(axax1, colorblack, linewidth1) ax1.annotate(植被健康区, xy(0.1, 0.9), xycoordsaxes fraction, colordarkgreen)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], maskedTrue) nbr_post calculate_nbr(modis_post.rio.clip(clip_box) * 0.0001) # 计算dNBR dnbr nbr_pre - nbr_post4.2 烧伤严重度分级根据USGS标准划分烧伤等级dNBR范围烧伤程度RGB颜色值0.66高重度烧伤#FF00000.44-0.66中度烧伤#FFA5000.27-0.44低度烧伤#FFFF000.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, cmapburn_cmap, title冷泉山火烧伤严重度分级, cbarFalse ) # 添加自定义图例 legend_labels { 未烧伤/恢复区: #00FF00, 低度烧伤: #FFFF00, 中度烧伤: #FFA500, 高重度烧伤: #FF0000 } patches [plt.plot([],[], markers, ms10, ls, colorcolor, labellabel)[0] for label, color in legend_labels.items()] plt.legend(handlespatches, 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, dtypefloat32) # 保存分类结果 rxr.DataArray(burn_classes).rio.to_raster(burn_severity_classes.tif)通过这个完整案例我们不仅掌握了MODIS数据的处理技巧更建立起了一套可复用的山火影响评估方法。在实际项目中可以进一步结合气象数据、地形因素等进行多维分析为生态恢复决策提供更全面的数据支持。