GEE计算遥感生态指数(RSEI)--Landsat 8为例RSEI原理代码运行结果感谢

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
233
234
235
236
237
238
239
240
241
242
243

// roi和table不是一个geometry，运行的话，注意一下
var roi = table.geometry().bounds()

Map.centerObject(roi, 5);

function remove_water(img) {
    var year = img.get('year')
    // jrc_year: JRC Yearly Water Classification History, v1.3
    var Mask = jrc_year.filterMetadata('year', 'equals', year)
                       .filterBounds(roi)
                       .first()
                       .select('waterClass')
                       .reproject('EPSG:4326',null,1000)
                       // 掩膜掉大片永久水体
                       .eq(3)
    // 此时Mask中value值有1 ， 0 ， masked，把masked转化为 0
    Mask = Mask.unmask(0).not()
    return img.updateMask(Mask)
}

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')
           .filterBounds(roi)
           .filterDate('2014-01-01', '2019-12-31')
           .filterMetadata('CLOUD_COVER', 'less_than',50)
           .map(function(image){
                    return image.set('year', ee.Image(image).date().get('year'))                           
                  })
           .map(removeCloud)


var L8imgList = ee.List([])
for(var a = 2014; a < 2020; a++){
   var img = L8.filterMetadata('year', 'equals', a).median().clip(roi)
   var L8img = img.set('year', a)
   L8imgList = L8imgList.add(L8img)
 }

var L8imgCol = ee.ImageCollection(L8imgList)
                 .map(function(img){
                      return img.clip(roi)
                   })
                 .map(remove_water)
                 
L8imgCol = L8imgCol.map(function(img){
  
  // 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'])
     })   
  img = img.addBands(Wet.rename('WET'))
  
  
  // NDVI
  var ndvi = img.normalizedDifference(['B5', 'B4']);
  img = img.addBands(ndvi.rename('NDVI'))
  
  
  // lst 直接采用MODIS产品，表示看不懂 LST 公式
  var lst = ee.ImageCollection('MODIS/006/MOD11A1').map(function(img){
                return img.clip(roi)
           })
           .filterDate('2014-01-01', '2019-12-31')
  
  var year = img.get('year')
  lst=lst.filterDate(ee.String(year).cat('-01-01'),ee.String(year).cat('-12-31')).select(['LST_Day_1km', 'LST_Night_1km']);
      
  // reproject主要是为了确保分辨率为1000
  var img_mean=lst.mean().reproject('EPSG:4326',null,1000);
  //print(img_mean.projection().nominalScale())
  
  img_mean = img_mean.expression('((Day + Night) / 2)',{
      'Day': img_mean.select(['LST_Day_1km']),
      'Night': img_mean.select(['LST_Night_1km']),
       })
  img = img.addBands(img_mean.rename('LST'))
  
  
  // ndbsi = ( ibi + si ) / 2
  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)
  return img.addBands(ndbsi.rename('NDBSI'))
})


var bandNames = ['NDVI', "NDBSI", "WET", "LST"]
L8imgCol = L8imgCol.select(bandNames)

// 归一化
var img_normalize = function(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()).set('year', year);
        return normalize;
}
var imgNorcol  = L8imgCol.map(img_normalize);


// 主成分
var pca = function(img){
      
      var bandNames = img.bandNames();
      var region = roi;
      var year = img.get('year')
      // Mean center the data to enable a faster covariance reducer
      // and an SD stretch of the principal components.
      var meanDict = img.reduceRegion({
            reducer:  ee.Reducer.mean(),
            geometry: region,
            scale: 1000,
            maxPixels: 10e13
        });
      var means = ee.Image.constant(meanDict.values(bandNames));
      var centered = img.subtract(means).set('year', year);
      
      
      // This helper function returns a list of new band names.
      var getNewBandNames = function(prefix, bandNames){
            var seq = ee.List.sequence(1, 4);
            //var seq = ee.List.sequence(1, bandNames.length());
            return seq.map(function(n){
                  return ee.String(prefix).cat(ee.Number(n).int());
              });      
        };
      
      // This function accepts mean centered imagery, a scale and
      // a region in which to perform the analysis.  It returns the
      // Principal Components (PC) in the region as a new image.
      var getPrincipalComponents = function(centered, scale, region){
            var year = centered.get('year')
            var arrays = centered.toArray();
        
            // Compute the covariance of the bands within the region.
            var covar = arrays.reduceRegion({
                  reducer: ee.Reducer.centeredCovariance(),
                  geometry: region,
                  scale: scale,
                  bestEffort:true,
                  maxPixels: 10e13
              });
            
            // Get the 'array' covariance result and cast to an array.
            // This represents the band-to-band covariance within the region.
            var covarArray = ee.Array(covar.get('array'));
            
            // Perform an eigen analysis and slice apart the values and vectors.
            var eigens = covarArray.eigen();
        
            // This is a P-length vector of Eigenvalues.
            var eigenValues = eigens.slice(1, 0, 1);
            // This is a PxP matrix with eigenvectors in rows.
            var eigenVectors = eigens.slice(1, 1);
        
            // Convert the array image to 2D arrays for matrix computations.
            var arrayImage = arrays.toArray(1)
            // Left multiply the image array by the matrix of eigenvectors.
            var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
        
            // Turn the square roots of the Eigenvalues into a P-band image.
            var sdImage = ee.Image(eigenValues.sqrt())
            .arrayProject([0]).arrayFlatten([getNewBandNames('SD',bandNames)]);
        
            // Turn the PCs into a P-band image, normalized by SD.
            return principalComponents
            // Throw out an an unneeded dimension, [[]] -> [].
            .arrayProject([0])
            // Make the one band array image a multi-band image, [] -> image.
            .arrayFlatten([getNewBandNames('PC', bandNames)])
            // Normalize the PCs by their SDs.
            .divide(sdImage)
            .set('year', year);
        }
        
        // Get the PCs at the specified scale and in the specified region
        img = getPrincipalComponents(centered, 1000, region);
        return img;
  };
  
var PCA_imgcol = imgNorcol.map(pca)
 
Map.addLayer(PCA_imgcol.first(), {"bands":["PC1"]}, 'pc1')

// 利用PC1，计算RSEI，并归一化
var RSEI_imgcol = PCA_imgcol.map(function(img){
        img = img.addBands(ee.Image(1).rename('constant'))
        var rsei = img.expression('constant - pc1' , {
             constant: img.select('constant'),
             pc1: img.select('PC1')
         })
        rsei = img_normalize(rsei)
        return img.addBands(rsei.rename('rsei'))
    })
print(RSEI_imgcol)

var visParam = {
    palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +
        '3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
 };

Map.addLayer(RSEI_imgcol.first().select('rsei'), visParam, 'rsei')


Map.addLayer(table, null, 'river_buffer')

//var Foo = require('users/2210902141/RSEI:downloadModules.js');

//Foo.foo(rsei_imgcol.first(), roi);