MODIS数据综合处理软件V1.0：从HDF4到NDVI的工程化实践

简介：面向遥感与 GIS 分析人员的 MODIS 数据综合处理软件 V1.0 安装包及配套使用手册，主要解决 NDVI/EVI、ET/PET、LST、LAI、GPP/NPP 等陆面产品批量读取、统计与可视化问题。软件支持多文件批量导入、多核并行计算，能够明显提升大批量时序数据的处理效率；内置可视化模块可快速出图，便于用户观察趋势变化与异常特征。资源打包为 zip 压缩包，共 2000 个文件，约 176.28MB，以 JavaScript、JSON、HTML/CSS 等界面与配置类文件为主，同时包含 Python 脚本、PDF 手册及 Markdown/TXT 说明文档，下载解压后即可参照使用手册完成安装与操作。资源已申请软件著作权，未经允许禁止搬运。目前已有 138 人学习下载，适合从事地理空间数据分析、植被遥感监测及生态模型应用的学生与科研人员使用，可作为 MODIS 产品处理流程的落地工具和参考模板。

1. MODIS 数据综合处理软件 V1.0 的边界与最小闭环

做 MODIS 数据综合处理的工程师，最难受的往往不是算法，而是那两三万个 HDF 文件。MODIS 数据免费、重访周期短、36 个波段，但交付的是 HDF4-EOS 格式，正弦投影按瓦片存，反射率带缩放因子，云和水汽混在信号里，一张能用的 NDVI 也要先过定标、掩膜、大气校正和重投影这几关。

MODIS 数据综合处理软件 V1.0 不是发明新算法，而是把这条链路收敛成可重复执行、批量可跑、参数可配的工具链。V1.0 的含义是边界：输入产品、输出格式、参数约定、失败判据都要先定死。功能可以少，结果必须可复现，三个月后换台机器还能跑出同一组数。

下文按产品选型、核心实现、批量调度、质量控制、验证对照的顺序展开，基于 GDAL 加 Python 这套最常见的开源组合，照着搭就能交付一个能用的 V1.0。

2. MODIS 数据综合处理的链路设计：选 L1B 还是 L2G

搭 V1.0 之前先回答一个问题：处理起点放在哪一级产品。这个决定影响后续要不要写大气校正、要不要处理扫描条带几何，甚至决定整条链路的稳定性。选错起点，后面每一步都在为前面买单。

2.1 产品分层：L1B、L2G、L3 各自承担什么

MODIS 归档产品按处理级别划分。L1B 是辐射定标后的条带观测，几何上保留扫描结构，定位信息在单独的 MOD03 里；L2G 是把逐日观测重采样到固定正弦网格并完成大气校正的格网产品；L3 是在 L2G 之上按时间合成的统计产品。先分清这三层，再看 V1.0 的输入放哪层。

从 MOD021KM 起步，等于把定标、几何定位、条带拼接全扛在自己身上，光是扫描条带边缘的弓形重叠就要写一套去重逻辑；从 MOD09GA 起步，这些环节官方已经处理过一轮，V1.0 只负责「反射率到指数」这一段。

2.2 推荐起点：MOD09GA，理由有三

第一，MODIS 卫星数据大气校正本身是一个独立工程。MOD09 内置的校正算法以 6S 为基础，配合 MOD04 气溶胶和 MOD07 水汽逐像元输入。自己复刻，要把 AOD、水汽、臭氧逐像元组装成 6S 输入再做查找表，工作量足够开第二个项目。第二，L1B 条带边缘像元重叠，处理顺序错了会出现重复扫描线；L2G 是格网数据，天然绕开这个问题。第三，V1.0 做植被应用，MOD09GA 的蓝、红、近红外、短波红外四个波段齐全，NDVI、EVI、NDWI 都不缺输入。

2.2.1 什么时候才需要回到 L1B 起点

做自定义波段组合、夜间火点、BRDF 反演时才需要 L1B。V1.0 不碰这些，就把 MOD09GA 定为唯一反射率输入。这样整条链路少一半不确定性，排错范围也小很多。

2.3 V1.0 的输入输出约定写进配置，不写进代码

代码里写死路径和分辨率，是脚本 V0.9 的典型症状。V1.0 把约定收敛到一个 yaml 文件，程序只读配置不读记忆。

