今天来简单分享下如何在GEE中批量提取生物量组分指数(BCI)
1.BCI指数介绍：
缨帽变换是由Kauth提出的，用于识别农业谷物长势的一种方法。由于它具有压缩光谱信息和突出不同土地覆被类型光谱特征的优点，近年来被广泛应用。在城市遥感上，TC分量广泛用于机器学习方法估算不透水层，该指数与不透水层丰度呈现明显正相关。不同的卫星传感器的缨帽变换系数各不相同，变换参数可查阅相关的文献资料。
BCI指数的思想来源于Ridd提出的V-I-S模型（Vegetation-Imperious-Soil），该模型将城市地表看作由植被、不透水面和土壤三种基本组分组成，是一种常用于城市环境中区分不同地物类型的指数。该指数是基于缨帽变换（TCtransformation）生成的亮度（TC1）、绿度（TC2）和湿度（TC3）三个分量所构建的。在城市遥感上，TC分量广泛用于机器学习方法估算不透水层。该指数与不透水层丰度呈现明显正相关，而与植被丰度呈负相关。通过BCI指数运算，所有地物的数值都分布在-1到1间。
式中H为高反射率即归一化TC1分量；V为植被即归一化TC2分量；L为低反射率即归一化TC3分量。
三个因子的计算公式如下：
其中TCi(i=1,2,3)是前三个TC分量；TCimin和TCimax分别是第i个TC分量的最小值和最大值。
计算水体掩膜的方式选择了徐涵秋提出的调整的归一化水体指数（Modification of normalized difference water index, MNDWI），该指数在以建筑用地为背景的水体提取上，较归一化水体指数更有优势。
其中，Green和NIR分别对应着遥感影像中的绿波段和近红外波段。
计算BI指数，选择MNDWI阈值，去除大于阈值的水体。BI指数利用土壤在红波段和短波红外波段的反射率较其它地表覆被类型高，而在中红外波段和蓝波段较其它地表覆被类型偏低的特征，可以有效提取土壤信息。它的取值在0到200之间，BI值越大，则土壤的概率越大，这样就可以通过设置阈值来提取土壤信息。BI值计算公式如下：
其中，BLUE、RED、NIR、SWIR1分别对应着遥感影像中的蓝、红、近红外和短波红外波段。
BCI指数能够更好地区分土壤、水域和植被。
2.GEE实现代码：
首先确定研究区和使用的数据集
我选择的研究区为徐州市，使用的数据集为LANDSAT/LC08/C02/T1_L2
目标是实现徐州市2020到2023年的BCI指数提取
以下是实现代码:
//将研究区缩放至合适大小
Map.centerObject(roi,7)
var styling = {color:"red",fillColor:"00000000"};
Map.addLayer(roi.style(styling),{},"geometry")
//landsat8的去云函数
function maskL8sr(image) {
// Bit 0 - Fill
// Bit 1 - Dilated Cloud
// Bit 2 - Cirrus
// Bit 3 - Cloud
// Bit 4 - Cloud Shadow
var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);
var saturationMask = image.select('QA_RADSAT').eq(0);
// Apply the scaling factors to the appropriate bands.
var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
// Replace the original bands with the scaled ones and apply the masks.
return image.addBands(opticalBands, null, true)
.addBands(thermalBands, null, true)
.updateMask(qaMask)
.updateMask(saturationMask);
}
//计算landsat8的NDVI
var getL8NDVI=function(image){
var ndvi = image.expression(
'(B5-B4)/(B5+B4)',{
'B5':image.select('SR_B5'),
'B4':image.select('SR_B4')
})
.clamp(-1,1)//去除异常值
.rename('ndvi')
.set('system:time_start', image.get('system:time_start'));
return (ndvi);
}
//利用mndwi指数将水体进行掩膜
var MNDWI_mask=function(image){
var water_mask = image.expression(
'(B3-B6)/(B3+B6)',{
'B3':image.select('SR_B3'),
'B6':image.select('SR_B6')
}).rename('mndwi')
var water=water_mask.lte(0).rename('mndwi'); //mndwi>0为水体
return (image.updateMask(water));
}
//缨帽变换
var getL8BCI=function(tif){
//
var tcb = tif.expression(
"B2*(0.3029)+B3*(0.2786)+B4*(0.4733)+B5*(0.5599)+B6*(0.5082)+B7*(0.1872)",
{
"B2": tif.select(["SR_B2"]),
"B3": tif.select(["SR_B3"]),
"B4": tif.select(["SR_B4"]),
"B5": tif.select(["SR_B5"]),
"B6": tif.select(["SR_B6"]),
"B7": tif.select(["SR_B7"]),
}).rename("TCB")
var tcg = tif.expression(
"B2*(-0.2941)+B3*(-0.2435)+B4*(-0.5424)+B5*(0.7276)+B6*(0.0713)+B7*(-0.1608)",
{
"B2": tif.select(["SR_B2"]),
"B3": tif.select(["SR_B3"]),
"B4": tif.select(["SR_B4"]),
"B5": tif.select(["SR_B5"]),
"B6": tif.select(["SR_B6"]),
"B7": tif.select(["SR_B7"]),
}).rename("TCG")
var tcw = tif.expression(
"B2*(0.1511)+B3*(0.1973)+B4*(0.3283)+B5*(0.3406)+B6*(-0.7112)+B7*(-0.4559)",
{
"B2": tif.select(["SR_B2"]),
"B3": tif.select(["SR_B3"]),
"B4": tif.select(["SR_B4"]),
"B5": tif.select(["SR_B5"]),
"B6": tif.select(["SR_B6"]),
"B7": tif.select(["SR_B7"]),
}).rename("TCW")
var bci=(((tcb.add(tcw)).divide(2)).subtract(tcg)).divide(((tcb.add(tcw)).divide(2)).add(tcg))
return(bci.rename('BCI'))
}
//Get annual composite images
var start_year = 2020;
var end_year = 2023;
var yearList = ee.List.sequence(start_year, end_year);
var yearImgs = yearList.map(function(year) {
year = ee.Number(year).toInt();
var _sdate = ee.Date.fromYMD(year, 1, 1);
var _edate = ee.Date.fromYMD(year.add(1), 1, 1);
var tempCol = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.filterDate(_sdate, _edate)
.filterBounds(roi)
.map(maskL8sr)
.map(function(image){
return(image.clip(roi).clamp(-1,1))
})
var ndvi_mask=tempCol.map(getL8NDVI).reduce(ee.Reducer.variance()).lte(0.055);
var yearImg = tempCol.median().updateMask(ndvi_mask);
//var yearImg = tempCol.median();
yearImg = yearImg.set("count", tempCol.size());
yearImg = yearImg.set("year", year);
yearImg = yearImg.set("system:index", ee.String(year));
return yearImg ;
});
var yearImgCol = ee.ImageCollection.fromImages(yearImgs)
.filter(ee.Filter.gt("count", 0));
print("yearImgCol", yearImgCol);
var l8BCI = yearImgCol.map(MNDWI_mask)
.map(getL8BCI)
print("l8BCI", l8BCI);
Map.addLayer(ee.Image(l8BCI.toList(4).get(0)),{min: 0, max: 1, palette: ['0000FF','FF0000']},'BCI_2023')
Map.addLayer(ee.Image(l8BCI.toList(1).get(0)),{min: 0, max: 1, palette: ['0000FF','FF0000']},'BCI_2020')
Map.addLayer(ee.Image(l8BCI.toList(2).get(0)),{min: 0, max: 1, palette: ['0000FF','FF0000']},'BCI_2021')
Map.addLayer(ee.Image(l8BCI.toList(3).get(0)),{min: 0, max: 1, palette: ['0000FF','FF0000']},'BCI_2022')
///----------------------数据批量输出函数-----------------------//
function exportImageCollection(imgCol, scale, roi, taskName, fileName) {
var indexList = l8BCI.reduceColumns(ee.Reducer.toList(), ["system:index"]).get("list");
indexList.evaluate(function (indexs) {
for (var i = 0; i < indexs.length; i++) {
var image = imgCol.filter(ee.Filter.eq("system:index", indexs[i])).first();
image = image.float(); //
Export.image.toDrive({
image: image,
description: taskName + "_" + indexs[i],
fileNamePrefix: fileName + "_" + indexs[i],
region: roi,
folder: 'BCI',
scale: scale,
crs: "EPSG:4326",
maxPixels: 1e13
});
}
});
}
exportImageCollection(l8BCI, 30, roi, "BCI", "BCI");结果显示：
以下图中红色为不透水面的主要分布地区
2020
2021
2022
2023
点击run all.全部导出至云盘
完整代码请在公众号后台私信“10.31批量计算BCI”
感谢关注，欢迎转发！
声明：仅供学习使用！
