利用GEE计算遥感生态指数(WBEI)——Landsat 8为例

何处惹尘埃？

已于 2022-07-17 16:21:05 修改

阅读量5.9k

点赞数 16

分类专栏：生态环境文章标签： java 大数据云计算

于 2022-07-12 20:16:29 首次发布

本文链接：https://blog.csdn.net/qq_35490698/article/details/125750578

版权

生态环境专栏收录该内容

2 篇文章 12 订阅

订阅专栏

文章目录

前言
第一步：定义研究区，自行更换自己的研究区
第二步：加载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