GEEer成长日记八：Landsat8_SR计算NDVI逐年时序变化，并通过影像判断城市扩张

　　前几期我们挨个介绍了Modis、Landsat、Sentinel-2产品和数据在逐日和逐月时间序列方面的研究。还介绍了Whittaker Smoother在时间序列研究的应用。本期我们将介绍年尺度的时间序列变化，并通过NDVI的影像分析城市扩张。

　　如果对大家有一点点的帮助，记得文末点个赞哦

---

　　话不多说，我们继续搞代码(前几期也没有为大家讲解代码，后续的研究我们会慢慢增加一些注释)：

//还是老样子哈，以广东省2020年为目标
var geometry = ee.FeatureCollection('users/ZhengkunWang/guangdongsheng')
Map.centerObject(geometry,6)
//去云的算法 
var cloudMaskL457 = function(image) {
  var qa = image.select('pixel_qa');
  // If the cloud bit (5) is set and the cloud confidence (7) is high
  // or the cloud shadow bit is set (3), then it's a bad pixel.
  var cloud = qa.bitwiseAnd(1 << 5)
                  .and(qa.bitwiseAnd(1 << 7))
                  .or(qa.bitwiseAnd(1 << 3));
  // Remove edge pixels that don't occur in all bands
  var mask2 = image.mask().reduce(ee.Reducer.min());
  return image.updateMask(cloud.not()).updateMask(mask2);
};
var years = ee.List.sequence(2000, 2011);
var L5_COL = ee.ImageCollection(years
  .map(function(y) {
    var start = ee.Date.fromYMD(y, 1, 1);
    var end = start.advance(12, 'month');
    var ndvi = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
                  .filterDate(start, end)
                  .map(cloudMaskL457)
                  .map(function(image){
                    var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI')
                    return image.addBands(ndvi)
                  })
    return ndvi.reduce(ee.Reducer.median()).clip(roi)
                  .set('Year',y)
}));
print (L5_COL);

var years = ee.List.sequence(2012, 2012);
var L7_COL = ee.ImageCollection(years
  .map(function(y) {
    var start = ee.Date.fromYMD(y, 1, 1);
    var end = start.advance(12, 'month');
    var ndvi = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')
                  .filterDate(start, end)
                  .map(cloudMaskL457)
                  .map(function(image){
                    var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI')
                    return image.addBands(ndvi)
                  })
    return ndvi.reduce(ee.Reducer.median()).clip(roi)
                  .set('Year',y)
}));
print (L7_COL);

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.
  var 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 years = ee.List.sequence(2013, 2020);
var L8_COL = ee.ImageCollection(years
  .map(function(y) {
    var start = ee.Date.fromYMD(y, 1, 1);
    var end = start.advance(12, 'month');
    var ndvi = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
                  .filterDate(start, end)
                  .map(maskL8sr)
                  .map(function(image){
                    var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI')
                    return image.addBands(ndvi)
                  })
    return ndvi.reduce(ee.Reducer.median()).clip(roi)
                  .set('Year',y)
}));
print (L8_COL);

var data = ee.ImageCollection(L5_COL.merge(L7_COL).merge(L8_COL))
print(data)
// Create a time series chart.
var yearlychart = ui.Chart.image.series({
    imageCollection: data.select('NDVI_median'),
    region: roi,
    reducer: ee.Reducer.mean(),
    scale: 500,
    xProperty:'Year'
    }).setOptions({
      interpolateNulls: true,
      lineWidth: 2,
      title: 'NDVI Daily Time Seires',
      vAxis: {title: 'NDVI', viewWindow: {max: 0.8,min: 0.4,}},
      hAxis: {title: 'Date'},
      trendlines: { 0: {title: 'NDVI_trend',type:'linear', showR2: true,  color:'red', visibleInLegend: true}}
    });
print(yearlychart);

　　我们使用2000-2020年Landsat5、7、8的SR影像进行NDVI的反演，当然是在去云之后进行计算，结果使用研究区域内的NDVI_median进行展示。结果如图：

　　在2011-2013之间有一个明显的上升，对于这么大的区域来说，变化算是很大了。然后呢，我们继续通过影像判断一下发生的变化。因为广东实在太大了，而且2000和2020年的影像云量较多，于是我换了个大家都听说过的好地方！

这个区域的NDVI是这样的：

　　看图的话，确实感觉不出它有什么变化，波澜不惊的那种。于是小编只能通过影像进行分析一波。

// 我就看门见山了，直接选取L5和L8在2000和2020的影像
var landsat5 = ee.ImageCollection("LANDSAT/LT05/C01/T1_SR")
                                  .map(function(image){
                                      var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI')
                                      return image.addBands(ndvi);
                                  });
var landsat8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
                                  .map(function(image){
                                    var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI')
                                    return image.addBands(ndvi);
                                  });
// filterBounds() 代表只过滤研究区域内的影像，其他区域就不加载了
// filterDate() 那就是过滤研究时间一目了然
// sort()也就是按照云量排序，选择云量最少的那一幅
// clip()则是将筛选的影像按照研究区域的边界进行裁剪
var image_ls5 = ee.Image(landsat5
  .filterBounds(geometry)
  .filterDate('2000-01-01', '2000-12-31')
  .sort('CLOUD_COVER') 
  .first()
  .clip(geometry));

var image_ls8 = ee.Image(landsat8
  .filterBounds(geometry)
  .filterDate('2020-01-01', '2020-12-31')
  .sort('CLOUD_COVER') 
  .first()
  .clip(geometry));
//接下来是调色
// 显示影像的真实颜色
var trueColor_ls5 = {bands: ['B3', 'B2', 'B1'], min: 0, max: 3000};
Map.addLayer(image_ls5, trueColor_ls5, 'Landsat5 yr2000 true color');

var trueColor_ls8 = {bands: ['B4', 'B3', 'B2'], min: 0, max: 3000};
Map.addLayer(image_ls8, trueColor_ls8, 'Landsat8 yr2020 true color');

// 显示NDVI的假彩色
var vegPalette = ['darkblue', 'blue', 'red', 'orange', 'yellow', 'green', 'darkgreen'];
Map.addLayer(image_ls5.select('NDVI'), {min: -0.4, max: 0.9, palette: vegPalette}, 'Landsat5 yr2000 NDVI');
Map.addLayer(image_ls8.select('NDVI'), {min: -0.4, max: 0.9, palette: vegPalette}, 'Landsat8 yr2020 NDVI');
Map.centerObject(geometry, 11.5)

 2000年的原始影像：

　2020年的原始影像：

 

　　应该是传感器的原因吧，有种焕然一新的感觉，不过这幅图已经能看出城市扩张的影子了。

2000年的NDVI：

2020年的NDVI：

 

　　对比两幅图，我们一起来找不同～今天的日记就到这里。代码已经更新，获取年平均NDVI时间序列的代码请在公众号回复：123001。

　　获取城市扩张判断的代码请在公众号回复123002。

　　如果真的可以帮到你，记得给小编点个赞哦~

更多精彩内容请关注：