太湖水体浑浊度计算

基于PIE-Engine获取2013-2020年覆盖太湖水域的Landsat 8 Surface Reflectance遥感影像数据，Landsat 8 SR数据集是OLI/TIRS传感器数据经过大气校正的表面反射率数据，分辨率为30米，重放周期为16天。影像包含5个可见光波段、1个近红外（NIR）波段和2个短波红外（SWIR）波段，经过正射校正后的表面反射率，以及两个经过正射校正后的亮度温度热红外（TIR）波段。

相关研究指出，Landsat 8数据的近红外波段和红光波段分别与蓝藻密度和浑浊度具有较高的相关系数，但相关程度都不高，主要原因有：① 水华与水草区域对水质反演的影响；② 水表面光滑度和微波随时间和空间的干扰；③ 陆地光谱对近岸边水体的干扰；④ 水体中较高浓度悬浮物的高反射干扰。在去除影响区域和异常数据，并利用普通水体与不同波段的组合构建相关性模型后发现，蓝藻密度与B5/B4以及浑浊度与(B4-B3)/(B4+B3)的波段组合具有最高的相关系数。相关性模型如下：

1) 在地图中加载并显示太湖边界矢量

//输入太湖边界矢量
var taihu = pie.FeatureCollection('user/101/public/shape/TaiHuShuiYu')
               .first()
               .geometry();
print(taihu);
Map.centerObject(taihu, 8);
Map.addLayer(taihu,{color:'FFFF00FF',width:2,fillColor:'00000000'},"taihu");复制错误复制成功

运行加载结果如下：

2) 筛选影像数据

根据时间、区域以及云量筛选Landsat8 SR数据集。

//加载Landsat 8 SR影像数据 
var l8Col = l8.filterDate("2018-01-01", "2020-11-01")
               .filterBounds(taihu)
               .filter(pie.Filter.lte('cloudCover', 5))
               .map(function(image) {
                   return image.set("id", image.id());
               });复制错误复制成功

3) 计算太湖水体浑浊度

根据公式循环计算太湖水体浑浊度。

// 设置图层显示属性
var visParam_t = {
    min: 0, 
    max: 180,
    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 dataInfo = l8Col.reduceColumns(pie.Reducer.toList(2), ["date", "id"]);

dataInfo.getInfo(function(dataObj) {
    print("影像列表信息", dataObj);
    var datas = dataObj.list;
    var dates = [];
    var t_images = [];
    //循环计算
    for (var i=0; i<datas.length; i++) {
        var date = datas[i].date;
        var id = datas[i].id;
        var image = pie.Image(id);
        var b3 = image.select("B3");
        var b4 = image.select("B4");
        // 根据浊度公式进行计算
        var turbidity = (((b4.subtract(b3)).divide(b4.add(b3)).power(2)).multiply(3117.4))
                       .add(((b4.subtract(b3)).divide(b4.add(b3))).multiply(1083.6)).add(106.17)
                       .clip(taihu);
        Map.addLayer(turbidity, visParam_t, date, false);
        var t_result = turbidity.reduceRegion(pie.Reducer.mean(), taihu, 30);
        dates.push(date);
        t_images.push(t_result);
    }  复制错误复制成功

4) 设置图表属性

//设置图表属性  
var line_t = {
    title: '太湖水域浑浊度动态变化',  
              legend: ['浑浊度'],  
              xAxisName: "日期",  
              yAxisName: "浑浊度",  
              chartType: "line",   
              smooth: true
};  
//显示折线图  
ChartImage(tbd_images, dates, line_t);  
//动画显示  
Map.playLayersAnimation(dates, 0.5, -1); 复制错误复制成功

5) 添加图例

添加太湖水体浑浊度渐变图例。

var data = {
    title: "浑浊度",
    colors: [
        '#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'
    ],
    labels: ["清澈","轻度浑浊","中度浑浊","重度浑浊"],
    step: 30
};
var style = {
    right: "150px",
    bottom: "10px",
    height: "70px",
    width: "350px"
};
var legend = ui.Legend(data, style);  
Map.addUI(legend);  复制错误复制成功

运行结果如下所示：

动态变化展示如下：

6) 完整代码

/**
 * @File    :   PIETaiHuTurbidity
 * @Time    :   2020/7/21
 * @Author  :   piesat
 * @Version :   1.0
 * @Contact :   400-890-0662
 * @License :   (C)Copyright 航天宏图信息技术股份有限公司
 * @Desc    :   太湖水体浑浊度计算
 */
// 导入landsat8
var l8 = pie.ImageCollection("LC08/01/T1_SR");
Map.setCenter(120.14928092430432,31.232446551838038, 8);

// 过滤范围&日期&云
var taihu = pie.FeatureCollection('user/101/public/shape/TaiHuShuiYu')
               .first()
               .geometry();
var l8Col = l8.filterDate("2018-01-01", "2020-11-01")
               .filterBounds(taihu)
               .filter(pie.Filter.lte('cloudCover', 5))
               .map(function(image) {
                   return image.set("id", image.id());
               });;
               
// 设置图层显示属性
var visParam_t = {
    min: 0, 
    max: 180,
    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 dataInfo = l8Col.reduceColumns(pie.Reducer.toList(2), ["date", "id"]);
dataInfo.getInfo(function(dataObj) {
    print("影像列表信息", dataObj);
    var datas = dataObj.list;
    var dates = [];
    var t_images = [];
    //计算
    for (var i=0; i<datas.length; i++) {
        var date = datas[i].date;
        var id = datas[i].id;
        var image = pie.Image(id);
        var b3 = image.select("B3");
        var b4 = image.select("B4");
        // 计算浊度
        var turbidity = (((b4.subtract(b3)).divide(b4.add(b3)).power(2)).multiply(3117.4))
                       .add(((b4.subtract(b3)).divide(b4.add(b3))).multiply(1083.6)).add(106.17)
                       .clip(taihu);
        Map.addLayer(turbidity, visParam_t, date, false);
        var t_result = turbidity.reduceRegion(pie.Reducer.mean(), taihu, 30);
        dates.push(date);
        t_images.push(t_result);
    }

    // 设置图表属性
    var line_t = {
    	title: '太湖水体浊度均值',
        legend: ['浊度'],
        xAxisName: "日期",
        yAxisName: "浊度",
        chartType: "line",
        // yMin:0,
        // yMax:100,
        smooth: true
    };

    // 显示折线图
    ChartImage(t_images, dates, line_t);

    // 动画显示
    Map.playLayersAnimation(dates, 0.5, -1);

    // 图例
    var data = {
        title: "浑浊度",
        colors: [
            '#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'
        ],
        labels: ["清澈","轻度浑浊","中度浑浊","重度浑浊"],
        step: 30
    };
    var style = {
        right: "150px",
        bottom: "10px",
        height: "70px",
        width: "350px"
    };
    var legend = ui.Legend(data, style);
    Map.addUI(legend); 
});