前言
基于GEE平臺，實現(xiàn)顧及水效益的生態(tài)環(huán)境質(zhì)量評價方法。
《Water Benefit-based Ecological Index for Urban Ecological Environment Quality Assessments》
論文鏈接：https://ieeexplore.ieee.org/document/9492814
運行結(jié)果：
第一步：定義研究區(qū)，自行更換自己的研究區(qū)
代碼如下（示例）：
var geometry0 =
/* color: #98ff00 */
/* displayProperties: [
{
"type": "rectangle"
}
] */
ee.Geometry.Polygon(
[[[120.1210075098537, 35.975189051414006],
[120.1210075098537, 35.886229778229115],
[120.25764996590839, 35.886229778229115],
[120.25764996590839, 35.975189051414006]]], null, false);
Map.centerObject(geometry0);
//圖例，紅-藍，1-0
var visParam = {
palette: '040274, 040281, 0502a3, 0502b8, 0502ce, 0502e6, 0602ff, 235cb1, 307ef3, 269db1, 30c8e2, 32d3ef, 3be285, 3ff38f, 86e26f, 3ae237, b5e22e, d6e21f, fff705, ffd611, ffb613, ff8b13, ff6e08, ff500d, ff0000, de0101, c21301, a71001, 911003'
};
//定義降采樣參數(shù)，便于計算
var reduceScale = 500;
// 定義合成影像開始結(jié)束時間
var date_start = '2019-07-01';
var date_end = '2019-10-30';
第二步：加載Landsat 8，LST影像數(shù)據(jù)集合
代碼如下（示例）：
//第二步:加載Landsat 8影像數(shù)據(jù)集合
var img1 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(geometry0)
.filterDate(date_start, date_end)
.filterMetadata('CLOUD_COVER', 'less_than',30)
.select(["B2", "B3", "B4", "B5", "B6", "B7", "B10"])
//.map(removeCloud)
.mean()
.clip(geometry0);
//第二步:加載LST影像數(shù)據(jù)集合
var LSTimg = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
.filterBounds(geometry0)
.filterDate(date_start, date_end)
.map(maskL8sr)
.median()
.clip(geometry0);
第三步：構(gòu)造SPWI指數(shù)函數(shù): (B5 - B7 + B2)/(B5 + B7 + B2)；構(gòu)造NDLI指數(shù)函數(shù): (B3 - B4)/(B3 + B4 + B6)；構(gòu)造地表溫度LST函數(shù)；構(gòu)造RVI指數(shù)函數(shù):B5/B4；構(gòu)造NDSI指數(shù)函數(shù): (B6 - B5)/(B6 + B5)
代碼如下（示例）：
//第三步：構(gòu)造SPWI指數(shù)函數(shù): (B5 - B7 + B2)/(B5 + B7 + B2)
function SPWI(image,roi,reduceScale) {
var nir = image.select('B5');
var swir2 = image.select('B7');
var blue = image.select('B2');
var SPWI0 = (nir.subtract(swir2).add(blue)).divide(nir.add(swir2).add(blue)).rename("SPWI")
var dic = SPWI0.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("SPWI_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("SPWI_p95"));
var SPWI = SPWI0.where(SPWI0.select('SPWI').gt(p95), p95).where(SPWI0.select('SPWI').lt(p5), p5);
var nSPWI = SPWI.subtract(p5).divide(ee.Number(p95).subtract(p5));
return nSPWI;
}
//第三步：構(gòu)造NDLI指數(shù)函數(shù): (B3 - B4)/(B3 + B4 + B6)
function NDLI(image,roi,reduceScale) {
var green = image.select('B3');
var red = image.select('B4');
var swir1 = image.select('B6');
var NDLI0 = (green.subtract(red)).divide(green.add(red).add(swir1)).rename("NDLI")
var dic = NDLI0.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("NDLI_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("NDLI_p95"));
var NDLI = NDLI0.where(NDLI0.select('NDLI').gt(p95), p95).where(NDLI0.select('NDLI').lt(p5), p5);
var nNDLI = NDLI.subtract(p5).divide(ee.Number(p95).subtract(p5));
print(nNDLI)
return nNDLI;
}
//第三步：構(gòu)造地表溫度LST函數(shù)
function maskL8sr(image) {
// var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0); //這里是去云
// var saturationMask = image.select('QA_RADSAT').eq(0); //這里是去云
var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2); // 這里就是對以SR開頭的波段進行一個計算，讓它回到一個正常的值，以便我們使用。
var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0); // 這里就是對我們要得到的地表溫度波段進行計算。
return image.addBands(opticalBands, null, true)
.addBands(thermalBands, null, true)
// .updateMask(qaMask)
// .updateMask(saturationMask);
}
function LST(image,roi,reduceScale) {
var img = image.select("ST_B10")
var lst = img.expression(
'B1-273.15',
{
B1:img.select('ST_B10'),
}).rename("LST");
var dic = lst.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("LST_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("LST_p95"));
print('最低溫度為:',p5);
print('最高溫度為:',p95);
var LST0 = lst.where(lst.select('LST').gt(p95), p95).where(lst.select('LST').lt(p5), p5);
var nLST = LST0.subtract(p5).divide(ee.Number(p95).subtract(p5));
return nLST;
}
//第三步：構(gòu)造RVI指數(shù)函數(shù):B5/B4
function RVI(image,roi,reduceScale) {
var NIR = image.select(['B5']);
var R = image.select(['B4']);
var rvi0 = NIR.divide(R).rename('RVI');
var dic = rvi0.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("RVI_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("RVI_p95"));
var RVI = rvi0.where(rvi0.select('RVI').gt(p95), p95).where(rvi0.select('RVI').lt(p5), p5);
var nRVI = RVI.subtract(p5).divide(ee.Number(p95).subtract(p5));
return nRVI;
}
//第三步：構(gòu)造NDSI指數(shù)函數(shù): (B6 - B5)/(B6 + B5)
function NDSI(image,roi,reduceScale) {
var NDSI0 = image.select("B6").subtract(image.select("B5"))
.divide(image.select("B6").add(image.select("B5"))).rename("NDSI")
var dic = NDSI0.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("NDSI_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("NDSI_p95"));
var NDSI = NDSI0.where(NDSI0.select('NDSI').gt(p95), p95).where(NDSI0.select('NDSI').lt(p5), p5);
var nNDSI = NDSI.subtract(p5).divide(ee.Number(p95).subtract(p5));
return nNDSI;
}
第四步：構(gòu)造指標熵計算函數(shù)
代碼如下（示例）：
function Ex(xi,roi,reduceScale) {
var xs = xi.reduceRegion(ee.Reducer.sum(), roi, reduceScale).get("XI").getInfo();
var f = xi.divide(xs);
var fs = (f.multiply(f.log())).reduceRegion(ee.Reducer.sum(), roi, reduceScale).get("XI").getInfo();
var n = xi.reduceRegion(ee.Reducer.count(), roi, reduceScale).get("XI").getInfo();
var lnn = Math.log(n);
var e = ee.Number(-fs).divide(lnn);
return e;
}
第五步：主函數(shù)開始，引用指標計算函數(shù)，導(dǎo)入數(shù)據(jù)，計算各生態(tài)指標
代碼如下（示例）：
var spwi0 = SPWI(img1, geometry0,reduceScale).rename("XI");
var ndli0 = NDLI(img1, geometry0,reduceScale).rename("XI");
var lst0 = LST(LSTimg, geometry0,reduceScale).rename("XI");
var rvi0 = RVI(img1, geometry0,reduceScale).rename("XI");
var ndsi0 = NDSI(img1, geometry0,reduceScale).rename("XI");
第六步：引用指標熵計算公式，導(dǎo)入指標，計算指標信息熵
代碼如下（示例）：
var e01 = Ex(spwi0,geometry0,reduceScale);
var e02 = Ex(ndli0,geometry0,reduceScale);
var e03 = Ex(lst0,geometry0,reduceScale);
var e04 = Ex(rvi0,geometry0,reduceScale);
var e05 = Ex(ndsi0,geometry0,reduceScale);
第七步：利用指標信息熵，求取各指標權(quán)重
代碼如下（示例）：
var w01 = ee.Number(1).subtract(e01).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
var w02 = ee.Number(1).subtract(e02).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
var w03 = ee.Number(1).subtract(e03).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
var w04 = ee.Number(1).subtract(e04).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
var w05 = ee.Number(1).subtract(e05).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
print('SPWI,NDLI,LST,RVI,NDSI,指標權(quán)重為：',w01,w02,w03,w04,w05);
第八步：加權(quán)融合，獲取最終的生態(tài)環(huán)境質(zhì)量結(jié)果
代碼如下（示例）：
var wbei0 = spwi0.multiply(w01)
.add(ndli0.multiply(w02))
.subtract(lst0.multiply(w03))
.add(rvi0.multiply(w04))
.subtract(ndsi0.multiply(w05))
.rename("WBEI0");
var dic = wbei0.reduceRegion(ee.Reducer.percentile([5, 95]), geometry0, reduceScale); //計算過程中將30m降采樣到100m，減少計算
var p5 = ee.Number(ee.Dictionary(dic).get("WBEI0_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("WBEI0_p95"));
var WBEI = wbei0.where(wbei0.select('WBEI0').gt(p95), p95).where(wbei0.select('WBEI0').lt(p5), p5);
var nWBEI = WBEI.subtract(p5).divide(ee.Number(p95).subtract(p5));
第九步：結(jié)果顯示
代碼如下（示例）：
//圖例，紅-藍，0-1
var visParam1 = {
min: 0.0, max: 1.0 ,palette: '911003, a71001, c21301, de0101, ff0000, ff500d, ff6e08, ff8b13,ffb613, ffd611, fff705, d6e21f, b5e22e, 3ae237, 86e26f, 3ff38f,3be285, 32d3ef, 30c8e2, 269db1, 307ef3,235cb1, 0602ff, 0502e6, 0502ce,0502b8,0502a3, 040281, 040274 '
};
Map.addLayer(nWBEI, visParam1, "nWBEI");
完整代碼
代碼如下（示例）：文章來源：http://www.zghlxwxcb.cn/news/detail-518946.html
//《Water Benefit-based Ecological Index for Urban Ecological Environment Quality Assessments》
// 顧及水效益的生態(tài)環(huán)境質(zhì)量評價方法，論文鏈接：https://ieeexplore.ieee.org/document/9492814
// 第一步:定義研究區(qū)，自行更換自己的研究區(qū)
var geometry0 =
/* color: #98ff00 */
/* displayProperties: [
{
"type": "rectangle"
}
] */
ee.Geometry.Polygon(
[[[120.1210075098537, 35.975189051414006],
[120.1210075098537, 35.886229778229115],
[120.25764996590839, 35.886229778229115],
[120.25764996590839, 35.975189051414006]]], null, false);
Map.centerObject(geometry0);
//圖例，紅-藍，1-0
var visParam = {
palette: '040274, 040281, 0502a3, 0502b8, 0502ce, 0502e6, 0602ff, 235cb1, 307ef3, 269db1, 30c8e2, 32d3ef, 3be285, 3ff38f, 86e26f, 3ae237, b5e22e, d6e21f, fff705, ffd611, ffb613, ff8b13, ff6e08, ff500d, ff0000, de0101, c21301, a71001, 911003'
};
//定義降采樣參數(shù)，便于計算
var reduceScale = 500;
// 定義合成影像開始結(jié)束時間
var date_start = '2019-07-01';
var date_end = '2019-10-30';
//第二步:加載Landsat 8影像數(shù)據(jù)集合
var img1 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterBounds(geometry0)
.filterDate(date_start, date_end)
.filterMetadata('CLOUD_COVER', 'less_than',30)
.select(["B2", "B3", "B4", "B5", "B6", "B7", "B10"])
//.map(removeCloud)
.mean()
.clip(geometry0);
//第二步:加載LST影像數(shù)據(jù)集合
var LSTimg = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
.filterBounds(geometry0)
.filterDate(date_start, date_end)
.map(maskL8sr)
.median()
.clip(geometry0);
//第三步：構(gòu)造SPWI指數(shù)函數(shù): (B5 - B7 + B2)/(B5 + B7 + B2)
function SPWI(image,roi,reduceScale) {
var nir = image.select('B5');
var swir2 = image.select('B7');
var blue = image.select('B2');
var SPWI0 = (nir.subtract(swir2).add(blue)).divide(nir.add(swir2).add(blue)).rename("SPWI")
var dic = SPWI0.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("SPWI_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("SPWI_p95"));
var SPWI = SPWI0.where(SPWI0.select('SPWI').gt(p95), p95).where(SPWI0.select('SPWI').lt(p5), p5);
var nSPWI = SPWI.subtract(p5).divide(ee.Number(p95).subtract(p5));
return nSPWI;
}
//第三步：構(gòu)造NDLI指數(shù)函數(shù): (B3 - B4)/(B3 + B4 + B6)
function NDLI(image,roi,reduceScale) {
var green = image.select('B3');
var red = image.select('B4');
var swir1 = image.select('B6');
var NDLI0 = (green.subtract(red)).divide(green.add(red).add(swir1)).rename("NDLI")
var dic = NDLI0.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("NDLI_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("NDLI_p95"));
var NDLI = NDLI0.where(NDLI0.select('NDLI').gt(p95), p95).where(NDLI0.select('NDLI').lt(p5), p5);
var nNDLI = NDLI.subtract(p5).divide(ee.Number(p95).subtract(p5));
print(nNDLI)
return nNDLI;
}
//第三步：構(gòu)造地表溫度LST函數(shù)
function maskL8sr(image) {
// var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0); //這里是去云
// var saturationMask = image.select('QA_RADSAT').eq(0); //這里是去云
var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2); // 這里就是對以SR開頭的波段進行一個計算，讓它回到一個正常的值，以便我們使用。
var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0); // 這里就是對我們要得到的地表溫度波段進行計算。
return image.addBands(opticalBands, null, true)
.addBands(thermalBands, null, true)
// .updateMask(qaMask)
// .updateMask(saturationMask);
}
function LST(image,roi,reduceScale) {
var img = image.select("ST_B10")
var lst = img.expression(
'B1-273.15',
{
B1:img.select('ST_B10'),
}).rename("LST");
var dic = lst.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("LST_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("LST_p95"));
print('最低溫度為:',p5);
print('最高溫度為:',p95);
var LST0 = lst.where(lst.select('LST').gt(p95), p95).where(lst.select('LST').lt(p5), p5);
var nLST = LST0.subtract(p5).divide(ee.Number(p95).subtract(p5));
return nLST;
}
//第三步：構(gòu)造RVI指數(shù)函數(shù):B5/B4
function RVI(image,roi,reduceScale) {
var NIR = image.select(['B5']);
var R = image.select(['B4']);
var rvi0 = NIR.divide(R).rename('RVI');
var dic = rvi0.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("RVI_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("RVI_p95"));
var RVI = rvi0.where(rvi0.select('RVI').gt(p95), p95).where(rvi0.select('RVI').lt(p5), p5);
var nRVI = RVI.subtract(p5).divide(ee.Number(p95).subtract(p5));
return nRVI;
}
//第三步：構(gòu)造NDSI指數(shù)函數(shù): (B6 - B5)/(B6 + B5)
function NDSI(image,roi,reduceScale) {
var NDSI0 = image.select("B6").subtract(image.select("B5"))
.divide(image.select("B6").add(image.select("B5"))).rename("NDSI")
var dic = NDSI0.reduceRegion(ee.Reducer.percentile([5, 95]), roi, reduceScale);
var p5 = ee.Number(ee.Dictionary(dic).get("NDSI_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("NDSI_p95"));
var NDSI = NDSI0.where(NDSI0.select('NDSI').gt(p95), p95).where(NDSI0.select('NDSI').lt(p5), p5);
var nNDSI = NDSI.subtract(p5).divide(ee.Number(p95).subtract(p5));
return nNDSI;
}
//第四步：構(gòu)造指標熵計算函數(shù)
function Ex(xi,roi,reduceScale) {
var xs = xi.reduceRegion(ee.Reducer.sum(), roi, reduceScale).get("XI").getInfo();
var f = xi.divide(xs);
var fs = (f.multiply(f.log())).reduceRegion(ee.Reducer.sum(), roi, reduceScale).get("XI").getInfo();
var n = xi.reduceRegion(ee.Reducer.count(), roi, reduceScale).get("XI").getInfo();
var lnn = Math.log(n);
var e = ee.Number(-fs).divide(lnn);
return e;
}
//第五步：主函數(shù)開始，引用指標計算函數(shù)，導(dǎo)入數(shù)據(jù)，計算各生態(tài)指標
var spwi0 = SPWI(img1, geometry0,reduceScale).rename("XI");
var ndli0 = NDLI(img1, geometry0,reduceScale).rename("XI");
var lst0 = LST(LSTimg, geometry0,reduceScale).rename("XI");
var rvi0 = RVI(img1, geometry0,reduceScale).rename("XI");
var ndsi0 = NDSI(img1, geometry0,reduceScale).rename("XI");
// 第六步，引用指標熵計算公式，導(dǎo)入指標，計算指標信息熵
var e01 = Ex(spwi0,geometry0,reduceScale);
var e02 = Ex(ndli0,geometry0,reduceScale);
var e03 = Ex(lst0,geometry0,reduceScale);
var e04 = Ex(rvi0,geometry0,reduceScale);
var e05 = Ex(ndsi0,geometry0,reduceScale);
// 第七步，利用指標信息熵，求取各指標權(quán)重
var w01 = ee.Number(1).subtract(e01).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
var w02 = ee.Number(1).subtract(e02).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
var w03 = ee.Number(1).subtract(e03).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
var w04 = ee.Number(1).subtract(e04).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
var w05 = ee.Number(1).subtract(e05).divide(ee.Number(5).subtract(e01).subtract(e02).subtract(e03).subtract(e04).subtract(e05));
print('SPWI,NDLI,LST,RVI,NDSI,指標權(quán)重為：',w01,w02,w03,w04,w05);
// 第八步，加權(quán)融合，獲取最終的生態(tài)環(huán)境質(zhì)量結(jié)果
var wbei0 = spwi0.multiply(w01)
.add(ndli0.multiply(w02))
.subtract(lst0.multiply(w03))
.add(rvi0.multiply(w04))
.subtract(ndsi0.multiply(w05))
.rename("WBEI0");
var dic = wbei0.reduceRegion(ee.Reducer.percentile([5, 95]), geometry0, reduceScale); //計算過程中將30m降采樣到100m，減少計算
var p5 = ee.Number(ee.Dictionary(dic).get("WBEI0_p5"));
var p95 = ee.Number(ee.Dictionary(dic).get("WBEI0_p95"));
var WBEI = wbei0.where(wbei0.select('WBEI0').gt(p95), p95).where(wbei0.select('WBEI0').lt(p5), p5);
var nWBEI = WBEI.subtract(p5).divide(ee.Number(p95).subtract(p5));
//第九步,結(jié)果顯示
//圖例，紅-藍，0-1
var visParam1 = {
min: 0.0, max: 1.0 ,palette: '911003, a71001, c21301, de0101, ff0000, ff500d, ff6e08, ff8b13,ffb613, ffd611, fff705, d6e21f, b5e22e, 3ae237, 86e26f, 3ff38f,3be285, 32d3ef, 30c8e2, 269db1, 307ef3,235cb1, 0602ff, 0502e6, 0502ce,0502b8,0502a3, 040281, 040274 '
};
Map.addLayer(nWBEI, visParam1, "nWBEI");
總結(jié)
提示：用完后記得給個反饋：
論文題目：《Water Benefit-based Ecological Index for Urban Ecological Environment Quality Assessments》
論文鏈接：https://ieeexplore.ieee.org/document/9492814文章來源地址http://www.zghlxwxcb.cn/news/detail-518946.html
到了這里，關(guān)于利用GEE計算遙感生態(tài)指數(shù)(WBEI)——Landsat 8為例的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
