太湖水域蓝藻密度计算
# 太湖水域蓝藻密度分析
基于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) 核心代码示例
/**
* @File : PIETaiHuAlg
* @Time : 2020/7/21
* @Author : piesat
* @Version : 1.0
* @Contact : 400-890-0662
* @License : (C)Copyright 航天宏图信息技术股份有限公司
* @Desc : 太湖水域蓝藻密度
*/
//输入太湖边界矢量
var taihu = pie.FeatureCollection('user/101/public/shape/TaiHuShuiYu')
.first()
.geometry();
Map.centerObject(taihu, 8);
//加载Landsat 8 SR
var l8Col = pie.ImageCollection("LC08/01/T1_SR")
.filterDate("2018-01-01", "2021-01-01")
.filterBounds(taihu)
.filter(pie.Filter.lte('cloudCover', 5))
.map(function(image) {
return image.set("id", image.id());
});
//设置图层显示属性
var visalg = {
min: 200,
max: 50000,
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 alg_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 b4 = image.select("B4");
var b5 = image.select("B5");
//计算蓝藻密度
var alg = ((((b5.divide(b4)).power(2)).multiply(1352)).subtract
((b5.divide(b4)).multiply(159.08))).add(192.87)
.clip(taihu);
//print("太湖水域蓝藻密度", date, alg);
Map.addLayer(alg, visalg, date, false);
var alg_mean = alg.reduceRegion(pie.Reducer.mean(), taihu, 30);
//print("太湖水域蓝藻密度平均值", date, alg_mean);
var alg_max = alg.reduceRegion(pie.Reducer.max(), taihu, 30);
//print("太湖水域蓝藻密度最大值", date, alg_max);
var alg_min = alg.reduceRegion(pie.Reducer.min(), taihu, 30);
//print("太湖水域蓝藻密度最小值", date, alg_min);
dates.push(date);
alg_images.push(alg_mean);
}
// 设置图表属性
var line_a = {title: '太湖水域蓝藻密度动态变化',
legend: ['蓝藻密度'],
xAxisName: "日期",
yAxisName: "蓝藻密度(万/L)",
chartType: "line",
yMin: 0,
yMax: 2500,
smooth: true};
// 显示折线图
ChartImage(alg_images, dates, line_a);
//动画显示
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: ["200(万个/L)", "2000(万个/L)"],
step: 30
};
var style = {
right: "150px",
bottom: "10px",
height: "70px",
width: "350px"
};
var legend = ui.Legend(data, style);
Map.addUI(legend);
})复制错误复制成功
具体的流程步骤如下:
a) 加载太湖边界矢量;
b) 筛选符合要求的覆盖太湖区域的Landsat8 SR数据集;
c) 设置图层显示属性;
d) 根据公式计算太湖水域蓝藻密度,及其平均值、最大值、最小值;
d) 设置折线图属性;
e) 动画显示太湖水域蓝藻密度变化;
f) 设置太湖水域蓝藻密度变化图例;
# 2) 案例中主要涉及到的算子
• pie.ImageCollection("xxx"):加载自己上传的影像数据集或者系统内置的影像数据集
• filterDate:对影像数据集进行时间过滤
• filterBounds:对影像数据集进行空间过滤
• reduceRegion:统计分析指定区域的信息结果
• reduceColumns :根据Reducer,对“date”, “id”属性值进行计算
• ChartArray:绘制图表
• playLayersAnimation:播放动画
• addUI:在地图上添加UI组件
# 3) 结果分析展示
使用PIE-Engine遥感计算云平台的实现效果如下图所示:
动态变化展示如下:
# 4) 完整代码
/**
* @File : PIETaiHuAlg
* @Time : 2020/7/21
* @Author : piesat
* @Version : 1.0
* @Contact : 400-890-0662
* @License : (C)Copyright 航天宏图信息技术股份有限公司
* @Desc : 太湖水域蓝藻密度
*/
//输入太湖边界矢量
var taihu = pie.FeatureCollection('user/101/public/shape/TaiHuShuiYu')
.first()
.geometry();
Map.centerObject(taihu, 8);
//加载Landsat 8 SR
var l8Col = pie.ImageCollection("LC08/01/T1_SR")
.filterDate("2018-01-01", "2021-01-01")
.filterBounds(taihu)
.filter(pie.Filter.lte('cloudCover', 5))
.map(function(image) {
return image.set("id", image.id());
});
//设置图层显示属性
var visalg = {
min: 200,
max: 50000,
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 alg_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 b4 = image.select("B4");
var b5 = image.select("B5");
//计算蓝藻密度
var alg = ((((b5.divide(b4)).power(2)).multiply(1352)).subtract
((b5.divide(b4)).multiply(159.08))).add(192.87)
.clip(taihu);
//print("太湖水域蓝藻密度", date, alg);
Map.addLayer(alg, visalg, date, false);
var alg_mean = alg.reduceRegion(pie.Reducer.mean(), taihu, 30);
//print("太湖水域蓝藻密度平均值", date, alg_mean);
var alg_max = alg.reduceRegion(pie.Reducer.max(), taihu, 30);
//print("太湖水域蓝藻密度最大值", date, alg_max);
var alg_min = alg.reduceRegion(pie.Reducer.min(), taihu, 30);
//print("太湖水域蓝藻密度最小值", date, alg_min);
dates.push(date);
alg_images.push(alg_mean);
}
// 设置图表属性
var line_a = {title: '太湖水域蓝藻密度动态变化',
legend: ['蓝藻密度'],
xAxisName: "日期",
yAxisName: "蓝藻密度(万/L)",
chartType: "line",
yMin: 0,
yMax: 2500,
smooth: true};
// 显示折线图
ChartImage(alg_images, dates, line_a);
//动画显示
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: ["200(万个/L)", "2000(万个/L)"],
step: 30
};
var style = {
right: "150px",
bottom: "10px",
height: "70px",
width: "350px"
};
var legend = ui.Legend(data, style);
Map.addUI(legend);
})