# config_v1.yaml：MODIS 数据综合处理软件 V1.0 处理约定 inputs: surface_reflectance: "MOD09GA" cloud_mask: "MOD35_L2" tiles: ["h26v05", "h27v06", "h28v07"] time_range: ["2020-06-01", "2020-08-31"] output: crs: "EPSG:4326" res_degree: 0.01 format: "COG" indices: ["NDVI", "EVI", "NDWI"]

tiles 用 h/v 瓦片编号，time_range 按天闭区间，输出统一 Cloud Optimized GeoTIFF。COG 的好处是后续 Web 出图可以直接 Range 请求，不用先转瓦片。配置化的好处是调参不动代码，跑批前 diff 一下 yaml 就能知道这次和上次差在哪。

3. MODIS 综合处理管线的核心实现：定标、大气校正与指数计算

3.1 读 HDF4 并完成反射率定标

读 MODIS 的 HDF4-EOS 文件有 pyhdf 和 GDAL 两条路。pyhdf 在 Python 3.10 以上环境编译经常失败，V1.0 直接用 GDAL 的 HDF4 驱动，少一个依赖。GDAL 把 HDF4-EOS 暴露成子数据集，每个波段对应独立子集，读取时按名字匹配。

from osgeo import gdal import numpy as np def read_sur_refl(hdf_path, band="sur_refl_b01"): ds = gdal.Open(hdf_path) # 在子数据集列表中按波段名定位，避免依赖固定顺序 subs = [s for s in ds.GetSubDatasets() if band in s[0]] if not subs: raise ValueError(f"{hdf_path} 中找不到波段 {band}") sds = gdal.Open(subs[0][0]) arr = sds.ReadAsArray().astype(np.float64) # MOD09GA 反射率缩放系数固定为 0.0001，无效值 -28672 scale = float(sds.GetMetadataItem("scale_factor")) fill = float(sds.GetMetadataItem("_FillValue")) ref = np.where(arr != fill, arr * scale, np.nan) ref[(ref < 0) | (ref > 1)] = np.nan # 越界反射率一律置空 return ref

先读子数据集再乘缩放系数，顺序不能反，否则 -28672 乘 0.0001 会变成 -2.8672，混进反射率后污染后续指数。把越界值置 np.nan 是因为部分云边界像元反射率会大于 1，不处理的话 NDVI 会出现超过 1 的假值。注意 scale_factor 在 MOD021KM 这类多波段产品里是按波段存的数组，GDAL 读出来是逗号分隔字符串，多波段时要按波段索引 split 取对应项。

提示：GDAL 的 HDF4 驱动不是所有发行版都自带，报 "not recognized as a supported file format" 时先检查 GDAL 是否编译了 HDF4 支持，而不是换文件。

3.2 大气校正：MODIS 卫星数据大气校正的最小闭环

如果非要对 MODIS 卫星数据大气校正，最稳的办法是直接取 MOD09GA，官方大气校正已经做完。自己动手的场景通常是项目要求统一处理链，或需要输出指定波段的表观反射率。自行校正的常见做法分两档：DOS（暗目标）适合快速建链，6S 适合追求绝对精度。

V1.0 用的是简化 DOS。思路是假设影像里存在地表反射率接近 1%~2% 的暗目标（清洁水体或浓密植被阴影），这些像元的表观反射率近似等于大气程辐射，用它们做基线把整幅影像往下拉。

def dos_correction(toa_stack, dark_band_idx): """简化 DOS 大气校正：逐波段估计程辐射并扣除""" out = np.full_like(toa_stack, np.nan) dark_ref = toa_stack[dark_band_idx] # 取低分位像元作为暗目标，分位点可配 thr = np.percentile(dark_ref[np.isfinite(dark_ref)], 2) mask = dark_ref < thr for b in range(toa_stack.shape[0]): rpath = np.mean(toa_stack[b][mask]) # 程辐射近似 out[b] = (toa_stack[b] - rpath) / (1 - rpath) return out

