Google earth engine学习笔记 2: 通过分析提取植被物候的代码学习gee 1
通过一篇 根据 sentinel-2 影像计算植被指数并提取植被物候 的文章所附的代码学习gee,（原作者的github）,写这个东西是为了自己学习，也希望国内的科研圈合作交流开源的氛围更好，所以如果有用到下文中的代码，希望可以引用原作者发表的文章。刚开始学习gee，难免有错误，希望指出。
首先简单介绍代码的主要组成部分，主要使用过创建和打印字符串变量来实现
var textSections = '' +
"\n SECTION 1 - Define parameters and setup " +
"\n SECTION 2 - Define function for cubic interpolation " +
"\n SECTION 3 - Prepare Sentinel-2 data " +
"\n SECTION 4 - Generate composites " +
"\n SECTION 5 - Increase window size for empty composites " +
"\n SECTION 6 - Interpolate time series " +
"\n SECTION 7 - Estimate phenology metrics " +
"\n SECTION 8 - Plot results " +
"\n SECTION 9 - Display legends "
print(textSections)
-
var str_var = ‘str’： 创建字符串变量
-
print(str_var)： 打印
第一部分，定义相关参数，通过earth engine的中的相关函数设置地图，通过创建字符串和数值变量定义参数，通过earth engine中的相关函数创建点和兴趣区这些空间对象
Map.setOptions('satellite')
var geom = /* color: #ff0000 */ee.Geometry.Point([91.38, 30.55]);
var point = geom // inspect results for this location
var roi = point.buffer(100000).bounds() // display results around the point of interest
var vegIndex = 'evi' // Specify the vegetation index (VI): 'ndvi', 'evi', 'gcc', or 'ndpi'.
var th = 0.5 // Define the percentage of amplitude for the estimation of the threshold
var threshMin = 0.2 // minimum NDVI value for the reclassification of snow values
var scale = 50 // scale of the analysis
var year1 = 2019 // year of processing
var lag = 20 // temporal resolution of the composites
var startDate = year1+'-04-15'
var endDate = year1+'-12-20'
Map.centerObject(roi,7)
-
Map.setOptions(‘satellite’) : 设置地图底图为卫星图, 主要参数为地图类型（‘roadmap’,‘terrain’,'satellite’等），地图风格（一个字典）
-
Map.centerObject(): 使对象居中，数字参数为缩放等级
-
var geom = / color: #ff0000 */ee.Geometry.Point([91.38, 30.55])：使用*ee.Geometry.Point()*创建一个空间点对象,主要接收参数为点的坐标 / color: #ff0000 */ 设置点显示的颜色
-
var roi = point.buffer(100000).bounds()：根据点对象创建roi，buffer() 函数根据对象建立缓冲区，主要参数为距离，bounds()为对象添加边界
第二部分，定义插值函数
var cubicInterpolation = function(collection,step){
var listDekads = ee.List.sequence(1, collection.size().subtract(3), 1);
var colInterp = listDekads.map(function(ii){
var ii = ee.Number(ii);
var p0 = ee.Image(collection.toList(10000).get(ee.Number(ii).subtract(1)));
var p1 = ee.Image(collection.toList(10000).get(ii));
var p2 = ee.Image(collection.toList(10000).get(ee.Number(ii).add(1)));
var p3 = ee.Image(collection.toList(10000).get(ee.Number(ii).add(2)));
var diff01 = ee.Date(p1.get('system:time_start')).difference(ee.Date(p0.get('system:time_start')), 'day');
var diff12 = ee.Date(p2.get('system:time_start')).difference(ee.Date(p1.get('system:time_start')), 'day');
var diff23 = ee.Date(p3.get('system:time_start')).difference(ee.Date(p2.get('system:time_start')), 'day');
var diff01nor = diff01.divide(diff12);
var diff12nor = diff12.divide(diff12);
var diff23nor = diff23.divide(diff12);
var f0 = p1;
var f1 = p2;
var f0p = (p2.subtract(p0)).divide(diff01nor.add(diff12nor));
var f1p = (p3.subtract(p1)).divide(diff12nor.add(diff23nor));
var a = (f0.multiply(2)).subtract(f1.multiply(2)).add(f0p).add(f1p);
var b = (f0.multiply(-3)).add(f1.multiply(3)).subtract(f0p.multiply(2)).subtract(f1p);
var c = f0p;
var d = f0;
/
var xValues = ee.List.sequence(0,diff12.subtract(1),step); !!!!!!!!!!!!!!!
var xDates = ee.List.sequence(p1.get('system:time_start'),p2.get('system:time_start'),86400000);
//print(xDates)
var interp = (xValues.map(function(x){
var im = ee.Image(ee.Number(x).divide(diff12));
return (im.pow(3)).multiply(a).add((im.pow(2)).multiply(b)).add(im.multiply(c)).add(d)
.set('system:time_start',ee.Number(xDates.get(x)));
}));
return interp
})
var colInterp = ee.ImageCollection(colInterp.flatten());
return colInterp
}
-
var cubicInterpolation = function(collection,step){ … }： 自定义函数
-
var listDekads = ee.List.sequence(1, collection.size().subtract(3), 1)：使用 ee.List.sequence() 函数创建数列，参数主要为开始，结束，和步长(或者数量)，其中**collection.size()**为影像集的大小
-
var colInterp = listDekads.map(function(ii){ … }: 通过map()实现循环
-
ee.Number(): 创建数字对象
-
collection.toList(10000)： 将数据集转为列表，数字参数为列表长度上限
-
var p0 = ee.Image(collection.toList(10000).get(ee.Number(ii).subtract(1)))： 列表使用get() 为根据下标提取列表中的元素
-
ee.Date(p1.get(‘system:time_start’)).difference(ee.Date(p0.get(‘system:time_start’)), ‘day’)：ee.Date() 创建一个日期对象，p1.get(‘system:time_start’)获取影像的拍摄时间，difference() 计算两个时间对象间的差，可以给定结果的单位（‘year’,'day’等）
-
diff01.divide(diff12)：divide() 除法
-
(f0.multiply(-3)).add(f1.multiply(3)).subtract(f0p.multiply(2)).subtract(f1p)： subtract() 减法，add()加法， multiply() 乘法，pow() 乘方
-
var colInterp = ee.ImageCollection(colInterp.flatten())：flatten() 讲多维数组转为一维，压平
第三部分，遥感影像预处理
// Call Sentinel2 data.
var S2 = ee.ImageCollection("COPERNICUS/S2_SR")
.filterDate(startDate,endDate)
.filterBounds(roi)
.filterMetadata('CLOUD_COVERAGE_ASSESSMENT','less_than',70)
.filterMetadata('SNOW_ICE_PERCENTAGE','less_than',90);
// Mask non-valid observations and reclassify snow values
var S2 = S2.map(function(im){
var SCL = im.select('SCL')
var SCLmask = SCL.eq(1).or(SCL.eq(4)).or(SCL.eq(5)).or(SCL.eq(11)) // mask for non-valid observations
var snowMask = SCL.eq(11) // snow mask
var blue = im.select('B2').multiply(0.0001)
var green = im.select('B3').multiply(0.0001)
var red = im.select('B4').multiply(0.0001)
var nir = im.select('B8').multiply(0.0001)
var swir = im.select('B12').multiply(0.0001)
// Generate VIs
var ndvi = im.normalizedDifference(['B8','B4']).rename('ndvi')
var evi = (ee.Image(2.5).multiply(nir.subtract(red))).divide(nir.add(red.multiply(6)).subtract(blue.multiply(7.5)).add(1)).rename('evi')
var gcc = green.divide(green.add(red).add(blue)).rename('gcc')
var alpha = ee.Image(0.51) // alpha parameter in NDPI formula
var ndpi = (nir.subtract(alpha.multiply(red).add((ee.Image(1).subtract(alpha)).multiply(swir)))).divide(nir.add(alpha.multiply(red).add((ee.Image(1).subtract(alpha)).multiply(swir)))).rename('ndpi')
// select VI of interest
var bio = ndvi.addBands(evi).addBands(ndpi).addBands(gcc)
.select(vegIndex).rename('bio')
return bio.where(bio.lt(threshMin),threshMin) // force low values to threshMin
.where(snowMask,threshMin) // force snow values to threshMin
.updateMask(SCLmask) // mask non-valid observations
.set('system:time_start', im.get('system:time_start'))
});
-
var S2 = ee.ImageCollection(“COPERNICUS/S2_SR”)
.filterDate(startDate,endDate)
.filterBounds(roi)
.filterMetadata(‘CLOUD_COVERAGE_ASSESSMENT’,‘less_than’,70)
.filterMetadata(‘SNOW_ICE_PERCENTAGE’,‘less_than’,90); ： 通过 ee.ImageCollection(“COPERNICUS/S2_SR”) 调用遥感影像集，其中’COPERNICUS/S2_SR’为sentinel-2影像的在数据库中的名称，数据集名称可以通过检索获得，
filterDate() 根据日期范围筛选影像，filterBounds() 设置影像的空间范围，filterMetadata(‘CLOUD_COVERAGE_ASSESSMENT’,‘less_than’,70) 根据云量筛选， filterMetadata(‘SNOW_ICE_PERCENTAGE’,‘less_than’,90) 根据冰雪情况筛选
-
var SCL = im.select(‘SCL’) ：选择影像波段
-
var SCLmask = SCL.eq(1).or(SCL.eq(4)).or(SCL.eq(5)).or(SCL.eq(11))： 根据条件创建掩模图层
-
var ndvi = im.normalizedDifference([‘B8’,‘B4’]).rename(‘ndvi’)：使用 normalizedDifference() 函数计算ndvi，参数为波段的名称，rename() 对波段进行重命名
-
bio = ndvi.addBands(evi).addBands(ndpi).addBands(gcc)： addBands() 添加波段
-
bio.where(bio.lt(threshMin),threshMin)
.where(snowMask,threshMin)
.updateMask(SCLmask)
.set(‘system:time_start’, im.get(‘system:time_start’))：
where()进行条件替换，参数为条件和用来替换的值，updateMask() 更新掩膜， set(‘system:time_start’, im.get(‘system:time_start’)) 设置时间参数
第四部分，构建时间序列
// Create the list of dates of the central day of the window
var listDates = ee.List.sequence(ee.Date(startDate).millis(), ee.Date(endDate).millis(), 86400000*lag)
// map the dates
var colDekadsRaw = ee.ImageCollection(listDates.map(function(dd){
var date_window = ee.Date(ee.Number(dd)) // central day of the moving window
var date_startW = date_window.advance(-lag/2, 'days') // first day of the moving window
var date_endW = date_window.advance(lag/2, 'days') // last day of the moving window
var col_window = S2.filterDate(date_startW,date_endW) // filter the S2 collection for a period equal to lag
var out = col_window.reduce(ee.Reducer.mean()) // compute the average
return out
.set('system:time_start',date_window.millis()) // set time property
.set('empty', col_window.size().eq(0)) // flag dates without images
// .set('nobs',col_window.filterBounds(point).size())
// .set('date_start',date_start)
// .set('date_end',date_end)
// .set('date_window',ee.Date(date_window));
}).flatten())
// filter dates without images and fill them with the minimum NDVI value
var colDekadsEmpty = colDekadsRaw.filterMetadata('empty', 'equals', 1).map(function(im){
return ee.Image(threshMin).rename('bio_mean').copyProperties(im,['system:time_start'])
})
var colDekads = colDekadsRaw.filterMetadata('empty', 'equals', 0).merge(colDekadsEmpty)
// plot composite
var chart1 = ui.Chart.image.series(colDekads.select('bio_mean'), point, ee.Reducer.first(), scale)
.setOptions({title: 'Moving average every '+lag+' days',
lineWidth: 0,
pointSize: 4})
print(chart1)
-
millis()：将日期转化成自1970-01-01T00:00:00Z的毫秒值
-
var date_startW = date_window.advance(-lag/2, ‘days’) ：advance() 时间对象的加减
-
var out = col_window.reduce(ee.Reducer.mean())：reduce(ee.Reducer.mean()) 计算均值, reduce()中可以给各种Reducer来完成各种运算
-
ee.Image(threshMin).rename(‘bio_mean’).copyProperties(im,[‘system:time_start’]): copyProperties() 复制属性
-
var colDekads = colDekadsRaw.filterMetadata(‘empty’, ‘equals’, 0).merge(colDekadsEmpty)： merge() 合并影像集
-
var chart1 = ui.Chart.image.series(colDekads.select(‘bio_mean’), point, ee.Reducer.first(), scale)
.setOptions({title: ‘Moving average every ‘+lag+’ days’,
lineWidth: 0,
pointSize: 4})
print(chart1)：
ui.Chart.image.series() 根据影像集数据画图，主要参数为 所用影像集(波段)，所选空间范围(点，roi), 所选择的 Reducer，用来计算y轴数值，空间尺度(单位为米)， ee.Reducer.first() 选择输入中的第一个数值， setOptions() 图标相关参数的设置，通常为一个字典
