【GEE】基于PCA的LANDSAT 8计算遥感生态指数（RSEI）

var table = ee.FeatureCollection("projects/ee-1261423515/assets/aibihu");
Map.centerObject(table, 8);
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);

}

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;

}

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;

}

function NDVI_cal(img) {

var ndvi_temp = img.normalizedDifference(["SR_B5","SR_B4"]);

//print(ndvi_temp);

return ndvi_temp;

}

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;

}

function NDWI_cal(img) {

var nir = img.select("SR_B5");

var green = img.select("SR_B3");

var ndwi_temp = green.subtract(nir).divide(green.add(nir));

return ndwi_temp;

}

function img_normalize(img){

var minMax = img.reduceRegion({

reducer:ee.Reducer.minMax(),

geometry: table,

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());
    return normalize;

} ;

// Map the function over one year of data.

var dataset = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2")

.filterDate('2020-06-01', '2020-10-31')

.map(maskL8sr)

.mean()

var dataset1 = ee.ImageCollection("MODIS/006/MOD11A2")

.filterDate('2020-06-01', '2020-10-31')

.mean()

var ndwi=NDWI_cal(dataset);

var NDWI_mask = ndwi.lt(0.1);

var dataset_no_water = dataset.updateMask(NDWI_mask).clip(table);

var dataset1_no_water = dataset1.updateMask(NDWI_mask).clip(table);
var lst = dataset1_no_water.expression(

'a*0.02-273.15',

{

a:dataset1_no_water.select('LST_Day_1km'),

});

var SI = SI_cal(dataset_no_water);

var IBI = IBI_cal(dataset_no_water);

var NDBSI =(SI.add(IBI)).divide(2.0);

var ndvi = NDVI_cal(dataset_no_water);

var wet = Wet_cal(dataset_no_water);

var visParams1 = {

palette: '313695,4575b4,74add1,abd9e9,e0f3f8,ffffbf,fee090,fdae61,f46d43,d73027,a50026'

};

var visParams = {

min: -1,

max: 1,

palette: ['0000FF','ff0000','00ff00']

};

var visParams2 = {

palette: 'a50026,d73027,f46d43,fdae61,fee090,ffffbf,e0f3f8,abd9e9,74add1,4575b4,313695'

};

var unit_ndvi = img_normalize(ndvi);

dataset_no_water=dataset_no_water.addBands(unit_ndvi.rename('NDVI').toFloat())

var unit_NDBSI = img_normalize(NDBSI);

dataset_no_water=dataset_no_water.addBands(unit_NDBSI.rename('NDBSI').toFloat())

var unit_Wet = img_normalize(wet);

dataset_no_water=dataset_no_water.addBands(unit_Wet.rename('Wet').toFloat())

var unit_lst = img_normalize(lst);

dataset_no_water=dataset_no_water.addBands(unit_lst.rename('LST').toFloat())

//print(dataset_no_water);

var bands = ["Wet","NDVI","NDBSI","LST"]

var sentImage =dataset_no_water.select(bands)
var region = table;

var image =  sentImage.select(bands);

var scale = 1000;

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);
var rsei_un_unit = pcImage.expression(

'constant - pc1' ,

{

constant: 1,

pc1: pcImage.select('pc1')

});

var rsei = img_normalize(rsei_un_unit);
//Map.addLayer(table,{color:'FFFF00'});

Map.addLayer(rsei, {}, 'PCA');

//Map.addLayer(dataset_no_water, {bands: ['SR_B4', 'SR_B3', 'SR_B2'], min: 0, max: 0.3});

//Map.addLayer(lst, visParams1, "LST");

//Map.addLayer(NDBSI, {}, "NDBSI");

//Map.addLayer(ndvi, visParams, "NDVI");

//Map.addLayer(wet, visParams2, "Wet");