dark_band_idx 在 MODIS 上通常选 band 7（2.13μm 短波红外）或 band 3（0.47μm 蓝光），前者对气溶胶更不敏感。分位点 2% 是经验值，城区或雪地场景要调到 1%。DOS 输出的绝对精度不如 6S，但相对变化够用，且实现只有十几行，适合作为 V1.0 的占位实现，后续替换 6S 查找表时接口不用动。

3.3 云掩膜取位与 NDVI 计算

3.3.1 从 MOD35_L2 提取晴空掩膜

云掩膜不能只信 MOD09GA 自带的 state 简标，我一般直接读 MOD35_L2。MOD35 的 Cloud_Mask 按 48 bit 记录每像元状态，第一个字节里的 bit 0-1 是云置信度编码：00 云、01 未定、10 基本晴空、11 高置信晴空。取位用位与操作，比字符串比较快一个量级。

def extract_clear_mask(mod35_path): ds = gdal.Open(mod35_path) subs = [s for s in ds.GetSubDatasets() if "Cloud_Mask" in s[0]] cm = gdal.Open(subs[0][0]).ReadAsArray() status = cm & 0b00000011 # 只保留最低两位 clear = (status == 0b10) | (status == 0b11) return clear.astype(np.uint8)

MOD35 是 1km 分辨率，MOD09GA 是 500m，两者做像元级联合前要把云掩膜重采样对齐到反射率网格。重采样选 nearest 保持掩膜语义，不要用 bilinear 插出灰色地带。

3.3.2 NDVI 计算的参数约定

拿到晴空掩膜后计算 NDVI：

def calc_ndvi(red, nir, clear): red = np.where(clear, red, np.nan) nir = np.where(clear, nir, np.nan) return (nir - red) / (nir + red + 1e-6)

分母加 1e-6 是防止晴空雪地像元 NIR+RED 接近 0 时除零。EVI 在此基础上加蓝光波段即可，系数固定为 G=2.5、C1=6、C2=7.5、L=1。V1.0 的指数计算统一走同一入口，便于后续加 NDWI（绿光与近红外组合）。

4. MODIS 批量综合处理 V1.0 的分块调度与质量控制

4.1 瓦片系统与重投影参数

MODIS 陆地产品按正弦投影分成 h/v 瓦片，每片 10°×10°，h 从西到东 0~35，v 从北到南 0~17，编号像 h27v06 这样按行列排布。相邻瓦片外缘有约 1% 的重叠，跨瓦片拼接前先统一重投影，再按地理位置镶嵌。

早年的 MRT（MODIS Reprojection Tool）依赖 32 位 Java 和 HDF4-EOS 旧 SDK，新环境基本跑不起来，gdalwarp 的 sinu 投影串可以完全替代，还支持 COG 直出。

# 单瓦片从 MODIS 正弦网格投影到 WGS84，输出 COG gdalwarp -overwrite \ -s_srs "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs" \ -t_srs EPSG:4326 \ -tr 0.01 0.01 -r near -of COG \ HDF4_EOS:EOS_GRID:"MOD09GA.A2020185.h27v06.061.2020202112345.hdf":MODIS_Grid_500m_2D:sur_refl_b01 \ sur_refl_b01_h27v06.tif

-s_srs 是 MODIS Sinusoidal 的投影串，地球半径 6371007.181 米，GDAL 3 以后也可以写 EPSG:6842，但完整串兼容性更好；-tr 0.01 表示输出 0.01 度分辨率；重采样选 near 保证类别数据不变形，反射率数据建议按应用换 bilinear 或 cubic；-of COG 依赖 GDAL 3.1+，老版本先输出 GTiff 再转。

注意：漏写 -s_srs 会导致 HDF4-EOS 网格被当作 WGS84 去投影，输出影像整体拉扁错位。批量脚本里这条参数建议从配置读取，不要手敲。

4.2 批量并行与断点续跑

MODIS 综合处理的批量任务天然适合进程级并行：单瓦片单文件互不依赖。用 ProcessPoolExecutor 比 multiprocessing 的 starmap 好维护，异常能透传，进度可打点。关键设计是输出存在才跳过，配合先写临时文件再原子改名，中断重启不产生半成品。

