GEE东拼西凑之随机森林分类_gee随机森林

这里记录一下使用landsat5做随机森林分类的代码，理一下思路。很多内容都是到处找教程东拼西凑的，十分感谢各位大佬。

导入研究区、制作标签

首先加载研究区边界，查看需要分类时间的原影像。在影像上添加标签（目视解译）。

点击左边这个像小气球似的地方，修改名称，选择feature，添加properties。我是添加了两个一个是label，是分类名，另一个是lc，也就是landcover，用数字做区分。

 

 

制作完成之后就是这样的，然后需要将这些feature集合为collection。这里我去掉了bareland分类。

var regions = new ee.FeatureCollection([sand,water,crop,wetland,road_canal,construction]);

 值得注意的是，这些标签可以在其他文件中复用。点击imports部分右边蓝色本样子的图标，复制其中的代码到新的file中，会显示convert，转换一下之后就和上面界面一样了。

var regions = new ee.FeatureCollection([sand,water,crop,wetland,road_canal,construction]);
var table = ee.FeatureCollection("projects/ee-lzying1999/assets/Boundary-wgs");
Map.centerObject(table, 9);
Map.addLayer(ee.Image().paint(table, 0, 2), {}, 'roi');
var date = ee.DateRange('2000-05-01','2000-10-01')

制作landsat影像集——尺度转换、添加光谱指数 

 这里是做一个预处理，主要是考虑到GEE中landsat影像集的缩放。以及后续分类要用到的光谱指数，在这里都处理好。在bands那里可以选择不同的波段以确定分类用到的特征变量。不过要注意的是，波段数太多可能超过GEE平台限制（80MiB），所以建议找最合适的变量。

//product:LANDSAT/LT05/C02/T1_L2
function input_landsat5L2(region,date){
//缩放 
  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);
  }
 
//NDVI 
  function NDVI(image) { 
    return image.addBands( 
     image.normalizedDifference(["SR_B4", "SR_B3"]) 
           .rename("NDVI")); 
  } 
 
//EVI
  var EVI = function(image) {
   var evi =  image.expression(
       '2.5 * float(nir - red)/ (nir + 6*red -7.5*blue +1 )',
       {
         'nir': image.select('SR_B4'),
          'red': image.select('SR_B3'),
          'blue':image.select('SR_B1')
       }).rename('EVI');
    return image.addBands(evi);
  };
 
//NDWI 
  function NDWI(image) { 
    return image.addBands( 
      image.normalizedDifference(["SR_B2", "SR_B4"]) 
          .rename("NDWI")); 
  } 
 
//MNDWI 
  function MNDWI(image) { 
    return image.addBands( 
      image.normalizedDifference(["SR_B2", "SR_B5"]) 
          .rename("MNDWI")); 
  } 
  
//NDBI 
  function NDBI(image) { 
    return image.addBands( 
      image.normalizedDifference(["SR_B5", "SR_B4"]) 
          .rename("NDBI")); 
  } 
 
//SI
  var SI = function(image) {
   var si =  image.expression(
       'float(swir + red-nir-blue)/ (swir +red +blue +nir )',
       {
         'swir':image.select('SR_B5'),
         'nir': image.select('SR_B4'),
          'red': image.select('SR_B3'),
          'blue':image.select('SR_B1')
       }).rename('SI');
    return image.addBands(si);
  };
//SAVI
  var SAVI = function(image) {
   var savi =  image.expression(
       '1.5*float(nir-red)/ (red +nir+0.5 )',
       {
         'nir': image.select('SR_B4'),
          'red': image.select('SR_B3'),
       }).rename('SAVI');
    return image.addBands(savi);
  };
//DEM
  function DEM(image){
    var strm=ee.Image("USGS/SRTMGL1_003");
    var dem=ee.Algorithms.Terrain(strm);
    var elevation=dem.select('elevation');
    var slope=dem.select('slope');
    return image.addBands(elevation.rename("ELEVATION")
      ).addBands(slope.rename("SLOPE"))
  };
 
//remove cloud by select bitwise（在网上找的去云方法）
  function bitwiseExtract(value, fromBit, toBit) {
    if (toBit === undefined) toBit = fromBit;
    var maskSize = 1 + toBit - fromBit; //位数
    var mask = ee.Number(1 << maskSize).subtract(1); 
    return value.rightShift(fromBit).bitwiseAnd(mask);
  }
 
  function cloudfree_landsat(image){
    var qa = image.select('QA_PIXEL')
    var cloudState = bitwiseExtract(qa, 3)    //5
    var cloudShadowState = bitwiseExtract(qa, 4)  //3
    var mask = cloudState.eq(0) // Clear
    .and(cloudShadowState.eq(0)) // No cloud shadow
    return image.updateMask(mask)  
  }
//加载landsat影像
  var collection=
ee.ImageCollection("LANDSAT/LT05/C02/T1_L2").filterDate(date).filterBounds(region);
  var preinput = collection.filterMetadata('CLOUD_COVER','less_than',3.0).sort('CLOUD_COVER', true)
              .map(cloudfree_landsat).map(applyScaleFactors) 
              .map(NDVI) 
              .map(NDWI) 
              .map(NDBI)
              .map(MNDWI)
              .map(EVI).map(DEM).map(SI).map(SAVI);
