本期也是GEE的时间，细心的朋友会发现，开了赞赏功能，每天都是干货，还不赏我一瓶啤酒？那么，本期分享如何用GEE基于Landsat 8数据反演绿度/热度/湿度/干度，并计算生态遥感指数，代码较长，如有不妥之处，后台私信即可。
代码语言：javascript复制//导入自己的研究区，将其定义为roi
var star_date = '2018-01-01'//定义起始时间
var end_date = '2018-12-30'//定义终止时间
var L8_SR = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")//加载L8_SR影像
var img_SR = L8_SR.filterBounds(roi)
.filterDate(star_date, end_date)
.filterMetadata('CLOUD_COVER', 'less_than',50)
.mean()//按照条件进行筛选
//归一化函数
var img_normalize = function(img){
var minMax = img.reduceRegion({
reducer:ee.Reducer.minMax(),
geometry: roi,
scale: 100,
maxPixels: 10e13,
tileScale: 16
})
var year = img.get('year')
var normalize = ee.ImageCollection.fromImages(
img.bandNames().map(function(name){
name = ee.String(name);
var band = img.select(name);
return band.unitScale(ee.Number(minMax.get(name.cat('_min'))), ee.Number(minMax.get(name.cat('_max'))));
})
).toBands().rename(img.bandNames());
return normalize;
}
// Load the image of landsat 8 data
var vizParams = {bands: ['B5', 'B4', 'B3'], min: 50, max: 1500, gamma: [0.95,1.1,1]};
var images = img_SR.clip(roi);
var thermal = images.select('B10').multiply(0.1)
var Emissivity = function(image, nir, red, min, max, de) {
min = min || 0.2
max = max || 0.5
de = de || 0.005
var ndvi_e = "(nir-red)/(nir red)"
var ndvi = ee.Image(0).expression(ndvi_e, {'nir': image.select(nir), 'red': image.select(red)}).rename('ndvi')
var pv = ee.Image(0).expression("((ndvi-min)/(max-min))**2", {ndvi: ndvi, min: min, max: max})
var exp = "ndvi < 0.2 ? 0.979 - (0.046 * b4) : 0.2 <= ndvi <= 0.5 ? (0.987 * pv) 0.971 * (1 - pv) de : 0.987 de"
return ee.Image(0).expression(exp, {ndvi: ndvi, b4: image.select(red).multiply(0.0001), pv: pv, de: de}).rename('EM')
}
var EMM1 = Emissivity(img_SR, 'B5', 'B4')
var EMM = EMM1.clip(roi);
//LST计算
var LST = thermal.expression(
'(Tb/(1 (0.0010904 * (Tb / 1.438))*log(Ep)))-273.15', {
'Tb': thermal.select('B10'),
'Ep': EMM.select('EM'),
});
LST = img_normalize(LST)
img_SR = img_SR.addBands(LST.rename('LST').toFloat())
var SR_LST= img_SR.select('LST')
// 去云函数
function removeCloud(image){
var qa = image.select('BQA')
var cloudMask = qa.bitwiseAnd(1 << 4).eq(0)
var cloudShadowMask = qa.bitwiseAnd(1 << 8).eq(0)
var valid = cloudMask.and(cloudShadowMask)
return image.updateMask(valid)
}
var L8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
var img = L8.filterBounds(roi)
.filterDate(star_date, end_date)
.filterMetadata('CLOUD_COVER', 'less_than',50)
.map(removeCloud).mean()
//计算NDVI
var ndvi = img.normalizedDifference(['B5', 'B4']);
img = img.addBands(ndvi.rename('NDVI'))
//计算WET
//WET
var Wet = img.expression('B*(0.1509) G*(0.1973) R*(0.3279) NIR*(0.3406) SWIR1*(-0.7112) SWIR2*(-0.4572)',{
'B': img.select(['B2']),
'G': img.select(['B3']),
'R': img.select(['B4']),
'NIR': img.select(['B5']),
'SWIR1': img.select(['B6']),
'SWIR2': img.select(['B7'])
})
Wet = img_normalize(Wet)
img = img.addBands(Wet.rename('WET').toFloat())
//计算NDBSI
var ibi = img.expression('(2 * SWIR1 / (SWIR1 NIR) - (NIR / (NIR RED) GREEN / (GREEN SWIR1))) / (2 * SWIR1 / (SWIR1 NIR) (NIR / (NIR RED) GREEN / (GREEN SWIR1)))', {
'SWIR1': img.select('B6'),
'NIR': img.select('B5'),
'RED': img.select('B4'),
'GREEN': img.select('B3')
})
var si = img.expression('((SWIR1 RED) - (NIR BLUE)) / ((SWIR1 RED) (NIR BLUE))', {
'SWIR1': img.select('B6'),
'NIR': img.select('B5'),
'RED': img.select('B4'),
'BLUE': img.select('B2')
})
var ndbsi = (ibi.add(si)).divide(2)
ndbsi= img_normalize(ndbsi)
img = img.addBands(ndbsi.rename('NDBSI'))
var L8_img = img.select(["NDVI","WET","NDBSI"])
L8_img = SR_LST.addBands(L8_img)
//掩膜，有些研究区可能需要使用到
// var mask = L8_img.select('NDVI').gte(0)
// var L8img = L8_img.updateMask(mask)
var bands = ["WET","NDVI","NDBSI","LST"]
var sentImage =L8_img.select(bands)
//输入基本信息
var region = roi;
var image = sentImage.select(bands);
var scale = 100;
var bandNames = image.bandNames();
// 图像波段重命名函数
var getNewBandNames = function(prefix) {
var seq = ee.List.sequence(1, bandNames.length());
return seq.map(function(b) {
return ee.String(prefix).cat(ee.Number(b).int());
});
};
//数据平均
var meanDict = image.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: region,
scale: scale,
maxPixels: 1e9
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image.subtract(means);
//主成分分析函数
var getPrincipalComponents = function(centered, scale, region) {
// 图像转为一维数组
var arrays = centered.toArray();
// 计算相关系数矩阵
var covar = arrays.reduceRegion({
reducer: ee.Reducer.centeredCovariance(),
geometry: region,
scale: scale,
maxPixels: 1e9
});
// 获取“数组”协方差结果并转换为数组。
// 波段与波段之间的协方差
var covarArray = ee.Array(covar.get('array'));
// 执行特征分析，并分割值和向量。
var eigens = covarArray.eigen();
// 特征值的P向量长度
var eigenValues = eigens.slice(1, 0, 1);
//计算主成分载荷
var eigenValuesList = eigenValues.toList().flatten()
var total = eigenValuesList.reduce(ee.Reducer.sum())
var percentageVariance = eigenValuesList.map(function(item) {
return (ee.Number(item).divide(total)).multiply(100).format('%.2f')
})
print(eigenValues ,'特征值')
print("贡献率", percentageVariance)
// PxP矩阵，其特征向量为行。
var eigenVectors = eigens.slice(1, 1);
// 将图像转换为二维阵列
var arrayImage = arrays.toArray(1);
//使用特征向量矩阵左乘图像阵列
var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
// 将特征值的平方根转换为P波段图像。
var sdImage = ee.Image(eigenValues.sqrt())
.arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
//将PC转换为P波段图像，通过SD标准化。
principalComponents=principalComponents
// 抛出一个不需要的维度，[[]]->[]。
.arrayProject([0])
// 使单波段阵列映像成为多波段映像，[]->image。
.arrayFlatten([getNewBandNames('pc')])
// 通过SDs使PC正常化。
.divide(sdImage);
return principalComponents
};
//进行主成分分析，获得分析结果
var pcImage = getPrincipalComponents(centered, scale, region);
//RESI
var img = img.addBands(ee.Image(1).rename('constant'))
var rsei = img.expression('constant - pc1' , {
constant: img.select('constant'),
pc1: pcImage.select('pc1')
})
rsei = img_normalize(resi)
rsei = pcImage.addBands(resi.rename('rsei'))
print(rsei,"RSEI")
var visParam = {
palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,'
'3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
};
//添加图层
Map.addLayer(rsei.select('rsei').clip(roi), visParam, 'rsei')
//计算RSEI平均值
var rsei_mean = image.reduceRegion({
reducer: ee.Reducer.mean(),
geometry:roi,
scale: 30,
maxPixels: 1e13
});
print(rsei_mean,'平均值')
var rsei_std = image.reduceRegion({
reducer: ee.Reducer.stdDev(),
geometry:roi,
scale: 30,
maxPixels: 1e13
});
print(resi_std,'标准差')
// var resimax = resi.select('rsei').reduceRegion({
// reducer: ee.Reducer.max(),
// geometry:roi,
// scale: 30,
// maxPixels: 1e13
// });
//下载函数
Export.image.toDrive({
// image: rsei.select(["rsei"]),
// description: 'rsei',
// scale: 30,
// region:roi,
// fileFormat: 'GeoTIFF',
// });
/////////////////下载掩膜过后的RSEI图像，只需要在下函数适用unmask(-1), 就能将掩膜掉的部分给转换为-1
// Export.image.toDrive({
// image: L8_img.select(["LST"]),
// description: 'LST',
// scale: 30,
// region:roi,
// fileFormat: 'GeoTIFF',
// });
本期“GEE基于Landsat 8数据反演绿度/热度/湿度/干度，并计算生态遥感指数代码”分享结束，感谢您的阅读！