from concurrent.futures import ProcessPoolExecutor, as_completed from pathlib import Path def process_one(hdf_path, out_dir, mod35_path): stem = Path(hdf_path).stem out_tif = Path(out_dir) / f"{stem}_ndvi.tif" if out_tif.exists(): return f"skip {stem}" # 断点续跑核心 red = read_sur_refl(str(hdf_path), "sur_refl_b01") nir = read_sur_refl(str(hdf_path), "sur_refl_b02") clear = extract_clear_mask(str(mod35_path)) ndvi = calc_ndvi(red, nir, clear) write_cog(ndvi, out_tif) return f"done {stem}" files = sorted(Path("/data/modis").glob("MOD09GA.*.hdf")) with ProcessPoolExecutor(max_workers=8) as pool: futures = [pool.submit(process_one, f, OUT_DIR, MOD35_DIR) for f in files] for fut in as_completed(futures): print(fut.result(), flush=True)

max_workers 一般取 CPU 物理核数，GDAL 内部开了多线程就不要叠太多进程，否则 IO 反而成瓶颈。每完成一个瓦片打印一行，配合日志做故障定位；跑批中途崩了，重跑一遍会自动跳过已完成瓦片。

4.3 质量控制指标与触发动作

批量处理最怕「看起来跑完了」。V1.0 把质控收敛成三个数值指标，写进每瓦片的 JSON 报告：有效像元比例、云覆盖率、NDVI 越界像元比例。

def qc_report(ndvi, clear, min_valid=0.5): valid = np.isfinite(ndvi) cloud_ratio = 1.0 - clear.mean() bad_ratio = ((ndvi < -1) | (ndvi > 1) & valid).mean() ok = valid.mean() >= min_valid and cloud_ratio <= 0.7 return { "valid_ratio": round(float(valid.mean()), 4), "cloud_ratio": round(float(cloud_ratio), 4), "out_of_range": round(float(bad_ratio), 4), "pass": bool(ok), }

4.3.1 三个指标怎么用

有效像元比例低于阈值，多是瓦片处于水域或云掩膜过严；云覆盖率偏高，说明该日影像不适合参与时间序列合成；越界像元比例大于 0.001，基本是定标阶段 fill 没拦干净。阈值按区域和应用调整，默认值写进配置，草地和林地场景差异很大。

4.3.2 日志与可复现性

md5sum config_v1.yaml >> run.log

每次跑批把配置散列值记进日志，配合瓦片输出时间戳，任何一次产物的参数来源都可回溯。这是 V1.0 可复现性最便宜的实现方式，比写十页技术文档管用。

5. MODIS 综合处理结果的验证与三个常见配置陷阱

5.1 与官方产品逐像素对比

验证不是肉眼看两张图差不多，而是把输出和官方 MOD13 对齐到同一坐标网格后做统计。官方 MOD13Q1 是 16 天合成，直接拿单日 NDVI 对比没有意义，正确做法是先把逐日 NDVI 做最大值合成，再和 MOD13A1（8 天合成）对比。

def compare_ndvi(ours, official, min_pixels=1000): mask = np.isfinite(ours) & np.isfinite(official) if mask.sum() < min_pixels: return {"skip": True, "n": int(mask.sum())} diff = ours[mask] - official[mask] return { "rmse": round(float(np.sqrt(np.mean(diff**2))), 4), "bias": round(float(np.mean(diff)), 4), "n": int(mask.sum()), }

经验阈值：RMSE 在 0.05 以内说明链路正常，0.05~0.1 需要检查云掩膜与重投影参数，超过 0.1 优先怀疑 fill 处理或投影串错误。

5.2 三个常见配置陷阱

陷阱一：忘记 _FillValue 就做算术。反射率数组里混着 -28672，NDVI 公式会给出一批接近 -0.99 的假像元，质控越界指标立刻报警。陷阱二：scale_factor 是多波段数组时只取第一项。MOD021KM 这类产品缩放系数按波段不同，取错波段就错一个数量级。陷阱三：重投影漏写 -s_srs，整景影像拉扁，目检难以发现，和官方产品对比时 bias 会系统性偏移。

自查办法是固定一个文件：取 h27v06 晴空影像，把云覆盖率、NDVI 直方图、与 MOD13 的 RMSE 三个数写进配置注释，后续每次改代码只跑这一个文件对照三数，有回归立刻发现，不必每轮全量重跑。

本文还有配套的精品资源，点击获取