谷歌地球引擎计算NDVI、NDWI年平均值

疯狂学习GIS

发布于 2024-12-19 15:35:33

8522

本文介绍在谷歌地球引擎（Google Earth Engine，GEE）中，基于哨兵2号遥感影像数据，计算某一指定区域NDVI、NDWI等指标的年平均值的方法。本文是谷歌地球引擎（Google Earth Engine，GEE）系列教学文章的第25篇。

首先，明确一下本文的需求。我们现在希望，基于哨兵2号遥感影像数据，计算上海市在2024年期间，NDVI与NDWI等2个指数的平均值；相当于最终希望得到2个结果遥感影像，分别表示上述2个指数的平均。

也很显然，通过本文的代码，我们可以实现基于任意遥感影像产品，计算任意时间与空间范围的任意波段指数的任意统计值——当然，这些都是建立在这些遥感影像数据在GEE中有（或者是我们手动将本地的数据上传到GEE中）的前提下。

本文所用代码如下。

var shanghai_coords = ee.Geometry.Rectangle([120.91, 30.67, 122.20, 31.86]);

var s2 = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
    .filterBounds(shanghai_coords)
    .filterDate('2024-01-01', '2024-12-31')
    .map(function(image) {
      var qa = image.select('QA60');
      var cloudBitMask = 1 << 10;
      var cirrusBitMask = 1 << 11;
      var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
                 .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
      return image.updateMask(mask);
    });

function addIndices(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  var ndwi = image.normalizedDifference(['B3', 'B8']).rename('NDWI');
  return image.addBands(ndvi).addBands(ndwi);
}

var s2_with_indices = s2.map(addIndices);

var mean_ndvi = s2_with_indices.select('NDVI').mean().rename('mean_NDVI_2024');
var mean_ndwi = s2_with_indices.select('NDWI').mean().rename('mean_NDWI_2024');

Map.centerObject(shanghai_coords, 10);
Map.addLayer(mean_ndvi, {min: -1, max: 1, palette: ['blue', 'white', 'green']}, 'Mean NDVI 2024');
Map.addLayer(mean_ndwi, {min: -1, max: 1, palette: ['red', 'yellow', 'blue']}, 'Mean NDWI 2024');

Export.image.toDrive({
  image: mean_ndvi,
  description: 'mean_NDVI_2024',
  scale: 10,
  region: shanghai_coords,
  maxPixels: 1e13
});
Export.image.toDrive({
  image: mean_ndwi,
  description: 'mean_NDWI_2024',
  scale: 10,
  region: shanghai_coords,
  maxPixels: 1e13
});

其中，我们首先创建一个表示研究区域的几何对象，用于后续的图像筛选；在这里，我的研究区域是上海市，且为了方便，就直接构建了一个矩形来表示上海地区。

随后，ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')加载Sentinel-2的地表反射率图像集合，并利用.filterBounds(shanghai_coords)获取覆盖上海地区的图像，并基于.filterDate('2024-01-01', '2024-12-31')选择2024年全年内的图像——虽然在写这个代码时，2024年还没有过完，但是代码这么写是没有问题的。

同时，这里用.map(...)进一步对每个图像进行处理——也就是去云。其中，qa = image.select('QA60')表示选择质量评估波段 QA60，cloudBitMask = 1 << 10和cirrusBitMask = 1 << 11定义云和卷云的二进制掩码位；mask = ...创建一个掩码，确保图像中没有云和卷云，并基于image.updateMask(mask)应用掩码更新图像，去除受云和卷云影响的像素。

在这里，之所以选用质量评估波段 QA60的第10位和第11位，是因为在Sentinel-2的地表反射率图像产品中，这两位就是分别表示有无云和卷云的比特位；如下图所示，在GEE的产品介绍中，就会有对其波段含义的详细介绍。如果大家用了其他的遥感影像产品，并且也需要做云掩膜，那就参考自己所用产品的介绍即可。

紧接着，通过addIndices函数为每张图像计算NDVI和NDWI指数。其中，normalizedDifference(['B8', 'B4'])表示使用近红外（B8）和红光（B4）波段计算NDVI，而normalizedDifference(['B3', 'B8'])表示使用绿光（B3）和近红外（B8）波段计算NDWI。随后，addBands表示将新计算的指数作为额外波段添加到原始图像中，s2_with_indices则表示应用addIndices函数到所有图像，生成带有NDVI和NDWI的图像集合。

在这里，具体遥感影像的哪一个波段分别表示哪一个波长，我们也可以在GEE的产品介绍中找到，如下图所示。

随后，通过select('NDVI')和select('NDWI')选择图像中的NDVI与NDWI波段，并计算选定波段的年度平均值；通过rename(...)重命名输出图像，以便更容易识别。

再接下来，Map.centerObject(shanghai_coords, 10)表示将地图中心定位到上海地区，缩放级别为10，并通过Map.addLayer(...)在地图上添加图层，显示NDVI和NDWI的年度平均值。这里的{min: -1, max: 1, palette: [...]}是设置颜色调色板，以可视化不同的指数值。可视化结果如下图所示。

最后，通过Export.image.toDrive(...)将图像导出到Google Drive。其中，image表示要导出的图像，description是导出文件的描述名称，scale则是输出图像的空间分辨率，region表示导出的地理区域，最后的maxPixels则表示允许的最大像素数，确保导出过程不会因过多像素而消耗太多资源（这个值建议大家设置大一些，就用我这个设置即可）。

随后，运行上述代码，可以看到Tasks列表中，已经有了我们的平均值下载任务了；如下图所示。