RSEI计算

https://zhuanlan.zhihu.com/p/664572483
https://www.cnblogs.com/jokerYk/p/15686064.html

//修改roi
var roi = table//center_group_clip;
Map.centerObject(roi, 8);
//Map.addLayer(roi,{},‘shp’);
//使用JRC/GSW1_4/YearlyHistory的30m分辨率数据进行水体掩膜，根据具体要求进行掩膜选择
//水体掩膜
function remove_water(img) {
var year = img.get(‘year’)
var jrc_year = ee.ImageCollection(‘JRC/GSW1_4/YearlyHistory’)
.filterDate(‘2020-01-01’, ‘2020-12-31’)
.first()
.clip(roi)
.select(‘waterClass’)
.reproject(‘EPSG:4326’,null,30)
// var Mask = jrc_year
// // 掩膜掉大片永久水体
// .eq(3)
var Mask = jrc_year
// 将永久水体和季节性水体都设为 1
.eq(3).or(jrc_year.eq(2))
// 此时Mask中value值有1 ， 0 ， masked，把masked转化为 0
Mask = Mask.unmask(0).not()
return img.updateMask(Mask)
}
//定义去云函数

//去云
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);
}
//SI
function SI_cal(img) {
var blue = img.select(“SR_B2”);
var red = img.select(“SR_B4”);
var nir = img.select(“SR_B5”);
var swir1 = img.select(“SR_B6”);
var swir2 = img.select(“SR_B7”);
var SI_temp =((swir1.add(red)).subtract(blue.add(nir)))
.divide((swir1.add(red)).add(blue.add(nir)))
//print(SI_temp);
return SI_temp;
}

//IBI
function IBI_cal(img) {
var green = img.select(“SR_B3”);
var red = img.select(“SR_B4”);
var nir = img.select(“SR_B5”);
var swir1 = img.select(“SR_B6”);
var IBI_temp =(((swir1.multiply(2.0)).divide(swir1.add(nir)))
.subtract((nir.divide(nir.add(red))).add(green.divide(green.add(swir1)))))
.divide((((swir1.multiply(2.0)).divide(swir1.add(nir)))
.add((nir.divide(nir.add(red))).add(green.divide(green.add(swir1))))))
//print(IBI_temp);
return IBI_temp;
}

//NDVI
function NDVI_cal(img) {
var ndvi_temp = img.normalizedDifference([“SR_B5”,“SR_B4”]);
//print(ndvi_temp);
return ndvi_temp;
}
//Wet
function Wet_cal(img) {
var blue = img.select(“SR_B2”);
var green = img.select(“SR_B3”);
var red = img.select(“SR_B4”);
var nir = img.select(“SR_B5”);
var swir1 = img.select(“SR_B6”);
var swir2 = img.select(“SR_B7”);
var wet_temp =blue.multiply(0.1511)
.add(green.multiply(0.1973))
.add(red.multiply(0.3283))
.add(nir.multiply(0.3407))
.add(swir1.multiply(-0.7117))
.add(swir2.multiply(-0.4559))
//print(wet_temp);
return wet_temp;
}

//归一化函数的定义

//img_normalize
function img_normalize(img){
var minMax = img.reduceRegion({
reducer:ee.Reducer.minMax(),
geometry: roi,
scale: 1000,
maxPixels: 10e13,
})
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());
    var normalized = normalize.min(1).max(0);
    return normalized;

}
//数据输入，Landsat以及MODIS数据，手动修改年份
//import data
var L8 = ee.ImageCollection(“LANDSAT/LC08/C02/T1_L2”)
.filterDate(‘2020-01-01’, ‘2020-12-31’)
.filterBounds(roi)
.filterMetadata(‘CLOUD_COVER’, ‘less_than’,20)
.map(maskL8sr)
.map(remove_water)
.mean()

//Map.addLayer(L8, {bands: [‘SR_B4’, ‘SR_B3’, ‘SR_B2’], min: 0, max: 0.3}, ‘L8_remove_water’)

var MODIS = ee.ImageCollection(“MODIS/006/MOD11A2”)
.filterDate(‘2020-04-01’, ‘2020-09-30’)
.map(remove_water)
.mean()

