基于GEE平台的Landsat8、Sentinel2、MODIS“去云”处理及FUI水色指数运算

文章目录

基于GEE平台的Landsat8、Sentinel2、MODIS“去云”处理及FUI水色指数运算
数据集
一、“去云”处理
- 1. Landsat 8
- 2. Sentinel-2
- 3. MODIS
二、FUI水色指数计算（以Landsat 8为例）
参考文章

数据集

Landsat 8 ：USGS Landsat 8 Surface Reflectance Tier 1 图集详情

Sentinel-2 ：Sentinel-2 MSI: MultiSpectral Instrument, Level-2A图集详情

MODIS ：MOD09A1.006 Terra Surface Reflectance 8-Day Global 5图集详情

一、“去云”处理

1. Landsat 8


//“ROI”为自定义兴趣区
Map.centerObject(ROI);// 时间序列,2018 TO 2020
// for(var i = 2018;i<=2020;i++){//   var yearn = i;
//   get_yearly_FUI(yearn);
// }//设置某一年份、月份
get(2019,10);function get(year,month){//设置需要提取的日期var startDate = ee.Date.fromYMD(year, month, 1);var endDate = ee.Date.fromYMD(year, month,30);//去云函数function maskL8sr(image) {// Bits 3 and 5 are cloud shadow and cloud, respectively.var cloudShadowBitMask = 1 << 3;var cloudsBitMask = 1 << 5;// Get the pixel QA band.avar qa = image.select('pixel_qa');// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).and(qa.bitwiseAnd(cloudsBitMask).eq(0));return image.updateMask(mask);}//筛选数据：数据集、研究范围、研究时间段、去云操作var bjL8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR").filterBounds(ROI)    .filterDate(startDate,endDate)//.filter(ee.Filter.lt('CLOUD_COVER', 20)).map(maskL8sr);//求得均值，裁剪出兴趣区                 var mean = bjL8.mean().clip(ROI);
//配色  var viz = {bands: ['B4', 'B3', 'B2'],min: 0,max:3000,gamma: 1.4,};//可视化,添加mean图层Map.addLayer(mean,viz,"BJ_RGB");
}

2. Sentinel-2

Map.centerObject(ROI);// 时间序列,2018 TO 2020
// for(var i = 2018;i<=2020;i++){//   var yearn = i;
//   get_yearly_FUI(yearn);
// }//设置某一年份、月份
get(2019,10);
function get(year,month){//设置需要提取的日期var startDate = ee.Date.fromYMD(year, month, 1);var endDate = ee.Date.fromYMD(year, month,30);//去云函数function maskS2(image) {var qa = image.select('QA60');// Bits 10 and 11 are clouds and cirrus, respectively.var cloudBitMask = 1 << 10;var cirrusBitMask = 1 << 11;// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));return image.updateMask(mask).divide(10000);}var bjS2 = ee.ImageCollection("COPERNICUS/S2_SR").filterBounds(ROI).filterDate(startDate,endDate)//.filter(ee.Filter.lt('CLOUD_COVER', 20)).map(maskS2);var mean = bjS2.mean().clip(ROI); var viz = {bands: ['B4', 'B3', 'B2'],min: 0,max: 0.3,};Map.addLayer(mean,viz,"BJ_RGB");
}

3. MODIS

Map.centerObject(ROI);
// 时间序列,2018 TO 2020
// for(var i = 2018;i<=2020;i++){//   var yearn = i;
//   get_yearly_FUI(yearn);
// }get_ym_FUI(2019,1);function get_ym_FUI(year,month){var startDate = ee.Date.fromYMD(year, month, 1);var endDate = ee.Date.fromYMD(year, month,30);function bitwiseExtract(value, fromBit, toBit){if (toBit === undefined) toBit = fromBitvar maskSize = ee.Number(1).add(toBit).subtract(fromBit)var mask = ee.Number(1).leftShift(maskSize).subtract(1)return value.rightShift(fromBit).bitwiseAnd(mask)}function maskMO(image){var qa = image.select('StateQA');var cloudState = bitwiseExtract(qa, 0, 1) ;var cloudShadowState = bitwiseExtract(qa, 2);var cirrusState = bitwiseExtract(qa, 8, 9);var mask = cloudState.eq(0) // Clear.and(cloudShadowState.eq(0)) // No cloud shadow.and(cirrusState.eq(0)); // No cirrusreturn image.updateMask(mask);} var bjMO = ee.ImageCollection("MODIS/006/MOD09A1").filterBounds(ROI).filterDate(startDate,endDate).map(maskMO);       var mean = bjMO.mean().clip(ROI);Map.addLayer(mean,{min: 0, max: 3500, bands: ['sur_refl_b01', 'sur_refl_b04', 'sur_refl_b03']},'img_m')}

