基于 GEE 的长时间序列 MODIS NDVI 年度均值批量提取与导出

最新推荐文章于 2025-12-17 01:19:38 发布

原创最新推荐文章于 2025-12-17 01:19:38 发布 · 744 阅读

18 ·

CC 4.0 BY-SA版权

文章标签：

#GEE #云计算 #云平台 #遥感大数据 #数据集 #MODIS #NDVI

Google Earth Engine（GEE）专栏收录该内容

154 篇文章

订阅专栏

该文章已生成可运行项目，

一、AOI（感兴趣区域）配置

二、全局参数配置

三、MODIS 数据加载与 NDVI 计算

若觉得代码对您的研究 / 项目有帮助，欢迎点击打赏支持！需要完整代码的朋友，打赏后可在后台私信（复制文章标题发给我），我会尽快发您完整可运行代码，感谢支持！

本代码基于 Google Earth Engine（GEE）平台，实现2004-2024 年 MODIS 卫星数据的 NDVI 年度均值批量计算与导出。NDVI（归一化植被指数）是反映植被覆盖度和生长状况的关键指标，代码通过处理 MODIS 反射率数据，筛选有效像元、按年度聚合均值，并将结果导出至 Google Drive，同时支持指定区域（AOI）的针对性分析。

一、AOI（感兴趣区域）配置

var aoi = table;
Map.centerObject(aoi, 8);
Map.addLayer(aoi, {color: 'red'}, 'AOI');

AOI 定义：var aoi = table 中，table 是 GEE 平台中已上传的矢量数据集（如行政区划、自定义研究区域），需用户对该数据集有访问权限；若无权限，需替换为公共矢量数据（如国家边界、经纬度范围生成的几何对象）。
地图定位与可视化：
- Map.centerObject(aoi, 8)：将地图视图中心定位到 AOI 区域，缩放级别为 8（GEE 缩放级别 0-20，数值越大视图越精细）。
- Map.addLayer(aoi, {color: 'red'}, 'AOI')：在地图上添加 AOI 图层，以红色显示，图层名称为 “AOI”，方便用户直观确认研究区域。

二、全局参数配置

var EXPORT_FOLDER = 'GEE_exports';
var YEAR_START = 2004;
var YEAR_END   = 2024;
var YEARS = ee.List.sequence(YEAR_START, YEAR_END).getInfo();

导出路径配置：EXPORT_FOLDER = 'GEE_exports' 定义导出文件在 Google Drive 中的存储文件夹名称，若文件夹不存在，GEE 会自动创建。
时间范围配置：
- YEAR_START 和 YEAR_END 分别指定时间序列的起始年份（2004）和结束年份（2024）。
- ee.List.sequence(...) 生成从起始年到结束年的连续年份列表（如 [2004,2005,...,2024]），getInfo() 将 GEE 服务器端的列表转换为客户端可迭代的数组，用于后续年度循环。

三、MODIS 数据加载与 NDVI 计算

var modisAll = ee.ImageCollection('MODIS/061/MOD09A1')
  .filterBounds(aoi)
  .map(function(img) {
    var nir = img.select('sur_refl_b02');  // 近红外波段
    var red = img.select('sur_refl_b01');  // 红波段
    
    // 计算NDVI并命名
    var ndvi = ee.Image(nir.subtract(red).divide(nir.add(red))).rename('NDVI');
    // 多步掩膜处理，保留有效像元
    ndvi = ndvi.updateMask(ndvi.gte(-0.2));
    ndvi = ndvi.updateMask(ndvi.lte(1));
    var validSum = nir.add(red).gt(0);
    ndvi = ndvi.updateMask(validSum);
    return ndvi.copyProperties(img, ['system:time_start']);
  });

