【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");

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232