//选择不同的bands作为特征变量，可以多尝试哪个效果更好
  var bands = ["SR_B1", "SR_B2", "SR_B3","SR_B5", "SR_B7", 
  "NDBI",  "NDVI" ,'MNDWI','NDWI',"SAVI",'EVI','ELEVATION'//,'SI','SLOPE',];
 
  var input = preinput.select(bands);  
  print(input);
//看看原影像是什么样的（用于制作标签）  
  Map.addLayer(preinput.mean().clip(region), {min:0, max:0.3, bands: ['SR_B3', 'SR_B2', 'SR_B1']} ,"rawscenecloud");
  return input
} 

划分数据集

这里我比较了mean、median、max、mosaic四种合成/镶嵌方式的效果，最后还是选择了mean。以mean作为分类的底图，也就是用各种波段、指数的平均值进行分类。然后7:3的比例划分数据集。

//调用上面处理好的landsat数据集
var input = input_landsat5L2(table,date)                            
//按不同方式合成
var mean = input.mean().clip(table)
var median = input.median().clip(table)
var max = input.max().clip(table)
var mosaic = input.mosaic().clip(table)
 
// 选取样本，把landcover属性赋予样本
var training = mean.sampleRegions({
  collection: regions,
  properties: ['lc'],
  scale: 30
});
print(training,'training')
 
//样本点随机排列
var withRandom = training.randomColumn('random');
//以7:3的比例划分训练集和测试集
var split = 0.7; 
var trainingPartition = withRandom.filter(ee.Filter.lt('random', split));
var testingPartition = withRandom.filter(ee.Filter.gte('random', split));

确定随机森林树的数目

这里再次感谢各路大佬。

//选取随机森林树数目
 var numTrees = ee.List.sequence(5, 50, 5); 
 var accuracies = numTrees.map(function(t)
 { 
   var classifier1 = ee.Classifier.smileRandomForest(t)
                     .train({
                     features: trainingPartition,
                     classProperty: 'lc',
                     inputProperties: mean.bandNames()
                     });
   return testingPartition
       .classify(classifier1)
       .errorMatrix('lc', 'classification')
       .accuracy();
 }); 
 print(ui.Chart.array.values({
   array: ee.Array(accuracies),
   axis: 0,
   xLabels: numTrees
 }));

然后就可以选择精度最高的了。

训练并分类

这里之前还弄出来笑话了，因为没仔细看大佬的代码，只分类了测试集，没有对整张图像分类，还到处搜为什么不出图。。

//训练
var classifier = ee.Classifier.smileRandomForest(20)
    .train({
      features: trainingPartition, //训练集
      classProperty: 'lc',  //标签
      inputProperties: mean.bandNames() //波段名，设置了输入图像的波段
    });
 
//将分类器应用于整张图    
var img_RF = mean.classify(classifier);
print(img_RF,'img_RF')
var imageVisParam = {min:1,max:6, palette:['yellow','blue','green','darkgreen','brown','grey']}/
Map.addLayer(img_RF,imageVisParam,'Classified_RF');
 
//用分类器对测试集进行分类
var Classified_RF = testingPartition.classify(classifier);
print(Classified_RF,'Classified_RF')

 精度分析

精度分析也是很全面了，论文中常见的系数都能计算。本质都是基于混淆矩阵计算的。

 
//混淆矩阵
var testAccuracy = Classified_RF.errorMatrix('lc', 'classification');
//总准确度
var accuracy = testAccuracy.accuracy();
//用户精度（多分误差）
var userAccuracy = testAccuracy.consumersAccuracy();
//制图精度（漏分误差）
var producersAccuracy = testAccuracy.producersAccuracy();
//Kappa系数，不仅考虑混淆矩阵对角线上精准分类的数量，也考虑了各种漏分、误分情况
var kappa = testAccuracy.kappa();
print('Validation error matrix:', testAccuracy);
print('Validation overall accuracy:', accuracy);
print('User acc:', userAccuracy);
print('Prod acc:', producersAccuracy);
print('Kappa:', kappa);

特征重要性分析

//几乎没有参数需要改
var dict = classifier.explain();
print('Explain:',dict);
 
var variable_importance = ee.Feature(null, ee.Dictionary(dict).get('importance'));
var chart_im =
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_im);

在应用的过程中，发现变量的重要性总是在变，不知道是不是正常的。。不过也总能找到一两个在任何尝试中重要性都比较低的变量，所以去掉了。

 

后记

发现CSDN上大佬真的太多了！虽然在尝试过程中闹了不少乌龙，但是最后总算能跑出来一个还挺不错的结果。我尝试了不同的合成方式、不同的特征变量、影像集不同时间段，期间还改了几次标签，感觉结果有点玄学。。标签真的太重要，精度很低的时候我发现好几个标签都打错了。。。更改之后有明显提升。但是对于长时间年份，每年都打几百上千的标签有点要命，所以打算把建筑、沙地、水体这种不会有明显变化的标签复用。还有就是GEE因为波段太多会停摆，在想如果在GEE上下载处理好的影像（mean）到本地，再用Python去跑随机森林代码效果会更好吗？苦于Python水平也很差，所以就先这样吧。。。