太湖水体浑浊度计算

基于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) 核心代码

def Turbidity(B3, B4):
    temp1 = B4.subtract(B3).divide(B4.add(B3)).power(2).multiply(3117.4)
    temp2 = B4.subtract(B3).divide(B4.add(B3)).multiply(1083.6)
    turbidity = temp1.add(temp2).add(106.17)
    return turbidity
# 选择影像计算浑浊度
image = Imgcol.getAt(1).select(["B4", "B3"])
date = image.get("date").getInfo()
b3 = image.select("B3")
b4 = image.select("B4")
turbidity = Turbidity(b3, b4)
Map.addLayer(turbidity.clip(taihu), visParam_t, date)
Map复制错误复制成功

流程步骤如下：

a) 加载影像数据；

b) 定义模型公式；

c) 选择影像计算浑浊度，并加载到地图上。

2) 案例中主要涉及的算子

• getAt（）：根据索引选择影像；

• clip（）：根据矢量裁剪影像。

3) 结果分析展示

调用PIE-Engine遥感计算云平台在jupyter notebook中的实现效果如下图所示：

4) 完整代码

注意： 需要将代码中用户名和密码替换为用户在平台上注册时的用户名/密码

"""
Program infomation：Python3.7
"""
# 太湖浑浊度分析
import pie
from pie import *
from pie.chart import PIEChart
import ipywidgets
from ipyleaflet import WidgetControl
import random

# 使用用户在平台上注册时的用户名/密码登录
pie.initialize()

# 调用地图
Map = pie.Map()

# 加载太湖边界矢量
taihu = pie.FeatureCollection('user/101/public/shape/TaiHuShuiYu').first().geometry()
Map.setCenter(120.14928092430432,31.232446551838038,  8)
Map.addLayer(taihu,{"color": "red"}, "taihu")

# 加载Landsat 8 SR
Imgcol = pie.ImageCollection("LC08/01/T1_SR") \
            .filterDate("2017-01-01", "2020-11-01") \
            .filterBounds(taihu) \
            .filter(pie.Filter.lte('cloudCover', 5))

#设置图层显示属性
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'
    ]
}

def Turbidity(B3, B4):
    temp1 = B4.subtract(B3).divide(B4.add(B3)).power(2).multiply(3117.4)
    temp2 = B4.subtract(B3).divide(B4.add(B3)).multiply(1083.6)
    turbidity = temp1.add(temp2).add(106.17)
    return turbidity

# 选择影像计算浑浊度
image = Imgcol.getAt(1).select(["B4", "B3"])
date = image.get("date").getInfo()
b3 = image.select("B3")
b4 = image.select("B4")
turbidity = Turbidity(b3, b4)
Map.addLayer(turbidity.clip(taihu), visParam_t, date)
Map