GEE:如何进行对MOD09GA数据集进行水体/云掩膜并计算NDVI将其导出至云盘?

简介: GEE:如何进行对MOD09GA数据集进行水体/云掩膜并计算NDVI将其导出至云盘?

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'})
目录
相关文章
|
数据挖掘 编译器 Python
忽略NoData计算多时相遥感影像各个像元的平均值:Python
本文介绍基于Python中whitebox模块,对大量长时间序列栅格遥感影像的每一个像元进行忽略NoData值的多时序平均值求取~
151 2
忽略NoData计算多时相遥感影像各个像元的平均值:Python
|
JSON 数据格式
PIE-engine 教程 ——利用NDWI加载青海湖三年水域影像和面积计算
PIE-engine 教程 ——利用NDWI加载青海湖三年水域影像和面积计算
351 0
PIE-engine 教程 ——利用NDWI加载青海湖三年水域影像和面积计算
GEE:如何批量处理并下载指定时间范围的月尺度NDVI数据集(MOD09GA为例)
GEE:如何批量处理并下载指定时间范围的月尺度NDVI数据集(MOD09GA为例)
659 0
|
6月前
|
算法 数据挖掘 定位技术
GEE 案例——如何计算sentinel-2中每一个单景影像的波段的DN值并绘制直方图
GEE 案例——如何计算sentinel-2中每一个单景影像的波段的DN值并绘制直方图
119 5
|
6月前
|
传感器
GEE——使用cart机器学习方法对Landsat影像条带修复以NDVI和NDWI为例(全代码)
GEE——使用cart机器学习方法对Landsat影像条带修复以NDVI和NDWI为例(全代码)
119 0
|
6月前
|
存储 定位技术 Python
Python中ArcPy实现对大量长时间序列栅格遥感影像批量逐像元求取像素平均值
Python中ArcPy实现对大量长时间序列栅格遥感影像批量逐像元求取像素平均值
PIE-engine 教程 ——NDVI和FVC计算时序分析以2013年-2020年福州市为例
PIE-engine 教程 ——NDVI和FVC计算时序分析以2013年-2020年福州市为例
625 0
PIE-engine 教程 ——NDVI和FVC计算时序分析以2013年-2020年福州市为例
Google Earth Engine ——LANDSAT 4系列增强型植被指数EVI指数8天/32天/年际合成数据集
Google Earth Engine ——LANDSAT 4系列增强型植被指数EVI指数8天/32天/年际合成数据集
198 0
Google Earth Engine ——LANDSAT 4系列增强型植被指数EVI指数8天/32天/年际合成数据集
|
大数据 atlas
Google Earth Engine——全球增强植被指数(EVI)产品的基础数据集是MODIS BRDF校正图像(MCD43B4)
Google Earth Engine——全球增强植被指数(EVI)产品的基础数据集是MODIS BRDF校正图像(MCD43B4)
610 0
Google Earth Engine——全球增强植被指数(EVI)产品的基础数据集是MODIS BRDF校正图像(MCD43B4)
|
编解码 区块链
Google Earth Engine——WWF/HydroSHEDS/30DIR该数据集的分辨率为30弧秒。30角秒的数据集是水文条件下的DEM、排水(流)方向和流量累积。1km分辨率DEM
Google Earth Engine——WWF/HydroSHEDS/30DIR该数据集的分辨率为30弧秒。30角秒的数据集是水文条件下的DEM、排水(流)方向和流量累积。1km分辨率DEM
146 0
Google Earth Engine——WWF/HydroSHEDS/30DIR该数据集的分辨率为30弧秒。30角秒的数据集是水文条件下的DEM、排水(流)方向和流量累积。1km分辨率DEM