ArcGIS Pro + Python：三步搞定MODIS ET数据的流域月均值计算（避坑MRT）

ArcGIS Pro Python三步搞定MODIS ET数据的流域月均值计算避坑MRT在GIS领域MODIS 16A2蒸散发ET数据是研究流域水文循环的重要数据源。然而传统处理流程往往需要在MRT、ArcMap和Python IDE之间反复切换不仅效率低下还容易出错。本文将介绍如何利用ArcGIS Pro内置的Python环境ArcPy一站式完成从原始HDF文件到流域月均值ET的全流程计算彻底告别繁琐的工具链切换。1. 环境准备与数据导入1.1 ArcGIS Pro环境配置确保已安装ArcGIS Pro 3.0及以上版本并启用以下扩展模块Spatial Analyst用于栅格计算和空间分析Image Analyst处理遥感影像数据在Python环境中需要检查以下关键库import arcpy from arcpy.sa import * import os import numpy as np1.2 MODIS数据直接读取技巧传统MRT处理的第一步是将HDF转换为GeoTIFF而ArcGIS Pro可以直接读取HDF中的科学数据集# 设置工作空间 arcpy.env.workspace 输入HDF文件夹路径 hdf_files arcpy.ListFiles(*.hdf) # 提取ET科学数据集 for hdf in hdf_files: # 获取HDF中的子数据集 subdatasets arcpy.ListSubDatasets(hdf) et_dataset [s for s in subdatasets if ET_500m in s][0] # 转换为临时栅格 temp_raster arcpy.MakeRasterLayer_management(et_dataset[0], ftemp_{hdf[:-4]})提示MOD16A2的ET数据存储在HDF4_EOS:EOS_GRID:文件名:MOD_Grid_MOD16A2:ET_500m路径下2. 核心处理流程2.1 时空基准统一处理MODIS数据需要统一到目标投影坐标系如WGS 84 UTM同时处理中国区域的特殊值# 定义投影转换参数 target_sr arcpy.SpatialReference(32650) # WGS 84 UTM Zone 50N # 批量投影转换与空值处理 output_rasters [] for raster in temp_rasters: # 投影转换 projected arcpy.ProjectRaster_management(raster, memory/projected, target_sr) # 处理特殊值城市、水体等 null_raster SetNull(projected, projected, VALUE 32760) # 应用比例因子0.1 for MOD16A2 scaled_raster Float(null_raster) * 0.1 output_rasters.append(scaled_raster)2.2 流域裁剪的优化方案传统方法需要先生成流域面数据而ArcGIS Pro可以直接使用栅格函数实现精确裁剪# 加载流域边界 basin_boundary 流域面数据.shp # 创建裁剪函数链 def clip_to_basin(raster): clipped ExtractByMask(raster, basin_boundary) # 自动适配流域范围 arcpy.env.extent basin_boundary arcpy.env.snapRaster raster return clipped basin_rasters [clip_to_basin(r) for r in output_rasters]3. 月均值计算的关键细节3.1 8天合成到月数据的转换处理时间序列时需要特别注意闰年和月末不足8天的情况# 创建日期映射字典 year 2023 # 示例年份 month_days { 1: range(1, 32), 2: range(32, 60) if year%4 else range(32, 61), # ...其他月份 } # 按月份分组栅格 monthly_groups {} for month, day_range in month_days.items(): monthly_groups[month] [r for r in basin_rasters if int(r.name.split(A)[1][:3]) in day_range]3.2 加权平均计算考虑不同时间段的天数权重特别是年末可能不足8天的情况# 计算月均值 monthly_results {} for month, rasters in monthly_groups.items(): # 计算每个栅格的权重 weights [8 if i len(rasters)-1 else (5 if (year%4 and month12) else 3) for i in range(len(rasters))] # 加权平均计算 weighted_sum sum([Raster(r) * w for r, w in zip(rasters, weights)]) monthly_mean weighted_sum / sum(weights) monthly_results[month] monthly_mean4. 成果输出与可视化4.1 批量输出与元数据记录# 创建输出文件夹 output_dir 月结果输出 os.makedirs(output_dir, exist_okTrue) # 保存结果并记录元数据 for month, raster in monthly_results.items(): output_path f{output_dir}/ET_{year}_{month:02d}.tif raster.save(output_path) # 添加元数据 arcpy.AddMetadata_management(output_path, [ (数据来源, MOD16A2), (处理方法, ArcGIS Pro直接处理), (时间分辨率, 月均值), (单位, mm/day) ])4.2 动态可视化技巧利用ArcGIS Pro的时间滑块功能实现动态展示将所有月结果加载到地图右键图层选择时间属性设置时间字段为文件名的年月信息启用时间滑块动画对于更专业的可视化可以使用以下Python代码生成时序图表import matplotlib.pyplot as plt # 提取流域平均ET值 monthly_values { month: float(arcpy.GetRasterProperties_management(raster, MEAN).getOutput(0)) for month, raster in monthly_results.items() } # 绘制趋势图 plt.figure(figsize(10, 5)) plt.plot(list(monthly_values.keys()), list(monthly_values.values()), o-) plt.title(f{year}年流域月均ET变化) plt.xlabel(月份) plt.ylabel(ET (mm/day)) plt.grid(True)在实际项目中这套方法将传统需要2-3天的手动处理流程缩短到1小时内完成。特别是在处理多年数据时只需简单修改年份参数即可批量运行大幅提升了科研效率。