概述
文章目录
- 前言
- 第一步:定义研究区,自行更换自己的研究区
- 第二步:加载Landsat 8,LST影像数据集合
- 第三步:构造SPWI指数函数: (B5 - B7 + B2)/(B5 + B7 + B2);构造NDLI指数函数: (B3 - B4)/(B3 + B4 + B6);构造地表温度LST函数;构造RVI指数函数:B5/B4;构造NDSI指数函数: (B6 - B5)/(B6 + B5)
- 第四步:构造指标熵计算函数
- 第五步:主函数开始,引用指标计算函数,导入数据,计算各生态指标
- 第六步:引用指标熵计算公式,导入指标,计算指标信息熵
- 第七步:利用指标信息熵,求取各指标权重
- 第八步:加权融合,获取最终的生态环境质量结果
- 第九步:结果显示
- 完整代码
- 总结
前言
基于GEE平台,实现顾及水效益的生态环境质量评价方法。
《Water Benefit-based Ecological Index for Urban Ecological Environment Quality Assessments》
论文链接:https://ieeexplore.ieee.org/document/9492814
运行结果:
第一步:定义研究区,自行更换自己的研究区
代码如下(示例):
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'
};
//定义降采样参数,便于计算
var reduceScale = 500;
// 定义合成影像开始结束时间
var date_start = '2019-07-01';
var date_end = '2019-10-30';
第二步:加载Landsat 8,LST影像数据集合
代码如下(示例):
//第二步:加载Landsat 8影像数据集合
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影像数据集合
var LSTimg = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
.filterBounds(geometry0)
.filterDate(date_start, date_end)
.map(maskL8sr)
.median()
.clip(geometry0);
第三步:构造SPWI指数函数: (B5 - B7 + B2)/(B5 + B7 + B2);构造NDLI指数函数: (B3 - B4)/(B3 + B4 + B6);构造地表温度LST函数;构造RVI指数函数:B5/B4;构造NDSI指数函数: (B6 - B5)/(B6 + B5)
代码如下(示例):
//第三步:构造SPWI指数函数: (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;
}
//第三步:构造NDLI指数函数: (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;
}
//第三步:构造地表温度LST函数
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;
}
//第三步:构造RVI指数函数: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;
}
//第三步:构造NDSI指数函数: (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;
}
第四步:构造指标熵计算函数
代码如下(示例):
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;
}
第五步:主函数开始,引用指标计算函数,导入数据,计算各生态指标
代码如下(示例):
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");
第六步:引用指标熵计算公式,导入指标,计算指标信息熵
代码如下(示例):
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);
第七步:利用指标信息熵,求取各指标权重
代码如下(示例):
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,指标权重为:',w01,w02,w03,w04,w05);
第八步:加权融合,获取最终的生态环境质量结果
代码如下(示例):
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));
第九步:结果显示
代码如下(示例):
//图例,红-蓝,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");
完整代码
代码如下(示例):
//《Water Benefit-based Ecological Index for Urban Ecological Environment Quality Assessments》
// 顾及水效益的生态环境质量评价方法,论文链接:https://ieeexplore.ieee.org/document/9492814
// 第一步:定义研究区,自行更换自己的研究区
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'
};
//定义降采样参数,便于计算
var reduceScale = 500;
// 定义合成影像开始结束时间
var date_start = '2019-07-01';
var date_end = '2019-10-30';
//第二步:加载Landsat 8影像数据集合
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影像数据集合
var LSTimg = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")
.filterBounds(geometry0)
.filterDate(date_start, date_end)
.map(maskL8sr)
.median()
.clip(geometry0);
//第三步:构造SPWI指数函数: (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;
}
//第三步:构造NDLI指数函数: (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;
}
//第三步:构造地表温度LST函数
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;
}
//第三步:构造RVI指数函数: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;
}
//第三步:构造NDSI指数函数: (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;
}
//第四步:构造指标熵计算函数
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;
}
//第五步:主函数开始,引用指标计算函数,导入数据,计算各生态指标
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");
// 第六步,引用指标熵计算公式,导入指标,计算指标信息熵
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);
// 第七步,利用指标信息熵,求取各指标权重
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,指标权重为:',w01,w02,w03,w04,w05);
// 第八步,加权融合,获取最终的生态环境质量结果
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));
//第九步,结果显示
//图例,红-蓝,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");
总结
提示:用完后记得给个反馈:
论文题目:《Water Benefit-based Ecological Index for Urban Ecological Environment Quality Assessments》
论文链接:https://ieeexplore.ieee.org/document/9492814
最后
以上就是满意日记本为你收集整理的利用GEE计算遥感生态指数(WBEI)——Landsat 8为例前言第一步:定义研究区,自行更换自己的研究区第二步:加载Landsat 8,LST影像数据集合第三步:构造SPWI指数函数: (B5 - B7 + B2)/(B5 + B7 + B2);构造NDLI指数函数: (B3 - B4)/(B3 + B4 + B6);构造地表温度LST函数;构造RVI指数函数:B5/B4;构造NDSI指数函数: (B6 - B5)/(B6 + B5)第四步:构造指标熵计算函数第五步:主函数开始,引用指标计算的全部内容,希望文章能够帮你解决利用GEE计算遥感生态指数(WBEI)——Landsat 8为例前言第一步:定义研究区,自行更换自己的研究区第二步:加载Landsat 8,LST影像数据集合第三步:构造SPWI指数函数: (B5 - B7 + B2)/(B5 + B7 + B2);构造NDLI指数函数: (B3 - B4)/(B3 + B4 + B6);构造地表温度LST函数;构造RVI指数函数:B5/B4;构造NDSI指数函数: (B6 - B5)/(B6 + B5)第四步:构造指标熵计算函数第五步:主函数开始,引用指标计算所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复