前言

城市生态与人类生活息息相关，快速、准确、客观地了解城市生态状况已成为生态领域的一个研究重点。基于遥感技术，提出一个完全基于遥感技术，以自然因子为主的遥感生态指数 (RSEI)来对城市的生态状况进行快速监测与评价。该指数利用主成分分析技术集成了植被指数、湿度分量、地表温度和建筑指数等 4个评价指标，它们分别代表了绿度、湿度、热度和干度等4大生态要素。
本文基于GEE平台，实现RSEI算法。
运行结果：

第一步：定义研究区，自行更换自己的研究区

代码如下（示例）：

// 第一步:定义研究区，自行更换自己的研究区
var roi = /* color: #98ff00 *//* displayProperties: [{"type": "rectangle"}] */ee.Geometry.Polygon([[[120.1210075098537, 35.975189051414006],[120.1210075098537, 35.886229778229115],[120.25764996590839, 35.886229778229115],[120.25764996590839, 35.975189051414006]]], null, false);Map.centerObject(roi);

第二步：加载数据集,定义去云函数

代码如下（示例）：

// 第二步：加载数据集,定义去云函数
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('2018-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 = 2018; 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)})L8imgCol = L8imgCol.map(function(img){// 湿度函数:Wetvar 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'))// 绿度函数:NDVIvar ndvi = img.normalizedDifference(['B5', 'B4']);img = img.addBands(ndvi.rename('NDVI'))// 热度函数:lst 直接采用MODIS产品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主要是为了确保分辨率为1000var 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 ) / 2var 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);

第四步：PCA融合，提取第一主成分

代码如下（示例）：

// 第四步:PCA融合，提取第一主成分
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 regionimg = getPrincipalComponents(centered, 1000, region);return img;};var PCA_imgcol = imgNorcol.map(pca)Map.addLayer(PCA_imgcol.first(), {"bands":["PC1"]}, 'pc1')

第五步：利用PC1，计算RSEI，并归一化

代码如下（示例）：

// 第五步：利用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: '040274, 040281, 0502a3, 0502b8, 0502ce, 0502e6, 0602ff, 235cb1, 307ef3, 269db1, 30c8e2, 32d3ef, 3be285, 3ff38f, 86e26f, 3ae237, b5e22e, d6e21f, fff705, ffd611, ffb613, ff8b13, ff6e08, ff500d, ff0000, de0101, c21301, a71001, 911003'};Map.addLayer(RSEI_imgcol.first().select('rsei'), visParam, 'rsei')

完整代码

代码如下（示例）：


// 第一步:定义研究区，自行更换自己的研究区
var roi = /* color: #98ff00 *//* displayProperties: [{"type": "rectangle"}] */ee.Geometry.Polygon([[[120.1210075098537, 35.975189051414006],[120.1210075098537, 35.886229778229115],[120.25764996590839, 35.886229778229115],[120.25764996590839, 35.975189051414006]]], null, false);Map.centerObject(roi);// 第二步：加载数据集,定义去云函数
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('2018-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 = 2018; 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)})L8imgCol = L8imgCol.map(function(img){// 湿度函数:Wetvar 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'))// 绿度函数:NDVIvar ndvi = img.normalizedDifference(['B5', 'B4']);img = img.addBands(ndvi.rename('NDVI'))// 热度函数:lst 直接采用MODIS产品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主要是为了确保分辨率为1000var 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 ) / 2var 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);// 第四步:PCA融合，提取第一主成分
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 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 = 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: '040274, 040281, 0502a3, 0502b8, 0502ce, 0502e6, 0602ff, 235cb1, 307ef3, 269db1, 30c8e2, 32d3ef, 3be285, 3ff38f, 86e26f, 3ae237, b5e22e, d6e21f, fff705, ffd611, ffb613, ff8b13, ff6e08, ff500d, ff0000, de0101, c21301, a71001, 911003'};Map.addLayer(RSEI_imgcol.first().select('rsei'), visParam, 'rsei')

总结

提示：用完代码记得反馈。

利用GEE计算城市遥感生态指数(RSEI)——Landsat 8为例相关推荐

ENVI基于Landsat影像构建郑州市2000-2019年遥感生态指数RSEI
一.数据介绍由于时间跨度较大,用到了三种不同的landsat传感器的影像,先分别介绍这三种传感器的波段信息: 1.1Landsat5TM产品说明 Landsat主题成像仪 (TM)是Landsat4 ...
GEE学习记录（一）基于GEE利用LANDSAT 8数据计算遥感生态指数（RSEI）
最近老师让看一下关于GEE的东西,实现大面积的反演.计算地表温度等,也算熟悉一下.参考网上很多大佬的文章,按照自己的思路和想法算出了RSEI,参考的文章都有列出来. 目录所用数据集影像数据矢量数 ...
Python实现遥感生态指数计算
最近在做一些遥感相关的图像处理项目,涉及到遥感生态指数的计算.由于项目要求Python实现,搜索互联网关于Python实现的遥感生态指数计算程序资料很少,于是就自己实现了一个并分享在这里,供需要的朋友 ...
部分常见遥感指数(RSEI)原理与计算方法，以及效果比较
目录前言指数使用前需知遥感指数植被指数归一化植被指数NDVI 归一化差值山地植被指数NDMVI 增强型植被指数EVI 比值植被指数SR/植被指数RVI 差值植被指数DVI 调节土壤的植被指数 ...
《区域生态环境变化的遥感评价指数》笔记
当前卫星遥感对地观测系统以其宏观.快速.实时的优点在生态环境领域得到广泛的应用.利用各种遥感指数来对森林.草地.城市.河流乃至整个流域的生态系统进行监测和评价,已经是生态环境保护领域的重要组成部分. ...
GEE 遥感特征指数谐波分析以Sentinel-2 NDVI为例（含完整代码链接）
谐波分析(Harmonic Regression)是常用的云缺失填补方法, 它是一种频域时序列分析,从原始时间序列数据的正弦和余弦函数重构周期/季节性波.对遥感特征指数使用谐波分析,有利于反映特征的周 ...
Google Earth Engine （GEE）——利用两种方式进行EVI指数（含函数的两种定义方式）
如何快速使用波段进行指数的计算,我们这里利用两种方式进行EVI指数计算,一种是利用expression的方式进行分析,虽然两种方法的结算结果都一样,但是代码有多有少,大家可以参考使用,但是两者的作用对 ...
使用GEE计算多种时序植被指数/建筑指数/水体指数
以Landsat5 2009年-2010年数据为源数据,使用 JavaScript 分别计算区域内逐月的 NDVI.EVI.RVI.DVI.ARVI.SAVI.NDBI.NDWI.MNDWI 九种指数 ...
hopper_如何利用卫星收集的遥感数据轻松对蚱hopper中的站点进行建模
hopper 建筑学与数据科学 (Architectonics and Data Science) Understanding the site and topography are crucial ...

利用GEE计算城市遥感生态指数(RSEI)——Landsat 8为例

文章目录

前言

第一步：定义研究区，自行更换自己的研究区

第二步：加载数据集,定义去云函数

第三步：主函数,计算生态指标

第四步：PCA融合，提取第一主成分

第五步：利用PC1，计算RSEI，并归一化

完整代码

总结

利用GEE计算城市遥感生态指数(RSEI)——Landsat 8为例相关推荐

最新文章

热门文章