|
GEE实战 | 各种指数分析（植被、水体、建筑）
指数植被指数(NDVI | SAVI | IVI)代码对比结论
水体指数(NDWI | MNDWI)代码对比结论
建筑指数(NDBI | BSI | EIBI | NDISI)代码对比结论
总代码（指数计算+可视化(dual) map）
该文章介绍下在GEE平台下，使用Landsat 8 SR波段信息提取各种指数，并可视化对比讨论分类效果。
指数
全部以Landsat 8 SR波段数据进行计算，用CLC作为土地覆盖分类参考数据库。 P.S. 各个指数的代码主要是计算指数，可视化见总代码~
植被指数(NDVI | SAVI | IVI)
比较了NDVI、SAVI、IVI，并以CLC做对比选择合适的用来提取相关地物的波段。
代码
var l8sr = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
function maskL8sr(image) {
var cloudShadowBitMask = 1
red:image.select('B4'),
NIR:image.select('B5')
}).float().rename('SAVI');
var ivi = image.expression(
'(savi-(ndbi+mndwi)/2)/(savi+(ndbi+mndwi)/2)',{
ndbi:image.select('NDBI'),
savi:image.select('SAVI'),
mndwi:image.select('MNDWI')
}).float().rename('IVI');
image = image.addBands([ndvi,savi,ivi]);
};
对比
先放上CLC数据库 注： 水体 - 蓝色 植被 - 绿色 作物 - 黄色 建筑 - 红色 湿地 - 灰色
   
