最近一直在做土地利用分类,尝试了非监督分类、随机森林,感觉制作标签好费事,就想通过光谱指数阈值偷懒。当然结果表明,这样分割只能看大致趋势!!实际分类结果不准!!对于一些趋势分析倒是勉强可以用用,而且确实省时省力。
选取典型土地利用类型,查看光谱特征
首先我在地图上选取了几种典型的土地利用类型,计算他们的NDVI、MNDWI、SI、NDBI等光谱指数,查看哪种光谱指数能分类得更好。结果发现,EVI和SAVI和NDVI结果相差不大,MNDWI效果没有NDWI效果好,IBI结果与NDBI效果差不多,随后就选择了NDVI、MNDWI、SI、NDBI这四种指标作为分割标准。
//将各土地利用类型合成一个featurecollection;table是导入的研究区边界
var regions = new ee.FeatureCollection([sand,water,crop,wetlands,city,road,barelands]);
Map.centerObject(table, 9);
Map.addLayer(ee.Image().paint(table, 0, 2), {}, 'roi');
var date = ee.DateRange('2022-06-01','2022-09-01')
//landsat影像处理与光谱指数计算
//product:LANDSAT/LT05/C02/T1_L2
function input_landsat8L2(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_B5", "SR_B4"])
.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_B5'),
'red': image.select('SR_B4'),
'blue':image.select('SR_B2')
}).rename('EVI');
return image.addBands(evi);
};
//NDWI
function NDWI(image) {
return image.addBands(
image.normalizedDifference(["SR_B3", "SR_B5"])
.rename("NDWI"));
}
//MNDWI
function MNDWI(image) {
return image.addBands(
image.normalizedDifference(["SR_B3", "SR_B6"])
.rename("MNDWI"));
}
//NDBI
function NDBI(image) {
return image.addBands(
image.normalizedDifference(["SR_B6", "SR_B5"])
.rename("NDBI"));
}
//SI
var SI = function(image) {
var si = image.expression(
'float(swir + red-nir-blue)/ (swir +red +blue +nir )',
{
'swir':image.select('SR_B6'),
'nir': image.select('SR_B5'),
'red': image.select('SR_B4'),
'blue':image.select('SR_B2')
}).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_B5'),
'red': image.select('SR_B4'),
}).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"))
}
//IBI
function IBI(image){
var ibiA = image.expression('2 * SWIR1 / (SWIR1 + NIR)', {
'SWIR1': image.select('SR_B6'),
'NIR' : image.select('SR_B5')
}).rename(['IBI_A']);
var ibiB = image.expression('(NIR / (NIR + RED)) + (GREEN / (GREEN + SWIR1))', {
'NIR' : image.select('SR_B5'),
'RED' : image.select('SR_B4'),
'GREEN': image.select('SR_B3'),
'SWIR1': image.select('SR_B6')
}).rename(['IBI_B']);
var ibiAB = ibiA.addBands(ibiB);
var ibi = ibiAB.normalizedDifference(['IBI_A', 'IBI_B']);
return image.addBands(ibi.rename("IBI"))
}
//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)
}
//load landsat
var collection = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").filterDate(date).filterBounds(region);//LANDSAT/LT05/C02/T1_L2
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).map(IBI)//.map(nightlight).qualityMosaic('NDVI').clip(region) //.mosaic() can be changed in different ways
var bands = [
"SR_B6", "SR_B2", "SR_B3","SR_B4","SR_B5", "SR_B7",
"NDBI", "NDVI" ,'MNDWI','NDWI',"SAVI",'EVI','ELEVATION','SI','SLOPE','IBI'
];
//get the base map
var input = preinput.select(bands);
print(input)
Map.addLayer(preinput.mean().clip(region), {min:0, max:0.3, bands: ['SR_B4', 'SR_B3', 'SR_B2']} ,"rawscenecloud");
return input
}
var input = input_landsat8L2(table,date)
//均值合成
var mean = input.mean().clip(table)
var ndvi_crop = mean.select('NDVI').reduceRegion({
reducer: ee.Reducer.mean(),
geometry: crop.geometry(),
scale: 30
});
print('ndvi_mean_crop:',ndvi_crop)
var ndvi_bare = mean.select('NDVI').reduceRegion({
reducer: ee.Reducer.mean(),
geometry: barelands.geometry(),
scale: 30
});
print('ndvi_mean_bareland:',ndvi_bare)
var ndvi_wetland = mean.select('NDVI').reduceRegion({
reducer: ee.Reducer.mean(),
geometry: wetlands.geometry(),
scale: 30
});
print('ndvi_mean_wetland:',ndvi_wetland)
var ndvi_sand = mean.select('SI').reduceRegion({
reducer: ee.Reducer.mean(),
geometry: sand.geometry(),
scale: 30
});
print('ndvi_mean_sand:',ndvi_sand)
由于影像质量每年都不确定,因此先根据每年平均值大致判断阈值(这里我为了省事直接人为设定阈值了,所以结果受主观影响很大,最好设计一套更符合每年情况的阈值确定方法)。
通过阈值进行分类
这里涉及的总思路是:(1)逐年计算6-9月样本点(田地、湿地、裸地)平均NDVI,作为阈值参考;(2)通过MNDWI分割水体(由于MNDWI对于水体分割效果很好,设置固定阈值),通过设置NDVI阈值划分湿地(阈值逐年根据NDVI均值变化);(3)通过设置NDVI阈值进一步划分田地(阈值逐年根据NDVI均值变化,考虑到夏季水面植物生长,用湿地、水域mask作为约束),通过SI阈值划分沙地(阈值不变,但用田地mask作为约束);(4)通过NDVI阈值划分裸地(阈值不变,使用水域、沙地、田地作为约束);(5)通过NDBI值划分不透水面(阈值不变,使用裸地、沙地、水域作为约束)。在各类型相互约束前提下,通过目视大致判断土地利用分布,调整阈值(湿地、田地),将各个土地类型分别赋值,不属于此类的像素mask,然后合成一张分类图,再按照灌域进行像素统计。这里还是强调,采用的方法太依赖人工经验,之后有机会的话再认真考虑阈值的选取。
//湿地——NDVI
var wetland = mean.select('NDVI').gte(0.73)
var wetland_mask = wetland.mask(wetland)
var WETLAND = wetland_mask.where(wetland_mask.select('NDVI').eq(1),1).rename('lc').set(ee.Dictionary(['lc',1])).toInt8()
var otherwl = mean.select('NDVI').lt(0.73)
//水体——MNDWI
var waterbodies = mean.select('MNDWI').gte(0.2)
var waterbodies_mask = waterbodies.mask(waterbodies)
var WATER = waterbodies_mask.where(waterbodies_mask.select('MNDWI').eq(1),3).rename('lc').set(ee.Dictionary(['lc',3])).toInt8()
var otherw = mean.select('MNDWI').lt(0.2)
//农田——NDVI
var field = mean.select('NDVI').gte(0.4)//0.35
var field_mask = field.mask(field).updateMask(otherwl.mask(otherwl)).updateMask(otherw.mask(otherw))
var FIELD = field_mask.where(field_mask.select('NDVI').eq(1),2).rename('lc').set(ee.Dictionary(['lc',2])).toInt8()
var otherf = mean.select('NDVI').lt(0.4)
//沙地——SI
var sands = mean.select('SI').gte(0.13)
var sands_mask = sands.mask(sands).updateMask(otherw.mask(otherw)).updateMask(otherf.mask(otherf))
var SANDS = sands_mask.where(sands_mask.select('SI').eq(1),4).rename('lc').set(ee.Dictionary(['lc',4])).toInt8()
var others = mean.select('SI').lt(0.13)
//var others = mean.select('SI').mask(mean.select('SI').lt(0.12)).unmask(0);
//裸地——NDVI
var bareland = mean.select('NDVI').gte(0.2)
var bareland_mask = bareland.mask(bareland).updateMask(otherf.mask(otherf))
.updateMask(otherw.mask(otherw)).updateMask(others.mask(others))
var BARELAND = bareland_mask.where(bareland_mask.select('NDVI').eq(1),5).rename('lc').set(ee.Dictionary(['lc',5])).toInt8()
var otherb = mean.select('NDVI').lt(0.2)
//不透水面——NDBI
var impervious = mean.select('NDBI').gt(-0.1)
var impervious_mask = impervious.mask(impervious).updateMask(others.mask(others))
.updateMask(otherb.mask(otherb)).updateMask(otherw.mask(otherw))
var IMPERVIOUS = impervious_mask.where(impervious_mask.select('NDBI').eq(1),6).rename('lc').set(ee.Dictionary(['lc',6])).toInt8();
//加载各分类情况分布
Map.addLayer(wetland_mask,{palette:["006400"]},'wetland')
Map.addLayer(field_mask,{palette:["7FFF00"]},'field')
Map.addLayer(bareland_mask,{palette:["8B4513"]},'bareland')
Map.addLayer(impervious_mask,{palette:["8A8A8A"]},'impervious')
Map.addLayer(sands_mask,{palette:["FFF68F"]},'sands')
Map.addLayer(waterbodies_mask,{palette:["1E90FF"]},'waterbodies')
这里用到了mask,mask会将符合条件的像素赋值为1,不符合条件的像素赋值为0,然后将值为0的像素掩膜掉,这里要注意不是删除原像素,只是将其掩膜。updateMask可以应用其他的mask,还有一个unmask()的命令,表示将掩膜掉的像素赋值为括号内的数字,如-9999等。
在上面的阈值分割过程中,使用了多种光谱指数对6种土地利用类型进行分割,需要注意的是不同土地利用类型之间不能有重叠,因此除了强特征的土地类型(湿地、沙地、水体),余下的(农田、裸地、不透水面)在分割的过程中需要对已分好的类型进行掩膜处理。这一步很关键,不然会影响后面面积统计的结果。类型之间的相互约束真的很重要。
在分类结束后,通过where()函数将掩膜出来的像素赋予不同的数字作为土地类型的区分,并将波段重命名为lc(landcover)。set函数是将每个影像添加一个property,名字为lc, 这是之前失败的尝试,实际上可加可不加。
合并各掩膜,进行面积统计
各土地类型分类完成后,将其集合为imagecollection,此时每张影像的波段都是lc。集合后,使用mean()函数将各个影像合成为一张土地分类图。如果之前土地掩膜没有弄好的话,这里合成就会出大问题。之前跳错说type不统一,因此在上面rename后将其转换为int(toInt8()),之后就不报错了。
//合成为一张土地分类图
var classification = ee.ImageCollection.fromImages([WETLAND,FIELD,WATER,SANDS,BARELAND,IMPERVIOUS])
//print(classification2)
Map.addLayer(classification2.mean(),{min:1,max:6,palette:["006400","7FFF00","1E90FF","FFF68F","8B4513","8A8A8A"]},'mosaic')
//提取研究区的小区域名称
var point = table.aggregate_array('Name1')
//for循环统计各小区域的土地类型面积
for(var i = 0; i<6 ; i++) {
var roi = table.filterMetadata('Name1','equals',point.get(i));
//count area
var areaImage = ee.Image.pixelArea().addBands(classification.mean())
var areas = areaImage.reduceRegion({
reducer: ee.Reducer.sum().group({
groupField: 1,
groupName: 'type',
}),
geometry: roi,
scale: 30,
maxPixels: 1e13
});
var classAreas = ee.List(areas.get('groups'))
var classAreaLists = classAreas.map(function(item) {
var areaDict = ee.Dictionary(item)
var classNumber = ee.Number(areaDict.get('type')).format()
var area = ee.Number(
areaDict.get('sum')).divide(1e6)//.round()
return ee.List([classNumber, area])
})
var result1 = ee.Dictionary(classAreaLists.flatten())
print('area:'+ i ,result1)
//}
}
//统计整个区域的各土地利用类型面积
var image_area = ee.Image.pixelArea().addBands(classification.mean())
var Area = image_area.reduceRegion({
'reducer': ee.Reducer.sum().group({
'groupField': 1,
'groupName': 'class'
}),
'geometry': table,
'scale': 30,
'maxPixels': 1e13
})
var classAreas1 = ee.List(Area.get('groups'))
var classAreaLists1 = classAreas1.map(function(item) {
var areaDict1 = ee.Dictionary(item)
var classNumber1 = ee.Number(areaDict1.get('class')).format()
var area1 = ee.Number(
areaDict1.get('sum')).divide(1e6)//.round()
return ee.List([classNumber1, area1])
})
var result2 = ee.Dictionary(classAreaLists1.flatten())
print(result2)
后记
这个原理虽然比较简单,但是也到处找了很久代码,在设置阈值时也做了很多的前处理分析,虽然最后表明人工阈值不太可行,但是在大致的趋势上面还是足够明显的。不过果然偷懒还是不太行,老老实实地回去打标签做随机森林分类了。。