GEE：如何进行对MOD09GA数据集进行水体/云掩膜并计算NDVI将其导出至云盘？_gee进行ndvi后掩膜

目录

01 为什么用GEE而不是传统的下载+ENVI+ArcGIS？

02 操作详解

---

01 为什么用GEE而不是传统的下载+ENVI+ArcGIS？

由于地理空间数据云中缺少2015年10月份的NDVI月合成影像，于是查看了地理空间数据云的NDVI数据集处理的一些介绍如下(地理空间数据云 (gscloud.cn))：

本打算去NASA下载的，本来下载链接都已经拿到了，但是一看8个G，这还仅仅只是下载，我还需要进行裁剪拼接、NDVI计算重采样等等操作，比较繁琐，最终还只是得到一个张月合成NDVI影像，这时间和精力不成正比：

 于是打算使用GEE平台进行NDVI的计算并导出。

02 操作详解

使用的数据集：MOD09GA.061 Terra Surface Reflectance Daily Global 1km and 500m

这里为了保证精度，对影像进行了云掩膜和水体掩膜。

首先，定义了起始日期和结束日期：

// 定义日期范围
var start_date = '2015-10-01';
var end_date = '2015-10-31';
'
运行
运行

接着我们定义一个水体和云掩膜函数：

// 定义云和水体掩膜函数
function maskCloudAndWater(image) {
  var QA = image.select('QC_500m');
  // 创建一个空的mask，初始值为1（即所有像素都不被掩膜覆盖）
  var mask = ee.Image.constant(1);
  
  // 遍历每个波段的数据质量标识
  for (var i = 0; i < 2; i++) {  // 因为我选取了两个波段进行ndvi的计算
    // 计算当前波段的数据质量标识的起始位(是从2开始)
    var startBit = 2 + i * 4;
    // 提取当前波段的数据质量标识
    var bandQuality = QA.rightShift(startBit).bitwiseAnd(15);
    // 如果数据质量标识为15，说明该像素可能被云或深海覆盖，需要被掩膜覆盖
    mask = mask.min(bandQuality.neq(15));  // min取两者间小的那个值,逐像元
  }
  
  // 应用掩膜
  return image.updateMask(mask);
}
'
运行
运行

这里的水体掩膜和云掩膜与一般的数据集不太一样，这里的QC波段是针对每一个波段影像都有4位二进制数进行标识，所以是对每一个波段进行mask的求取再将所有mask求一个类似的或运算。

以下是MOD09GA数据集AC_500m波段的一个介绍(来自GEE):

 

这里主要讲讲二进制掩膜的问题。

在我们的函数里，QA是一个32位的整数，也就是说它是由32个二进制数位组成的，每一个波段都有4位的数据质量标识，这4位被放在QA中的某一位置。例如：

我们举一个简单的例子，对于一个某一个像元，它的属性(32位的二进制整数)值为：

0010 1100 1011 0011 0101 1001 0110 1011

(注意：空格是我为了方便阅读加上的)

在上述二进制整数中，每4位表示一个波段的质量部分，那么假定我们关心的是最左侧的0010部分(这就是我们称之为数据质量的部分)。

那么我们可以将其往右边移动(右移操作)28位，得到：

0000 0000 0000 0000 0000 0000 0000 0010

接着我们就可以将其与掩膜值做比较,这里需要使用到位与操作：

位与操作是对两个二进制数进行比较，只有当两个相应的二进制位都为1时，结果的相应位才为1，否则为0。

通过此前的截图我们知道，Bits2~5表示第1个波段的数据质量部分，掩膜值为15(十进制)表示深海或者云层。

而我们知道，十进制的15转化为32位二进制为：

0000 0000 0000 0000 0000 0000 0000 1111

因此二者进行位与操作之后为：

0000 0000 0000 0000 0000 0000 0000 0010

其转化为十进制不等于15，说明该像元位置不是深海或者云层。

---

其他的代码部分由于时间原因就不一一说明，这里贴出完整代码：

// 定义日期范围
var start_date = '2015-10-01';
var end_date = '2015-10-31';
 
// 定义云和水体掩膜函数
function maskCloudAndWater(image) {
  var QA = image.select('QC_500m');
  // 创建一个空的mask，初始值为1（即所有像素都不被掩膜覆盖）
  var mask = ee.Image.constant(1);
  
  // 遍历每个波段的数据质量标识
  for (var i = 0; i < 2; i++) {  // 因为我选取了两个波段进行ndvi的计算
    // 计算当前波段的数据质量标识的起始位(是从2开始)
    var startBit = 2 + i * 4;
    // 提取当前波段的数据质量标识
    var bandQuality = QA.rightShift(startBit).bitwiseAnd(15);
    // 如果数据质量标识为15，说明该像素可能被云或深海覆盖，需要被掩膜覆盖
    mask = mask.min(bandQuality.neq(15));  // min取两者间小的那个值,逐像元
  }
  
  // 应用掩膜
  return image.updateMask(mask);
}
 
// 定义地理空间范围（四川省）
var geom = ee.FeatureCollection('projects/ee-chaoqiezione/assets/china_admin_province')
geom = geom.filter(ee.Filter.eq('省', '四川省'));
// 加载MODIS数据根据日期和地理范围进行筛选
var modis_ndvi = ee.ImageCollection('MODIS/006/MOD09GA')
  .filterDate(start_date, end_date)
  .filterBounds(geom)
  .select(['sur_refl_b02', 'sur_refl_b01', 'QC_500m'])
  .map(function (img) {
    img = maskCloudAndWater(img);  // 水体和云掩膜
    return img.normalizedDifference(['sur_refl_b02', 'sur_refl_b01']).rename('ndvi')  // 计算ndvi
  })
  .mean().clip(geom)
 
// 添加到地图上以便可视化
print(modis_ndvi)  // 命令面板输出简要信息
Map.addLayer(modis_ndvi.select('ndvi'), {min: 0, max: 1, palette: ['blue', 'white', 'green']}, 'MeanNDVI');
Map.centerObject(geom, 6)
 
// 导出至云盘
Export.image.toDrive({
  image: modis_ndvi.select('ndvi'),
  description: 'Mean_NDVI',
  region: geom,
  scale: 500,
  maxPixels: 1e13,
  fileFormat: 'GeoTIFF'})