这里记录一下使用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);
在应用的过程中,发现变量的重要性总是在变,不知道是不是正常的。。不过也总能找到一两个在任何尝试中重要性都比较低的变量,所以去掉了。
后记
大佬真的太多了!虽然在尝试过程中闹了不少乌龙,但是最后总算能跑出来一个还挺不错的结果。我尝试了不同的合成方式、不同的特征变量、影像集不同时间段,期间还改了几次标签,感觉结果有点玄学。。标签真的太重要,精度很低的时候我发现好几个标签都打错了。。。更改之后有明显提升。但是对于长时间年份,每年都打几百上千的标签有点要命,所以打算把建筑、沙地、水体这种不会有明显变化的标签复用。还有就是GEE因为波段太多会停摆,在想如果在GEE上下载处理好的影像(mean)到本地,再用Python去跑随机森林代码效果会更好吗?苦于Python水平也很差,所以就先这样吧。。。