数据加载：
- ee.ImageCollection('MODIS/061/MOD09A1') 加载 MODIS 卫星的 061 版本 MOD09A1 数据集，该数据集是 8 天合成的地表反射率产品，空间分辨率 500 米，包含多个波段（红波段、近红外波段等），适合 NDVI 计算。
- filterBounds(aoi) 筛选数据集在 AOI 区域内的影像，减少无关数据处理量，提升效率。
NDVI 计算逻辑：
- NDVI 核心公式：(近红外波段 - 红波段) / (近红外波段 + 红波段)，植被对近红外波段反射强、红波段吸收强，因此 NDVI 值越高，植被覆盖越好。
- nir = img.select('sur_refl_b02') 和 red = img.select('sur_refl_b01') 分别提取近红外波段（sur_refl_b02）和红波段（sur_refl_b01），这是 MOD09A1 数据集的固定波段标识。
- nir.subtract(red) 计算近红外与红波段的差值，nir.add(red) 计算两者的和，再通过 divide() 得到 NDVI 值，最后用 rename('NDVI') 为结果波段命名，方便后续处理。
掩膜（Mask）处理（关键步骤）：
- 目的：剔除无效 NDVI 值，保证结果准确性。
- 第一步：ndvi.updateMask(ndvi.gte(-0.2)) 保留 NDVI≥-0.2 的像元（NDVI 理论范围为 - 1~1，低于 - 0.2 多为水体、裸岩等无植被区域，或传感器噪声，需剔除）。
- 第二步：ndvi.updateMask(ndvi.lte(1)) 保留 NDVI≤1 的像元（超过 1 的数值为异常值，可能由数据误差导致）。
- 第三步：var validSum = nir.add(red).gt(0) 筛选 “近红外 + 红波段” 反射率之和大于 0 的像元，避免出现分母为 0 的计算错误（若两者之和为 0，说明无有效反射信号，属于无效像元）。
- 注：采用多步 updateMask() 而非 and() 逻辑，是为了避免 GEE 中波段类型推断错误（如 “gte 不是函数” 的报错），确保掩膜逻辑稳定执行。
属性继承：return ndvi.copyProperties(img, ['system:time_start']) 将原始影像的时间属性（system:time_start）复制到计算后的 NDVI 影像中，为后续按年度筛选影像提供时间依据。

四、诊断输出（调试用）

print('MODIS全集(前几条)：', modisAll.limit(5));
print('MODIS总影像数：', modisAll.size());

在 GEE 控制台输出关键信息，方便用户验证数据加载和处理结果是否正常：

输出处理后的数据集中前 5 条影像的信息，可查看 NDVI 波段是否存在、属性是否完整。
输出总影像数，若数量为 0，可能是 AOI 与数据集无交集或数据权限问题，需排查。

五、按年度计算均值与导出

YEARS.forEach(function(y) {
  var start = ee.Date.fromYMD(y, 1, 1);
  var end   = start.advance(1, 'year');
  var modisYear = modisAll.filterDate(start, end);
  var ndviMean  = modisYear.mean().clip(aoi);
  
  // 诊断输出年度数据量
  print('Year', y, 'image count:', modisYear.size(), 'mean:', ndviMean);
  
  // 预览首尾两年结果
  if (y === YEAR_START || y === YEAR_END) {
    Map.addLayer(
      ndviMean, 
      {min: 0, max: 1, palette: ['#f7fcf5','#e5f5e0',...,'#00441b']},
      'NDVI_' + y
    );
  }
  
  // 导出到Google Drive
  Export.image.toDrive({
    image: ndviMean,
    description: 'MODIS_NDVI_Mean_' + y,
    folder: EXPORT_FOLDER,
    fileNamePrefix: 'MODIS_NDVI_Mean_' + y,
    region: aoi.geometry(),
    scale: 500,
    crs: 'EPSG:4326',
    maxPixels: 1e13
  });
});

