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'})