结论
SAVI和NDVI有时候对于植被和农作物是相反的，迷惑.在没找到明显区别的情况下，认为NDVI更合适
水体指数(NDWI | MNDWI)
水体指数采用以下两种NDWI和MNDWI：
N
D
W
I
=
(
G
r
e
e
n
−
N
I
R
)
/
(
G
r
e
e
n
+
N
I
R
)
NDWI=(Green-NIR)/(Green+NIR)
NDWI=(Green−NIR)/(Green+NIR)
M
N
D
W
I
=
(
G
r
e
e
n
−
M
I
R
)
/
(
G
r
e
e
n
+
M
I
R
)
MNDWI=(Green-MIR)/(Green+MIR)
MNDWI=(Green−MIR)/(Green+MIR)
其中，Green为绿光波段，NIR为近红外波段，MIR为中红外波段。其水体指数的阈值是[0，1]。
代码
var l8sr = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
function maskL8sr(image) {
var cloudShadowBitMask = 1
var cloudShadowBitMask = 1
green:image.select('B3'),
NIR:image.select('B5'),
SWIR2:image.select('B7')
}).float().rename('NDBBI');
var savi = image.expression(
'((NIR-red)*1.5)/(NIR+red+0.5)',{
red:image.select('B4'),
NIR:image.select('B5')
}).float().rename('SAVI');
var bsi = image.expression(
'((SWIR1+red)-(NIR+blue))/((SWIR1+red)+(NIR+blue))',{
blue:image.select('B2'),
red:image.select('B4'),
NIR:image.select('B5'),
SWIR1:image.select('B6')
}).float().rename('BSI');
image = image.addBands([ndbi,mndwi,ndbbi,savi,bsi]);
var ebsi = image.expression(
'(bsi-mndwi)/(bsi+mndwi)',{
bsi:image.select('BSI'),
mndwi:image.select('MNDWI')
}).float().rename('EBSI');
image = image.addBands(ebsi);
var eibi = image.expression(
'(ndbbi-(4*ebsi+savi+mndwi)/6)/(ndbbi+(4*ebsi+savi+mndwi)/6)',{
ndbbi:image.select('NDBBI'),
ebsi:image.select('EBSI'),
savi:image.select('SAVI'),
mndwi:image.select('MNDWI')
}).float().rename('EIBI');
var ndisi = image.expression(
'(TIR-(mndwi+NIR+SWIR)/3)/(TIR+(mndwi+NIR+SWIR)/3)',{
NIR:image.select('B5'),
SWIR:image.select('B6'),
mndwi:image.select('MNDWI'),
TIR:image.select('B10'),
}).float().rename('NDISI');
return image.addBands([eibi,ndisi]);//,
};
对比
先放上CLC数据库 注： 水体 - 蓝色 植被 - 绿色 作物 - 黄色 建筑 - 红色 湿地 - 灰色   
结论
NDBI适合提取农作物和城市EIBI中，植被与农作物+建筑+水体对比明显可视化对比不明显下，选择简单的NDBI更合适
总代码（指数计算+可视化(dual) map）
var dataset = ee.Image("COPERNICUS/CORINE/V20/100m/2018"),
var l8sr = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
function maskL8sr(image) {
var cloudShadowBitMask = 1
green:image.select('B3'),
NIR:image.select('B5'),
SWIR1:image.select('B6'),
SWIR2:image.select('B7'),
}).float().rename('AWEI');
var ndwi = image.normalizedDifference(['B3','B5']).rename('NDWI');
var ndvi = image.normalizedDifference(['B5','B4']).rename('NDVI');
var ndbi = image.normalizedDifference(['B6','B5']).rename('NDBI');
var mndwi = image.normalizedDifference(['B3','B6']).rename('MNDWI');
var lswi = image.normalizedDifference(['B5','B6']).rename('LSWI');
var ndbbi = image.expression(
'(1.5*SWIR2-(NIR+green)/2)/(1.5*SWIR2+(NIR+green)/2)',{
green:image.select('B3'),
NIR:image.select('B5'),
SWIR2:image.select('B7')
}).float().rename('NDBBI');
var savi = image.expression(
'((NIR-red)*1.5)/(NIR+red+0.5)',{
red:image.select('B4'),
NIR:image.select('B5')
}).float().rename('SAVI');
var bsi = image.expression(
'((SWIR1+red)-(NIR+blue))/((SWIR1+red)+(NIR+blue))',{
blue:image.select('B2'),
red:image.select('B4'),
NIR:image.select('B5'),
SWIR1:image.select('B6')
}).float().rename('BSI');
image = image.addBands([awei,ndvi,ndbi,ndwi,mndwi,lswi,ndbbi,savi,bsi]);
var ebsi = image.expression(
'(bsi-mndwi)/(bsi+mndwi)',{
bsi:image.select('BSI'),
mndwi:image.select('MNDWI')
}).float().rename('EBSI');
image = image.addBands(ebsi);
var eibi = image.expression(
'(ndbbi-(4*ebsi+savi+mndwi)/6)/(ndbbi+(4*ebsi+savi+mndwi)/6)',{
ndbbi:image.select('NDBBI'),
ebsi:image.select('EBSI'),
savi:image.select('SAVI'),
mndwi:image.select('MNDWI')
}).float().rename('EIBI');
var ivi = image.expression(
'(savi-(ndbi+mndwi)/2)/(savi+(ndbi+mndwi)/2)',{
ndbi:image.select('NDBI'),
savi:image.select('SAVI'),
mndwi:image.select('MNDWI')
}).float().rename('IVI');
var ndisi = image.expression(
'(TIR-(mndwi+NIR+SWIR)/3)/(TIR+(mndwi+NIR+SWIR)/3)',{
NIR:image.select('B5'),
SWIR:image.select('B6'),
mndwi:image.select('MNDWI'),
TIR:image.select('B10'),
}).float().rename('NDISI');
return image.addBands([eibi,ivi,ndisi]);//,
};
var clcVis_5 = {bands:'landcover',min: 0.0,max: 4.0, palette: ["ff0000","ffff00","00ff00","CCCCFF","0000ff"]}; // 建筑 农业 植被 湿地 水体
var Vis={bands:['B4', 'B3', 'B2'],min: 0.0,max: 0.3};
var visParams = {min: -0.8, max: 0.8, palette: [ 'green','white','blue']};
var visParams2 = {min: -0.8, max: 0.8, palette: [ '0000ff', '00ffff', 'ffff00', 'ff0000', 'ffffff']};
var visParams3 = {min: -0.8, max: 0.8, palette: [ 'FFF0F5','FFB6C1','FFC0CB','DC143C']};
var ndviVis = {min: -0.2,max: 0.8,palette: ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901', '66A000', '529400']};
var vis_flag=''//'split'//'map';
if (vis_flag=='map'){
Map.setCenter(7.1013, 61.6904, 10);
Map.addLayer(image,Vis,'Land8');
Map.addLayer(landimg.select('NDVI'),ndviVis,'NDVI',false);
Map.addLayer(landimg.select('SAVI'),ndviVis,'SAVI',false);
Map.addLayer(landimg.select('IVI'),ndviVis,'IVI',false);
// Map.setCenter(12.8461, 63.4398, 9);
// Map.setCenter(7.5135, 59.8554, 9);
// Map.setCenter(10.756, 59.9114, 9);
// Map.addLayer(landimg.select('AWEI'),visParams,'AWEI',false);
// Map.addLayer(landimg.select('NDWI'),visParams,'NDWI',false);
// Map.addLayer(landimg.select('MNDWI'),visParams,'MNDWI',false);
// Map.addLayer(landimg.select('LSWI'),visParams,'LSWI',false);
// Map.addLayer(landimg.select('NDBI'),visParams2,'NDBI',false);
// Map.addLayer(landimg.select('BSI'),visParams3,'BSI',false);
// // Map.addLayer(landimg.select('NDBBI'),visParams2,'NDBBI',false);
// Map.addLayer(landimg.select('EIBI'),visParams3,'EIBI',false);
// // Map.addLayer(landimg.select('EBSI'),visParams3,'EBSI',false);
// Split panel可视化
} else if (vis_flag=='split'){
var images ={
// 'Landsat':image.visualize(Vis),
'CLC_5class':landcover_5.visualize(clcVis_5),
// 'NDWI':landimg.select('NDWI').visualize(visParams),
// 'MNDWI':landimg.select('MNDWI').visualize(visParams),
// 'LSWI':landimg.select('LSWI').visualize(visParams),
// 'NDBI':landimg.select('NDBI').visualize(visParams),
// 'BSI':landimg.select('BSI').visualize(visParams),
// 'EIBI':landimg.select('EIBI').visualize(visParams),
// 'NDISI':landimg.select('NDISI').visualize(visParams),
'NDVI':landimg.select('NDVI').visualize(ndviVis),
'SAVI':landimg.select('SAVI').visualize(ndviVis),
'IVI':landimg.select('IVI').visualize(ndviVis),
};
// Create the left map, and have it display layer 0.
var leftMap = ui.Map();
leftMap.setControlVisibility(false);
var leftSelector = addLayerSelector(leftMap, 0, 'top-left');
// Create the right map, and have it display layer 1.
var rightMap = ui.Map();
rightMap.setControlVisibility(false);
var rightSelector = addLayerSelector(rightMap, 1, 'top-right');
function addLayerSelector(mapToChange, defaultValue, position) {
var label = ui.Label('Choose an image to visualize');
function updateMap(selection) {
mapToChange.layers().set(0, ui.Map.Layer(images[selection]));
}
var select = ui.Select({items: Object.keys(images), onChange: updateMap});
select.setValue(Object.keys(images)[defaultValue], true);
var controlPanel =
ui.Panel({widgets: [label, select], style: {position: position}});
mapToChange.add(controlPanel);
}
var splitPanel = ui.SplitPanel({
firstPanel: leftMap,
secondPanel: rightMap,
wipe: true,
style: {stretch: 'both'}
});
ui.root.widgets().reset([splitPanel]);
var linker = ui.Map.Linker([leftMap, rightMap]);
leftMap.setCenter(12.5181, 55.7815, 8);
} else {
print()
}
|