从HDF到月尺度ET：基于MOD16A2的流域蒸散发数据处理全流程解析

1. MOD16A2数据简介与准备工作

MOD16A2是MODIS系列卫星产品中专门用于监测全球陆地蒸散发（Evapotranspiration, ET）的8天合成数据产品。对于水文和生态研究者来说，这套数据能提供从2000年至今、空间分辨率500米的连续观测记录。我第一次接触这个数据集是在研究黄河流域水循环时，当时就被它完整的时间覆盖和稳定的数据质量所吸引。

数据产品特点值得重点关注：每个HDF文件包含ET（实际蒸散发）、PET（潜在蒸散发）、LE（潜热通量）三个核心参数。文件命名规则像"A2018001.h24v05.006.2018012234656.hdf"这样的结构，其中"A2018001"表示2018年第1天开始的8天周期，"h24v05"对应MODIS特有的全球分块编号系统。中国大部分区域集中在h23v04到h30v06这些分块中。

在开始处理前，需要准备以下基础工具链：

• MRT（MODIS Reprojection Tool）：NASA官方提供的格式转换工具
• Python环境（推荐Anaconda）+ GDAL库
• GIS软件（QGIS或ArcGIS）
• 至少50GB的可用磁盘空间（原始数据和处理中间文件会很占空间）

提示：建议建立清晰的文件夹结构，例如：

• /raw_hdf 存放原始数据
• /reprojected 存放重投影结果
• /mosaic 存放拼接后的文件
• /final 存放最终成果

2. 数据获取与初步处理

2.1 高效下载策略

NASA的LAADS DAAC（https://ladsweb.modaps.eosdis.nasa.gov/）是官方数据源，但直接网页下载大范围长时间序列数据会很痛苦。我推荐两种更高效的方法：

方法一：使用R脚本批量下载

library(MODIStsp) MODIStsp( gui = FALSE, out_folder = "D:/MOD16A2", selprod = "ET_500m", bandsel = "ET_500m", user = "你的NASA账号", password = "你的密码", start_date = "2020.01.01", end_date = "2020.12.31", spatmeth = "tiles", h = c(26,27), v = c(04,05), download_server = "LAADS" )

方法二：Python+wget自动化

import requests from bs4 import BeautifulSoup base_url = "https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/6/MOD16A2/" years = ["2020"] tiles = ["h26v04","h26v05","h27v04","h27v05"] for year in years: for tile in tiles: url = f"{base_url}/{year}/001/" response = requests.get(url) soup = BeautifulSoup(response.text, 'html.parser') for link in soup.find_all('a'): if tile in link.get('href'): file_url = url + link.get('href') !wget --load-cookies ~/.urs_cookies --save-cookies ~/.urs_cookies --keep-session-cookies -P ./raw_hdf {file_url}

2.2 HDF文件解析技巧

下载得到的HDF文件是分层数据结构，可以用HDFView工具直观查看。以MOD16A2为例，关键数据层包括：

• ET_500m：8天实际蒸散发量
• ET_QC_500m：质量评估层
• PET_500m：潜在蒸散发量

用GDAL可以快速查看文件信息：

gdalinfo HDF4_EOS:EOS_GRID:"MOD16A2.A2020001.h26v05.006.2020006033444.hdf":MOD_Grid_MOD16A2:ET_500m

3. 数据预处理全流程

3.1 投影转换与拼接

MRT工具虽然界面老旧，但在处理MODIS数据上依然是最稳定的选择。我总结了一个高效批处理方案：

1. 先制作参数模板文件（.prm），示例内容：

INPUT_FILENAME = input.hdf SPECTRAL_SUBSET = ( 1 1 1 0 0 0 ) SPATIAL_SUBSET_TYPE = INPUT_LAT_LONG SPATIAL_SUBSET_UL_CORNER = ( 41.0 98.0 ) SPATIAL_SUBSET_LR_CORNER = ( 32.0 108.0 ) OUTPUT_FILENAME = output.tif RESAMPLING_TYPE = NEAREST_NEIGHBOR OUTPUT_PROJECTION_TYPE = UTM DATUM = WGS84 UTM_ZONE = 48 OUTPUT_PIXEL_SIZE = 500

2. 创建批处理脚本（Windows示例）：

@echo off set MRT_HOME=C:\MRT set PATH=%PATH%;%MRT_HOME%\bin for %%f in (.\raw_hdf\*.hdf) do ( java -jar %MRT_HOME%\bin\MRTBatch.jar -p modis.prm -d %%~nf.hdf -o .\reprojected\%%~nf.tif )

3.2 流域裁剪实战

当需要针对特定流域分析时，精确裁剪是关键步骤。我改进过的ArcPy脚本可以处理复杂情况：