//镶嵌数据，水体掩膜
var L8_no_water = L8.clip(roi);
var MODIS = MODIS.clip(roi);
Map.addLayer(L8_no_water, {bands: [‘SR_B4’, ‘SR_B3’, ‘SR_B2’], min: 0, max: 0.3}, ‘L8_no_water_clip’)

//指数计算，带入函数
//index
var ndvi = NDVI_cal(L8_no_water);
var wet = Wet_cal(L8_no_water);
var lst = MODIS.expression(
‘a*0.02-273.15’,
{
a:MODIS.select(‘LST_Day_1km’),
});
var SI = SI_cal(L8_no_water);
var IBI = IBI_cal(L8_no_water);
var NDBSI =(SI.add(IBI)).divide(2.0);

//指数归一化，波段添加
//normalize index & add bands
var re_ndvi = img_normalize(ndvi);
L8_no_water = L8_no_water.addBands(re_ndvi.rename(‘NDVI’).toFloat())
var re_NDBSI = img_normalize(NDBSI);
L8_no_water = L8_no_water.addBands(re_NDBSI.rename(‘NDBSI’).toFloat())
var re_Wet = img_normalize(wet);
L8_no_water = L8_no_water.addBands(re_Wet.rename(‘Wet’).toFloat())
var re_lst = img_normalize(lst);
L8_no_water = L8_no_water.addBands(re_lst.rename(‘LST’).toFloat())
//print(L8_no_water);
// 导出具体的每个指数，供自己选择

//导出数据
// Export.image.toDrive({
// image: re_ndvi,
// description: ‘NDVI_2020’,
// scale: 500,
// region: roi,
// crs: ‘EPSG:4326’,
// maxPixels: 1e13,
// fileFormat: ‘GeoTIFF’,
// });

// Export.image.toDrive({
// image: re_Wet,
// description: ‘Wet_2020’,
// scale: 500,
// region: roi,
// crs: ‘EPSG:4326’,
// maxPixels: 1e13,
// fileFormat: ‘GeoTIFF’,
// });

// Export.image.toDrive({
// image: re_lst,
// description: ‘LST_2020’,
// scale: 500,
// region: roi,
// crs: ‘EPSG:4326’,
// maxPixels: 1e13,
// fileFormat: ‘GeoTIFF’,
// });

// Export.image.toDrive({
// image: re_NDBSI,
// description: ‘NDBSI_2020’,
// scale: 500,
// region: roi,
// crs: ‘EPSG:4326’,
// maxPixels: 1e13,
// fileFormat: ‘GeoTIFF’,
// });

//主成分分析前的预准备
//select bands & bandnames
var bands = [“NDVI”,“Wet”,“LST”,“NDBSI”]
var select_Image =L8_no_water.select(bands)
var scale = 1000;
var bandNames = select_Image.bandNames();
//print(select_Image)
//print(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());
});
};

//mean img-mean
var meanDict = select_Image.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: roi,
scale: scale,
maxPixels: 1e9
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = select_Image.subtract(means);

//主成分分析函数
var getPrincipalComponents = function(centered, scale, roi) {
// 图像转为一维数组
var arrays = centered.toArray();

// 计算协方差矩阵
var covar = arrays.reduceRegion({
  reducer: ee.Reducer.centeredCovariance(),
  geometry: roi,
  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);
print("特征向量", eigenVectors)  

// 将图像转换为二维阵列
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 pca_Image = getPrincipalComponents(centered, scale, roi);
//print(pca_Image)
//Map.addLayer(pca_Image, {“bands”:[“pc1”]}, ‘pc1’)

//根据特征向量，选择RSEI_0，根据第一主成分湿度的特征向量的正负
//1-pc1 r(wet)<0
var rsei_0 = pca_Image.expression(
‘constant - pc1’ ,
{
constant: 1,
pc1: pca_Image.select(‘pc1’)
});
//pc1 r(wet)>0
//var rsei_0 = pca_Image.select(‘pc1’)
//print(rsei_0)
//输出结果啦！

//RSEI归一化
var rsei = img_normalize(rsei_0)
var visParam = {
palette: ‘FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,’ +
‘3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301’
};
Map.addLayer(rsei, visParam, ‘rsei’);

Export.image.toDrive({
image: rsei,
description: ‘RESI_2020’,
scale: 30,
region: roi,
crs: ‘EPSG:4326’,
maxPixels: 1e13,
fileFormat: ‘GeoTIFF’,
});