谷歌地球引擎计算DVI、DWI年平均值

谷歌地球引擎计算DVI、DWI年平均值

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

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

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

  本文所用代码如下。

var shanghai_coords = ee.Geometry.Rectangle([120.91, 0.67, 122.20, 1.86]);

var s2 = ee.ImageCollection('COPERICUS/S2_SR_HARMOIZED')
    .filterBounds(shanghai_coords)
    .filterDate('2024-01-01', '2024-12-1')
    .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 = (['B8', 'B4']).rename('DVI');
  var ndwi = (['B', 'B8']).rename('DWI');
  return image.addBands(ndvi).addBands(ndwi);
}

var s2_with_indices = (addIndices);

var mean_ndvi = s2_with_indices.select('DVI').mean().rename('mean_DVI_2024');
var mean_ndwi = s2_with_indices.select('DWI').mean().rename('mean_DWI_2024');

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

Export.({
  image: mean_ndvi,
  description: 'mean_DVI_2024',
  scale: 10,
  region: shanghai_coords,
  maxPixels: 1e1
});
Export.({
  image: mean_ndwi,
  description: 'mean_DWI_2024',
  scale: 10,
  region: shanghai_coords,
  maxPixels: 1e1
});

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

  随后，ee.ImageCollection('COPERICUS/S2_SR_HARMOIZED')加载Sentinel-2的地表反射率图像集合，并利用.filterBounds(shanghai_coords)获取覆盖上海地区的图像，并基于.filterDate('2024-01-01', '2024-12-1')选择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函数为每张图像计算DVI和DWI指数。其中，normalizedDifference(['B8', 'B4'])表示使用近红外（B8）和红光（B4）波段计算DVI，而normalizedDifference(['B', 'B8'])表示使用绿光（B）和近红外（B8）波段计算DWI。随后，addBands表示将新计算的指数作为额外波段添加到原始图像中，s2_with_indices则表示应用addIndices函数到所有图像，生成带有DVI和DWI的图像集合。

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

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

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

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

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

  随后，执行这2个任务，即可开始将结果文件导出至Google Drive，如下图所示。

  导出完毕后，就可以在Google Drive的对应路径下，看到我们刚刚导出到这里的结果图像文件，如下图所示。

  随后，在文件上右键，选择“Download”即可开始下载文件，如下图所示。

  等到文件下载完毕，我们就可以在本地打开这些结果图像文件了。

  至此，大功告成。