年度循环逻辑：
- YEARS.forEach(function(y) {...}) 遍历 2004-2024 年的每一个年份（y 为当前循环的年份）。
- var start = ee.Date.fromYMD(y, 1, 1) 定义当前年份的起始日期（1 月 1 日）。
- var end = start.advance(1, 'year') 定义当前年份的结束日期（次年 1 月 1 日），通过 advance(1, 'year') 实现年份递进，避免手动计算闰年等问题。
年度 NDVI 均值计算：
- modisAll.filterDate(start, end) 从处理后的 NDVI 数据集中，筛选出当前年份内的所有影像（8 天合成产品，每年约 46 幅影像）。
- modisYear.mean() 对当前年份的所有 NDVI 影像计算像素级均值，得到该年度的 NDVI 均值影像（反映全年植被覆盖平均水平）。
- clip(aoi) 对均值影像进行裁剪，仅保留 AOI 区域内的像元，确保导出数据与研究区域完全匹配。
结果预览（可选）：
- 仅对起始年（2004）和结束年（2024）的均值影像进行地图可视化，避免添加过多图层导致地图卡顿。
- 可视化参数：min: 0, max: 1 设定 NDVI 的显示范围（0~1，植被覆盖从低到高），palette 定义颜色梯度（从浅绿到深绿），直观反映植被覆盖差异。
批量导出配置：Export.image.toDrive(...) 是 GEE 中导出影像到 Google Drive 的核心函数，参数说明如下：
- image: ndviMean：指定导出的影像（当前年度的 NDVI 均值影像）。
- description 和 fileNamePrefix：导出任务的描述和文件前缀，均包含年份标识（如MODIS_NDVI_Mean_2004），方便区分不同年份的文件。
- folder: EXPORT_FOLDER：指定导出到 Google Drive 的目标文件夹（前文定义的GEE_exports）。
- region: aoi.geometry()：指定导出的地理范围（AOI 的几何对象）。
- scale: 500：设定导出影像的空间分辨率（500 米，与 MOD09A1 原始分辨率一致，避免数据重采样导致的精度损失）。
- crs: 'EPSG:4326'：设定坐标参考系（WGS84 经纬度坐标系，全球通用，方便后续 GIS 软件处理）。
- maxPixels: 1e13：设定最大导出像元数（1e13 为极大值，避免因 AOI 过大导致导出失败）。

六、任务完成提示

print('✅ 已创建 2004–2024 年度 NDVI 导出任务，请到 Tasks 面板运行。');

在 GEE 控制台输出明确的完成提示，告知用户批量导出任务已创建，需手动到 GEE 界面的 “Tasks” 面板中启动所有导出任务（GEE 默认不会自动运行导出任务，需用户确认）。

七、关键技术要点与注意事项

技术要点：

数据选择合理性：选择 MOD09A1 数据集（8 天合成反射率），兼顾时间分辨率（8 天更新）和空间分辨率（500 米），适合年度尺度的植被覆盖分析。
掩膜逻辑严谨性：通过 “NDVI 范围筛选 + 分母非零筛选”，剔除无效像元，避免异常值影响年度均值结果，是代码准确性的核心保障。
批量处理高效性：利用ee.List.sequence()生成年份列表，结合forEach()循环实现 21 年数据的自动化处理，无需手动逐年度操作，大幅提升工作效率。
属性继承必要性：复制原始影像的system:time_start属性，确保后续filterDate()能准确筛选出各年度的影像，是按年度聚合的前提。

注意事项：

权限与 AOI 问题：若var aoi = table报错，需替换为公共 AOI（如var aoi = ee.Geometry.Rectangle([x1,y1,x2,y2])，通过经纬度范围定义），或确保对table数据集有访问权限。
数据量与导出时间：21 年的 NDVI 均值影像导出任务可能需要较长时间（取决于 AOI 大小和 GEE 服务器负载），需耐心等待，避免重复提交任务。
参数可调整性：
- 时间范围：可修改YEAR_START和YEAR_END，适配其他时间段（如 1999-2020，需确保 MOD09A1 数据集有对应年份数据）。
- 空间分辨率：scale参数可调整（如 1000 米，加快导出速度；但不可低于 500 米，否则超出原始数据精度）。
- 坐标参考系：crs可改为其他坐标系（如 UTM 投影），需根据后续数据使用场景调整。
调试与验证：建议先运行代码查看控制台输出的 “总影像数”“各年度影像数”，确认数据加载正常后，再启动导出任务；若某年度影像数为 0，可能是该年度该区域无 MODIS 数据，需排查时间范围或 AOI。