首页 > 编程札记 > 编程

GEE计算遥感生态指数(RSEI)

阅读：评论：0

GEE计算遥感生态指数(RSEI)

RSEI原理

湿度指标（Wet）

绿度指标（NDVI）

热度指标（LST）

干度指标（NDBSI）

Landsat-8波段

归一化（normalization）

主成分分析( PCA )

PCA-->RSEI

Wet、NDVI、LST、NDBSI、PCA、RSEI之间的关系

代码

运行结果

感谢

RSEI原理

RSEI是一个完全基于遥感技术，以自然因子为主的遥感生态指数(RSEI) 来对城市的生态状况进行快速监测与评价的方法，由徐涵秋.城市遥感生态指数的创建及其应用[J].生态学报,2013,33(24):7853-7862.中提出。

该指数利用主成分分析技术集成了植被指数、湿度分量、地表温度和建筑指数等 4 个评价指标，它们分别代表了绿度、湿度、热度和干度等 4 大生态要素。

湿度指标（Wet）

WET = c1 B1 + c2 B2 + c3 B3 + c4 B4 + c5 B5 + c6 B6

B1-B6 分别代表蓝波段、绿波段、红波段、近红波段、中红外波段 1、中红外波段 2；

c1~c6是传感器参数。由于传感器的类型不同，参数也相应有所不同。

其中，TM 传感器，c1 ～ c6 分别为 0.0315、0.2021、0.3012、0.1594、-0.6806、-0.6109；

OLI 传感器，c1 ～ c6 分别为 0.1511、 0.1973、 0.3283、 0.3407、-0.7117、 -0.4559。

Landsat 8上携带的是陆地成像仪（Operational Land Imager ，OLI）

绿度指标（NDVI）

NDVI=（NIR-Red）/（NIR+Red）

热度指标（LST）

直接计算有点麻烦，故直接采用 MODIS 的LST，利用以下公式

LST = （ LST_Day_1km + LST_Night_1km ）/ 2

补充：现在才知道 Landsat 也有LST产品，且分辨率为30m，但是不改代码了，如何分辨率改为30m，下面代码的研究区太大了，会有点问题

LANDSAT/LC08/C02/T1_L2

介绍：

://le/earth-engine/datasets/catalog/LANDSAT_LC08_C02_T1_L2

干度指标（NDBSI）

IBI 是建筑指数

SI为土壤指数

Landsat-8波段

徐涵秋使用的是Landsat-7 ETM+遥感影像，但是本文采用的是Landsat-8遥感影像为例

归一化（normalization）

先使用ee.Reducer.minMax求出该幅遥感影像的最大值和最小值

在使用unitScale(low, high)使输入值的范围[低，高]变为[0，1]

主成分分析( PCA )

PCA就是一个有损的降维算法。

具体的请看代码的英文注释，不敢讲，请自行百度理解。

PCA-->RSEI

Wet、NDVI、LST、NDBSI、PCA、RSEI之间的关系

在 PC1 中，代表绿度的 NDVI 和代表湿度的 Wet 指标呈正值，说明它们对生态系统起正面的贡献，而代表热度和干度的 LST、NDBSI 则呈负值，这与实际情况相符。

综合来看，以 NDVI 为代表的植被和以 NDBSI 为代表的建筑用地对城市生态的影响力最大，且 NDVI 大于 NDBSI。NDBSI 所代表的建筑不透水面对城市地表温度LST具有正相关关系

RSEI 即为所建的遥感生态指数，其值介于［0，1］之间。ＲSEI 值越接近 1，生态越好，反之，生态越差。

好好品一下这个话：为使 PC1 大的数值代表好的生态条件，可进一步用 1 减去 PC1，获得初始的生态指数 RSEI0

PC1的值大，不代表生态好，PC1小的值，也不代表生态不好。只有RSEI0 = 1 - PC1才能代表生态

代码


// roi和table不是一个geometry，运行的话，注意一下
var roi = ry().bounds()Object(roi, 5);function remove_water(img) {var year = ('year')// jrc_year: JRC Yearly Water Classification History, v1.3var 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转化为 0Mask = 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){// Wetvar Wet = pression('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.ame('WET'))// NDVIvar ndvi = alizedDifference(['B5', 'B4']);img = img.ame('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 = ('year')lst=lst.filterDate(ee.String(year).cat('-01-01'),ee.String(year).cat('-12-31')).select(['LST_Day_1km', 'LST_Night_1km']);// reproject主要是为了确保分辨率为1000var img_mean&#an().reproject('EPSG:4326',null,1000);//print(img_mean.projection().nominalScale())img_mean = pression('((Day + Night) / 2)',{'Day': img_mean.select(['LST_Day_1km']),'Night': img_mean.select(['LST_Night_1km']),})img = img.addBands(ame('LST'))// ndbsi = ( ibi + si ) / 2var ibi = pression('(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 = pression('((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.ame('NDBSI'))
})var bandNames = ['NDVI', "NDBSI", "WET", "LST"]
L8imgCol = L8imgCol.select(bandNames)// 归一化
var img_normalize = function(img){var minMax = duceRegion({reducer:ee.Reducer.minMax(),geometry: roi,scale: 1000,maxPixels: 10e13,})var year = ('year')var normalize  = ee.ImageCollection.fromImages(img.bandNames().map(function(name){name = ee.String(name);var band = img.select(name);return band.unitScale(ee.(name.cat('_min'))), ee.(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 = ('year')// Mean center the data to enable a faster covariance reducer// and an SD stretch of the principal components.var meanDict = duceRegion({reducer:  an(),geometry: region,scale: 1000,maxPixels: 10e13});var means = stant(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 = ('year')var arrays = Array();// Compute the covariance of the bands within the region.var covar = duceRegion({reducer: edCovariance(),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'));// 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 = Array(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 urn 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 regionimg = 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 = pression('constant - pc1' , {constant: img.select('constant'),pc1: img.select('PC1')})rsei = img_normalize(rsei)return img.ame('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);