import arcpy, os from arcpy.sa import * arcpy.CheckOutExtension("Spatial") # 设置工作空间 input_folder = r"D:\reprojected" output_folder = r"D:\clipped" basin_shp = r"D:\boundary\yangtze.shp" # 创建空间参考对象确保坐标系一致 sr = arcpy.Describe(basin_shp).spatialReference arcpy.env.outputCoordinateSystem = sr arcpy.env.extent = basin_shp # 批量处理 for tif in arcpy.ListRasters("*.tif", "TIF", input_folder): out_name = os.path.join(output_folder, "clip_" + tif) # 先裁剪再处理无效值 clipped = ExtractByMask(tif, basin_shp) # 处理MOD16特有的填充值 null_raster = SetNull( (clipped == 32761) | (clipped == 32762) | (clipped == 32763) | (clipped == 32764) | (clipped == 32765) | (clipped == 32766), clipped ) # 转换为实际物理值 calibrated = null_raster * 0.1 calibrated.save(out_name) print(f"Processed: {tif}")

4. 时间尺度转换与月值合成

4.1 8天到月尺度转换

MOD16A2的8天合成数据需要特殊处理才能转为月值。这里有个容易踩的坑：每年最后一个周期可能是5天（平年）或6天（闰年），不能简单用8天系数换算。

可靠的计算方法：

1. 先计算每日ET值：

   • 常规周期：ET_daily = ET_8day * 0.1 / 8
   • 年末周期：ET_daily = ET_5day * 0.1 / 5（或6天）

2. 按月汇总：

import numpy as np import pandas as pd from osgeo import gdal def monthly_aggregate(tif_folder, year): # 创建日期索引 dates = pd.date_range(f"{year}-01-01", f"{year}-12-31") doy_ranges = [] # MODIS 8天周期划分 for start in range(1, 366, 8): end = min(start+7, 365) doy_ranges.append((start, end)) # 读取所有TIFF文件并匹配日期 monthly_data = {} for tif in os.listdir(tif_folder): if not tif.endswith(".tif"): continue # 解析文件名中的日期信息 parts = tif.split(".") doy = int(parts[1][1:]) # 找到对应的8天周期 period = next((i for i, (s,e) in enumerate(doy_ranges) if s <= doy <= e), None) if period is None: continue # 计算权重（年末特殊处理） if period == len(doy_ranges)-1: days = 5 if not calendar.isleap(year) else 6 else: days = 8 # 读取数据 ds = gdal.Open(os.path.join(tif_folder, tif)) band = ds.GetRasterBand(1) arr = band.ReadAsArray() # 累加到对应月份 month = dates[doy-1].month if month not in monthly_data: monthly_data[month] = np.zeros_like(arr) monthly_days = np.zeros_like(arr) valid_mask = arr != -9999 monthly_data[month][valid_mask] += arr[valid_mask] * days monthly_days[valid_mask] += days # 计算月平均值 for month in monthly_data: monthly_data[month] = np.where( monthly_days > 0, monthly_data[month] / monthly_days, -9999 ) return monthly_data

4.2 结果验证与质量控制

完成计算后必须进行合理性检查：

1. 数值范围验证：典型ET值应在0-10 mm/day之间
2. 时空连续性检查：相邻月份不应出现剧烈跳变
3. 与气象站观测数据对比（如有）

我常用的验证脚本：

import matplotlib.pyplot as plt def plot_monthly_stats(monthly_data): means = [np.nanmean(data[data != -9999]) for data in monthly_data.values()] stds = [np.nanstd(data[data != -9999]) for data in monthly_data.values()] plt.figure(figsize=(10,5)) plt.bar(range(1,13), means, yerr=stds) plt.xlabel("Month") plt.ylabel("ET (mm/day)") plt.title("Monthly ET Variation") plt.grid(True) plt.show()

5. 完整流程优化建议

经过多个项目实践，我总结出几个效率提升技巧：

1. 内存优化：处理大范围数据时，使用分块处理

   chunk_size = 1024 # 分块大小 for i in range(0, rows, chunk_size): for j in range(0, cols, chunk_size): window = ((i, min(i+chunk_size, rows)), (j, min(j+chunk_size, cols))) chunk = ds.ReadAsArray(j, i, window[1][1]-window[1][0], window[0][1]-window[0][0]) # 处理分块数据

2. 并行处理：利用多核加速

   from multiprocessing import Pool def process_tif(tif_path): # 单个文件处理逻辑 pass with Pool(processes=4) as pool: pool.map(process_tif, tif_files)

3. 自动化工作流：使用Apache Airflow或Prefect构建数据处理流水线，实现从数据下载到最终产品生成的全自动运行。

在实际项目中，我发现长江流域夏季ET值明显高于理论计算值，经过反复验证发现是水稻田灌溉导致的特殊现象。这种实地观察与遥感数据的相互印证，往往能发现更有价值的研究切入点。