【Google Earth Engine】利用GEE进行Landsat 8 SR数据土地利用分类

最近用学习到的知识进行了利用GEE和Landsat 8 SR数据进行土地利用分类的小实验，在这里进行一些学习记录。

一、数据导入

首先在GEE中上传要进行土地利用分类的行政区域边界，这里是以雄安新区为例。

二、遥感数据筛选

使用的数据是Landsat 8 OLI/TIRS传感器的SR数据集，SR数据利用QA波段进行影像去云处理，这里构造了去云函数便于后续调用；

筛选想要进行土地利用分类的时间，并用clip函数将研究区裁剪出来。

// Applies scaling factors.
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}

//L8 cloud_remove
function maskL8sr(image) {
  // 第3位和第5位分别是云影和云。
  var cloudShadowBitMask = 1 << 4;
  var cloudsBitMask = 1 << 3;
 
  // 获取pixel QA band.
  var qa = image.select('QA_PIXEL');
 
  // 明确条件，设置两个值都为0
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));
 
  // 更新掩膜云的波段，最后按照反射率缩放，在选择波段属性，最后赋值给影像
  return image.updateMask(mask)
      .select("SR_B[1-7]*")
      .copyProperties(image, ["system:time_start"]);
}

//Filter image collection for time window, spatial location, and cloud cover
var startDate = ee.Date('2016-01-01'); 
var endDate = ee.Date('2016-12-31'); 
var collection = l8
                  .filterDate(startDate, endDate)
                  .map(applyScaleFactors)//时间过滤
                  .map(maskL8sr)
                  .median()
                  ;
//print(remove_cloud(collection))
var image=collection.clip(roi);
var visualParam = {bands: ['SR_B4', 'SR_B3', 'SR_B2'], min:0.0,max: 0.3};//可视化参数1
print("l8Image",visualParam,image)
Map.centerObject(roi, 11); 
Map.addLayer(image, visualParam,"l8Image");

三、 分类特征集构建

可以构建NDVI/NDBI/MNDWI作为光谱特征；利用DEM数据的坡度和高程作为地形特征；然后将构建的这些特征作为影像的一个波段进行使用。

var mndwi = image.normalizedDifference(['SR_B3', 'SR_B6']).rename('MNDWI');//计算MNDWI
var ndbi = image.normalizedDifference(['SR_B6', 'SR_B5']).rename('NDBI');//计算NDBI
var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');//计算NDVI
var strm=ee.Image("USGS/SRTMGL1_003");
var dem=ee.Algorithms.Terrain(strm);
var elevation=dem.select('elevation');
var slope=dem.select('slope');

image=image
      .addBands(ndvi)
      .addBands(ndbi)
      .addBands(mndwi)
      .addBands(elevation.rename("ELEVATION"))
      .addBands(slope.rename("SLOPE"))

四、 样本的选取

选取样本的方法就是在视图窗口用point或polygon进行分类选取，这里要注意在样本的属性设置里要改变图层类型为featurecollection，并添加分类属性landcover与这个类对应的值。

** 五、样本的处理**

这里以选取的透水面和不透水面两类样本进行示例，将两类的全部样本进行混合，然后随机排列，将样本以7：3的比例分为训练集和测试集。

var classNames = butoushuimian2014.merge(toushuimian2014);

var bands = ['SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7','MNDWI','NDBI','NDVI','SLOPE', 'ELEVATION'];

var training = image.select(bands).sampleRegions({
  collection: classNames,
  properties: ['landcover'],
  scale: 30
});

// random uniforms to the training dataset.
var withRandom = training.randomColumn('random');
var split = 0.7; 
var trainingPartition = withRandom.filter(ee.Filter.lt('random', split));
var testingPartition = withRandom.filter(ee.Filter.gte('random', split));

六、选取合适的分类方法进行分类

下面就可以分类了，除了这里使用的smilerandomForest()随机森林法，还可以使用libsvm/smileCart() 等方法，具体可以查询GEE的函数文档。

var classProperty = 'landcover';
var classifier = ee.Classifier.smileRandomForest(30).train({
  features: trainingPartition,
  classProperty: 'landcover',
  inputProperties: bands
});
var classified = image.select(bands).classify(classifier);
Map.addLayer(classified,imageVisParam3);
print(classified)
var test = testingPartition.classify(classifier);

此外，在使用随机森林方法进行分类时，可以根据自己的需要来进行最佳决策树数量的选择和所使用特征的特征重要性计算。

//决策树数量选择
var numTrees = ee.List.sequence(5, 100, 5); 
var accuracies = numTrees.map(function(t)
{ 
  var classifier = ee.Classifier.smileRandomForest(t)
                    .train({
                features: trainingPartition,
  classProperty: 'landcover',
  inputProperties: bands
                    });
  return testingPartition
      .classify(classifier)
      .errorMatrix('landcover', 'classification')
      .accuracy();
}); 
print(ui.Chart.array.values({
  array: ee.Array(accuracies),
  axis: 0,
  xLabels: numTrees
}));
//随机森林特征重要性
var dict = classifier.explain();
print('Explain:',dict);

var variable_importance = ee.Feature(null, ee.Dictionary(dict).get('importance'));
var chart =
ui.Chart.feature.byProperty(variable_importance)
.setChartType('ColumnChart')
.setOptions({
title: 'Random Forest Variable Importance',
legend: {position: 'none'},
hAxis: {title: 'Bands'},
vAxis: {title: 'Importance'}
});

print(chart);

七、结果评价

使用混淆矩阵、OA、KAPPA系数进行分类效果的评价。

var confusionMatrix = test.errorMatrix('landcover', 'classification');
print('confusionMatrix',confusionMatrix);
print('overall accuracy', confusionMatrix.accuracy());
print('kappa accuracy', confusionMatrix.kappa());

八、各分类类型的面积计算

土地利用分类完成后通常要进行各种地类的面积计算，我这里使用的是pixelArea方式，此外还可以使用count方式。

var areaall = ee.Image.pixelArea().addBands(classified)
.reduceRegion({
        reducer: ee.Reducer.sum().group({
          groupField: 1, 
          groupName: "landcover"
        }),
        geometry: roi, 
        scale: 30,
        maxPixels:10e15,
})
print(areaall)

九、 结果输出

这里我是把结果导出到了Drive里，在Google Drive里下载下来就可以在Arcgis和ENVI里进行使用啦~

Export.image.toDrive({
        image:  classified,
        description: 'xiongan_RF',
        folder: 'xiongan_RF',
        scale: 30,
        region: roi,
        maxPixels:34e10
      });

展示一下用这样的方法得到的分类结果，效果还是可以的（大前提是样本的质和量都要够哈）~