二、FUI水色指数计算（以Landsat 8为例）

Map.centerObject(ROI);// 时间序列,2018 TO 2020
// for(var i = 2018;i<=2020;i++){//   var yearn = i;
//   get_yearly_FUI(yearn);
// }//设置某一年份、月份
get_ym_FUI(2021,1);function get_ym_FUI(year,month){//设置需要提取的日期var startDate = ee.Date.fromYMD(year, month, 1);var endDate = ee.Date.fromYMD(year, month,31);  //记得大小月！！！//筛选 Landsat 8 数据function maskL8sr(image) {// Bits 3 and 5 are cloud shadow and cloud, respectively.var cloudShadowBitMask = 1 << 3;var cloudsBitMask = 1 << 5;// Get the pixel QA band.avar qa = image.select('pixel_qa');// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).and(qa.bitwiseAnd(cloudsBitMask).eq(0));return image.updateMask(mask);}var bjL8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR").filterBounds(ROI).filterDate(startDate,endDate)//.filter(ee.Filter.lt('CLOUD_COVER', 20)).map(maskL8sr);var mean = bjL8.mean().clip(ROI);//NDWI指数提取水体，用掩膜对已筛选出的影像进行水体提取var NDWI = mean.normalizedDifference(['B3', 'B5']).rename('ndWi');var water_mask = NDWI.gte(0);var img_clip = mean.updateMask(water_mask);Map.addLayer(img_clip, {bands: ['B5', 'B4', 'B3'], min: 0, max: 3000, gamma: 1.4}, 'img_clip', true);             print('example img info: ', img_clip); //计算三刺激值var getFU = function(image){var X = image.expression('11.053 * b1 + 6.950 * b2 + 51.135 * b3 + 34.457 * b4',{b1:image.select('B1'),b2: image.select('B2'),b3: image.select('B3'),b4: image.select('B4')}).rename('X');var Y = image.expression('1.32 * b1 + 21.053 * b2 + 66.023 * b3 + 18.034 * b4',{b1:image.select('B1'),b2: image.select('B2'),b3: image.select('B3'),b4: image.select('B4')}).rename('Y');var Z = image.expression('58.038 * b1 + 34.931 * b2 + 2.606 * b3 + 0.016 * b4',{b1: image.select('B1'),b2: image.select('B2'),b3: image.select('B3'),b4: image.select('B4')}).rename('Z');//计算色度角var chrx = X.divide((X.add(Y).add(Z))).rename('chrx');var chry = Y.divide((X.add(Y).add(Z))).rename('chry');var atan2 = chrx.subtract(0.3333).atan2((chry).subtract(0.3333)).rename('atan2');var hue_pre = atan2.multiply(ee.Number(180)).divide(Math.PI).rename('hue_pre');var hue = hue_pre < 0 ? (hue_pre + 360) : hue_pre.rename('hue');var a = hue.divide(100).rename('a');//修正色度角var cor = a.pow(5).multiply(-52.16).add(a.pow(4).multiply(373.81)).add(a.pow(3).multiply(-981.3)).add(a.pow(2).multiply(1134.19)).add(a.multiply(-533.61)).add(76.72).rename('cor');var hueAngle = hue.add(cor).rename('hueAngle');//FUI值——21个等级配色 var FU = a.where(hueAngle.gt(19.000).and(hueAngle.lte(22.741)),21).where(hueAngle.gt(22.741).and(hueAngle.lte(26.337)),20).where(hueAngle.gt(26.337).and(hueAngle.lte(30.439)),19).where(hueAngle.gt(30.439).and(hueAngle.lte(34.906)),18).where(hueAngle.gt(34.906).and(hueAngle.lte(39.769)),17).where(hueAngle.gt(39.769).and(hueAngle.lte(45.129)),16).where(hueAngle.gt(45.129).and(hueAngle.lte(50.665)),15).where(hueAngle.gt(50.665).and(hueAngle.lte(56.435)),14).where(hueAngle.gt(56.435).and(hueAngle.lte(62.186)),13).where(hueAngle.gt(62.186).and(hueAngle.lte(67.957)),12).where(hueAngle.gt(67.957).and(hueAngle.lte(74.572)),11).where(hueAngle.gt(74.572).and(hueAngle.lte(83.346)),10).where(hueAngle.gt(83.346).and(hueAngle.lte(94.037)),9).where(hueAngle.gt(94.037).and(hueAngle.lte(109.054)),8).where(hueAngle.gt(109.054).and(hueAngle.lte(132.999)),7).where(hueAngle.gt(132.999).and(hueAngle.lte(163.084)),6).where(hueAngle.gt(163.084).and(hueAngle.lte(190.779)),5).where(hueAngle.gt(190.779).and(hueAngle.lte(209.994)),4).where(hueAngle.gt(209.994).and(hueAngle.lte(220.977)),3).where(hueAngle.gt(220.977).and(hueAngle.lte(227.168)),2).where(hueAngle.gt(227.168),1).rename('FU');return image.addBands(FU);}// 将计算的FUI可视化，添加到图层var visFU = {min: 1,max: 21,palette: ['2158BC', '316DC5', '327CBB', '4B80A0', '568F96', '6D9298','698C86', '759E72', '7BA654', '7DAE38', '95B645', '94B660','A5BC75', 'AAB86D', 'ADB55F', 'A8A965', 'AE9F5C', 'B3A053','AF8A44', 'A4600A', 'A45308']};// 文件名var year_folder = year+''+month;var clipedFU = getFU(img_clip).select('FU').toInt16();Map.addLayer(clipedFU,visFU,year_folder+'FU');// 减少地区，区域参数是Feature几何体var meanDictionary = clipedFU.reduceRegion({reducer: ee.Reducer.mean(),geometry:ROI,scale: 30,maxPixels: 1e13});// 打印结果print(meanDictionary,year_folder);//输出到文件夹：bj，文件名为：year_folder+"_FUI_Landsat"Export.image.toDrive({image:clipedFU,scale:30,crs:'EPSG:32650',maxPixels : 1e13,region:ROI,folder:"bj",description:year_folder+"_FUI_Landsat"});}

参考文章

【1】 ESA SNAP水色反演算法在GEE中的实现
【2】Google Earth Engine的影像“去云”处理（Landsat & MODIS）

Google Earth Engine学习笔记（一）相关推荐

Google Earth Engine 学习7---ASTER GED 地表发射率数据下载
Google Earth Engine 学习7-ASTER GED 地表发射率数据下载 The Advanced Spaceborne Thermal Emission and Reflection ...
GEE(Google Earth Engine)学习——常用筛选器Filter操作
目录一.筛选器Filter (1)关系比较型筛选器 (2)差值筛选器 (3)字符筛选器 (4)详细时间筛选 (5)筛选器叠加筛选二.Join配合Filter进行两数据集联合筛选 (1)仅保留左侧数 ...
GEE(Google Earth Engine) 最基础代码学习笔记二 —— JavaScript 语言
GEE(Google Earth Engine) 学习笔记二 Javascript 语言 1. 注释 print('Hello World!'); 如果要注释,则在代码前面加//,比如: // pri ...
GEE （Google Earth Engine）最基础代码学习笔记三
GEE (Google Earth Engine)代码学习笔记三本次学习核心为:将JavaScript objects and primitives放入Earth Engine 容器传到服务器,并处 ...
GEE （Google Earth Engine）最基础代码学习笔记四——循环
1. for 循环 Google Earth Engine 不推荐使用for循环,官方推荐使用map()实现for循环. 比如: // 产生一个1到10的list var myList = ee.Li ...
Google Earth Engine（GEE）最基础代码学习笔记6——计算坡度坡向
1.Google Earth Engine 计算坡度计算坡度坡向使用ee.Terrain包计算. // 加载 SRTM 影像. var srtm = ee.Image('CGIAR/SRTM90_V ...
全球最大的云计算平台GEE（Google Earth Engine）正确的学习姿势（1）
对于很多人如何再全球使用量最大的云平台上Google Earth Engine进行运算操作,可能到目前都依旧一无所知,不过今天给大家系统的介绍一下GEE平台正确的学习姿势! 首先,如果大家英文水平还可 ...
Google Earth Engine（GEE）对比显示不同城市的地表温度
今天来简单分享下如何在GEE中对比显示不同城市的地表温度.之前分享过如何利用landsat数据去计算热度. 数据介绍: 数据源为MODIS/061/MOD11A2,数据集里面有个波段名字叫做LST_D ...
再谈批量下载Modis数据之Google earth engine
再谈批量下载Modis数据之Google earth engine 最近在学习GEE,就拿Modis数据练个手.本文选用MOD13A1 NDVI产品为例,介绍如何通过gee来批量下载Modis数据到本 ...

Google Earth Engine学习笔记（一）