要使用 Python 与 GDAL 从 MOD13C2 数据中批量提取 NDVI 和 EVI 并保存为 GeoTIFF 格式，核心在于理解 HDF 文件结构、应用正确的缩放因子、处理无效值、并正确设置地理参考。以下是基于多个参考资料的完整方案和代码。 ### **核心流程与方案对比** 在进行批量转换前，我们需要明确 MOD13C2 数据的特点。MOD13C2 是全球月度的植被指数产品，其 HDF 文件内包含多个子数据集。下表对比了两种主要处理方案的核心差异： | 方案 | 核心工具库 | 适用场景 | 优点 | 缺点 | | :--- | :--- | :--- | :--- | :--- | | **纯 GDAL 方案** | `gdal`, `numpy` | 需要精细控制、定制化处理流程、或集成到大型自动化脚本中。 | 灵活性极高，可完全控制每个步骤（如缩放、掩膜、投影）。 | 代码量相对较大，需要用户对 GDAL 和 MODIS 数据结构有较好理解 [ref_4]。 | | **rasterio 简化方案** | `rasterio`, `numpy` | 追求代码简洁、现代 Python 风格，处理流程相对标准。 | 语法简洁，与 `numpy` 结合紧密，易于读写和操作数组数据。 | 对复杂投影变换或 HDF 多层结构的直接支持不如 GDAL 原生。 | 以下将详细介绍**纯 GDAL 方案**的完整实现，因为它提供了最全面的控制，是处理 MODIS HDF 数据的基石 [ref_4]。 ### **代码实现：基于 GDAL 的批量转换** 此代码实现了从指定文件夹读取所有 MOD13C2 HDF 文件，提取 NDVI 和 EVI 图层，应用缩放和无效值处理，并输出为 GeoTIFF 文件。 ```python import os import glob from osgeo import gdal, osr import numpy as np def mod13c2_hdf_to_geotiff(input_folder, output_folder, vi_type='NDVI'): """ 批量将 MOD13C2 HDF 文件中的 NDVI 或 EVI 转换为 GeoTIFF。 参数: input_folder (str): 存放输入 HDF 文件的文件夹路径。 output_folder (str): 输出 GeoTIFF 文件的文件夹路径。 vi_type (str): 要提取的植被指数类型，'NDVI' 或 'EVI'。 """ # 1. 创建输出文件夹 os.makedirs(output_folder, exist_ok=True) # 2. 定义 MOD13C2 中 NDVI 和 EVI 的子数据集名称 [ref_1] subdataset_dict = { 'NDVI': 'HDF4_EOS:EOS_GRID:"{file_path}":MOD_CMG_VI:CMG 0.05 Deg Monthly NDVI', 'EVI': 'HDF4_EOS:EOS_GRID:"{file_path}":MOD_CMG_VI:CMG 0.05 Deg Monthly EVI' } # 3. MODIS 植被指数的缩放因子和无效值定义 [ref_1] SCALE_FACTOR = 0.0001 FILL_VALUE = -3000 # 有效值范围（缩放前） VALID_RANGE = (-2000, 10000) # 4. 查找所有 HDF 文件 hdf_files = glob.glob(os.path.join(input_folder, '*.hdf')) if not hdf_files: print(f"在文件夹 {input_folder} 中未找到 HDF 文件。") return print(f"找到 {len(hdf_files)} 个 HDF 文件，开始处理...") for hdf_path in hdf_files: try: # 5. 构建完整的子数据集路径 subdataset_path = subdataset_dict[vi_type].format(file_path=hdf_path) # 6. 使用 GDAL 打开子数据集 [ref_4] src_ds = gdal.Open(subdataset_path, gdal.GA_ReadOnly) if src_ds is None: print(f"无法打开子数据集: {subdataset_path}") continue # 7. 读取数据为 NumPy 数组 band = src_ds.GetRasterBand(1) data = band.ReadAsArray().astype(np.float32) # 8. 应用缩放因子并处理无效值 [ref_1, ref_4] # 将填充值设为 NaN data[data == FILL_VALUE] = np.nan # 应用缩放因子得到真实的 NDVI/EVI 值 data = data * SCALE_FACTOR # 可选：根据有效范围进行额外掩膜 # data[(data < VALID_RANGE[0] * SCALE_FACTOR) | (data > VALID_RANGE[1] * SCALE_FACTOR)] = np.nan # 9. 准备输出文件信息 base_name = os.path.splitext(os.path.basename(hdf_path))[0] output_path = os.path.join(output_folder, f"{base_name}_{vi_type}.tif") # 10. 获取输入数据的仿射变换参数和投影信息 geo_transform = src_ds.GetGeoTransform() projection = src_ds.GetProjection() # 11. 创建输出 GeoTIFF 文件 [ref_2] rows, cols = data.shape driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create(output_path, cols, rows, 1, gdal.GDT_Float32) # 12. 设置地理参考和投影 out_ds.SetGeoTransform(geo_transform) out_ds.SetProjection(projection) # 13. 将处理后的数据写入输出文件 out_band = out_ds.GetRasterBand(1) out_band.WriteArray(data) out_band.SetNoDataValue(np.nan) # 设置 NoData 值为 NaN # 14. 设置波段描述和统计信息（可选） out_band.SetDescription(f'{vi_type} from MOD13C2') out_band.FlushCache() out_band.ComputeStatistics(False) # 不计算忽略 NaN 值的统计 # 15. 关闭数据集，释放资源 src_ds = None out_ds = None print(f"成功处理: {os.path.basename(hdf_path)} -> {os.path.basename(output_path)}") except Exception as e: print(f"处理文件 {hdf_path} 时出错: {e}") continue print("批量转换完成！") # ====== 使用示例 ====== if __name__ == "__main__": # 设置你的输入和输出文件夹路径 input_dir = r"D:\MODIS\MOD13C2\Raw" # 替换为你的 HDF 文件所在路径 output_dir_ndvi = r"D:\MODIS\MOD13C2\GeoTIFF_NDVI" # NDVI 输出路径 output_dir_evi = r"D:\MODIS\MOD13C2\GeoTIFF_EVI" # EVI 输出路径 # 批量提取 NDVI print("开始提取 NDVI...") mod13c2_hdf_to_geotiff(input_dir, output_dir_ndvi, vi_type='NDVI') # 批量提取 EVI print("\n开始提取 EVI...") mod13c2_hdf_to_geotiff(input_dir, output_dir_evi, vi_type='EVI') ``` ### **关键步骤与技术要点解析** 1. **子数据集路径构造**：MODIS HDF 文件采用 HDF-EOS 格式，数据存储在特定的“EOS_GRID”结构中。代码中的 `subdataset_dict` 定义了访问 MOD13C2 中 NDVI 和 EVI 栅格层的标准 GDAL 子数据集路径 [ref_4]。 2. **缩放与无效值处理**：这是保证数据科学性的核心。原始数据是带缩放因子的整型（如 `int16`）。根据 MODIS 官方文档（参考资料中提及），NDVI/EVI 的缩放因子（Scale Factor）为 0.0001，无效值（Fill Value）为 -3000 [ref_1]。代码先将无效值替换为 `np.nan`，再乘以缩放因子得到实际的浮点型指数值（范围通常在 -0.2 到 1.0 之间）。 3. **地理信息继承**：通过 `GetGeoTransform()` 和 `GetProjection()` 方法从源数据中获取仿射变换参数（定义像元大小和左上角坐标）和投影信息（MOD13C2 通常为全球正弦投影），并完整地写入到输出的 GeoTIFF 中，确保结果文件具有正确的地理定位 [ref_2]。 4. **批量与健壮性**：代码使用 `glob` 模块进行文件遍历，并使用了 `try-except` 块来捕获和处理单个文件处理过程中的异常（如文件损坏、路径错误），确保一个文件的失败不会中断整个批量任务。 ### **进阶处理与优化建议** * **内存优化**：对于非常大的全球数据，可以使用 GDAL 的分块读取功能（`ReadAsArray(xoff, yoff, xsize, ysize)`）来避免一次性将全部数据读入内存 [ref_2]。 * **投影转换**：如果你需要将数据转换为其他投影（如 WGS84），可以在写入输出文件前，使用 `gdal.Warp()` 函数进行重投影。 * **质量控制**：MOD13C2 同样提供质量控制（QC）图层。你可以扩展代码，读取 QC 层，根据位掩码（bitmask）过滤掉低质量（如云、雪、冰覆盖）的像元，生成更干净的产品 [ref_2]。 * **使用 HEG 工具**：对于追求稳定性和官方流程的用户，NASA 提供的 **HDF-EOS to GeoTIFF Conversion Tool (HEG)** 是另一个强大的选择。它可以图形化或命令行批处理，专门用于转换 MODIS 数据，并能方便地进行重投影和子区裁剪 [ref_6]。Python 可以调用其命令行工具实现自动化封装。 ### **环境配置与依赖** 运行上述代码需要安装以下 Python 库： ```bash pip install numpy gdal ``` 请注意，`gdal` 的安装在不同系统上可能略有不同，有时需要使用 `conda install gdal` 或寻找预编译的 whl 文件。此代码在正确配置 GDAL 环境后，能够高效、自动化地完成 MOD13C2 数据从 HDF 到 GeoTIFF 的批量转换任务，为后续的时间序列分析、制图或模型输入做好准备 [ref_4][ref_5]。