基于GEE使用Landsat 8和Landsat 5影像计算RSEI
一、使用2018年landsat 8影像数据计算遥感生态指数RSEI
//could masking,GEE上有现有的去云函数
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);
}
var dataset1 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.filterDate('2018-x-x', '2018-x-x')
.filterBounds(roi)
.filter(ee.Filter.lte('CLOUD_COVER',5))
.sort('CLOUD_COVER')
.map(maskL8sr)
//.mean()
print(dataset1)
var image1=dataset1.mosaic().clip(roi)
print(image1,'maskL8sr')
//Map.addLayer(image1,{bands:["SR_B4","SR_B3","SR_B2"],gamma:1.3,max:108,min:15},'Image1');
//Map.addLayer(roi,{color:'yellow',fillColor: "00000000", width: 1},'roi');
var jrc = jrcwater18.select('b1')
var image = image1.updateMask(jrc.eq(1))
//Map.addLayer(image,{bands:["SR_B4","SR_B3","SR_B2"],gamma:1.3,max:108,min:15},'Image');
print(image,'water and cloud')
//计算SI
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;
}
//计算IBI
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;
}
//计算湿度指数Wet
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;
}
//外部导入的NDVI重采样
var NDVI=ndvi.reproject('EPSG:4326',null,1000);
//normalization function
function img_normalize(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());
return normalize;
}
//使用MODIS的地表温度数据
var dataset2 = ee.ImageCollection('MODIS/061/MOD11A2')
.filterDate('2018-x-x', '2018-x-x')
.mean()
//var dataset_no_water = dataset.updateMask(NDWI_mask).clip(table);
var dataset2_no_water = dataset2.updateMask(jrc.gte(0.1)).clip(roi);
//将地表温度单位转换为摄氏度
var lst = dataset2_no_water.expression(
'a*0.02-273.15',
{
a:dataset2_no_water.select('LST_Day_1km'),
});
var SI = SI_cal(image);
print(SI,'SI');
var IBI = IBI_cal(image);
var NDBSI =(SI.add(IBI)).divide(2.0);
var wet = Wet_cal(image);
var unit_ndvi = img_normalize(NDVI);
image=image.addBands(unit_ndvi.rename('NDVI').toFloat())
var unit_NDBSI = img_normalize(NDBSI);
image=image.addBands(unit_NDBSI.rename('NDBSI').toFloat())
var unit_Wet = img_normalize(wet);
image=image.addBands(unit_Wet.rename('Wet').toFloat())
var unit_lst = img_normalize(lst);
image=image.addBands(unit_lst.rename('LST').toFloat())
print(image)
var bands = ["Wet","NDVI","NDBSI","LST"]
var sentImage =image.select(bands)
var region = roi;
var image2= sentImage.select(bands);
var scale = 1000;
var bandNames = image2.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 = image2.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: region,
scale: scale,
maxPixels: 1e9
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image2.subtract(means);
//主成分分析函数,GEE上有
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')
});
print(rsei_un_unit,'rsei_un_unit')
var rsei = img_normalize(rsei_un_unit);
//Map.addLayer(rsei, {}, 'PCA');
//Map.addLayer(lst, visParams1, "LST");
//Map.addLayer(NDBSI, {}, "NDBSI");
//Map.addLayer(ndvi, visParams, "NDVI");
//Map.addLayer(wet, visParams2, "Wet");
//print(wet)
print(rsei);
Map.addLayer(rsei, {}, "rsei");
//导出结果
Export.image.toDrive({
image: rsei,//分类结果
description: 'rsei',//文件名
folder: 'test',
scale: 30,//分辨率
region: roi,//区域
maxPixels:34e10//此处值设置大一些,防止溢出
});
二、使用2009年landsat 5影像数据计算遥感生态指数RSEI
//could masking
function maskL5sr(image) {
// Bit 0 - Fill
// Bit 1 - Dilated Cloud
// Bit 2 - Unused
// 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 thermalBand = image.select('ST_B6').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(thermalBand, null, true)
.updateMask(qaMask)
.updateMask(saturationMask);
}
var dataset1 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
.filterDate('2009-x-x', '2009-x-x')
.filterBounds(roi)
.filter(ee.Filter.lte('CLOUD_COVER',5))
.sort('CLOUD_COVER')
.map(maskL5sr)
//.mean()
print(dataset1)
var image1=dataset1.mosaic().clip(roi)
print(image1,'maskL5sr')
//Map.addLayer(image1,{bands:["SR_B3","SR_B2","SR_B1"],gamma:1.3,max:108,min:15},'image1');
var jrc = jrcwater09.select('b1')
var image = image1.updateMask(jrc.eq(1))
Map.addLayer(image,{bands:["SR_B3","SR_B2","SR_B1"],gamma:1.3,max:108,min:15},'Imagewc');
print(image,'water and cloud')
function SI_cal(img) {
var blue = img.select("SR_B1");
var red = img.select("SR_B3");
var nir = img.select("SR_B4");
var swir1 = img.select("SR_B5");
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_B2");
var red = img.select("SR_B3");
var nir = img.select("SR_B4");
var swir1 = img.select("SR_B5");
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 Wet_cal(img) {
var blue = img.select("SR_B1");
var green = img.select("SR_B2");
var red = img.select("SR_B3");
var nir = img.select("SR_B4");
var swir1 = img.select("SR_B5");
var swir2 = img.select("SR_B7");
var wet_temp =blue.multiply(0.0315)
.add(green.multiply(0.2021))
.add(red.multiply(0.3102))
.add(nir.multiply(0.1594))
.add(swir1.multiply(-0.6806))
.add(swir2.multiply(-0.6109))
//print(wet_temp);
return wet_temp;
}
var NDVI=ndvi2009.reproject('EPSG:4326',null,1000);
//var img_mean=lst.mean().reproject('EPSG:4326',null,1000);
//normalization function
function img_normalize(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());
return normalize;
}
var dataset2 = ee.ImageCollection('MODIS/061/MOD11A2')
.filterDate('2009-X-X', '2009-X-X')
.mean()
//var dataset_no_water = dataset.updateMask(NDWI_mask).clip(table);
var dataset2_no_water = dataset2.updateMask(jrc.gte(0.1)).clip(roi);
var lst = dataset2_no_water.expression(
'a*0.02-273.15',
{
a:dataset2_no_water.select('LST_Day_1km'),
});
var SI = SI_cal(image);
print(SI,'SI');
var IBI = IBI_cal(image);
var NDBSI =(SI.add(IBI)).divide(2.0);
var wet = Wet_cal(image);
var unit_ndvi = img_normalize(NDVI);
image=image.addBands(unit_ndvi.rename('NDVI').toFloat())
var unit_NDBSI = img_normalize(NDBSI);
image=image.addBands(unit_NDBSI.rename('NDBSI').toFloat())
var unit_Wet = img_normalize(wet);
image=image.addBands(unit_Wet.rename('Wet').toFloat())
var unit_lst = img_normalize(lst);
image=image.addBands(unit_lst.rename('LST').toFloat())
print(image)
var bands = ["Wet","NDVI","NDBSI","LST"]
var sentImage =image.select(bands)
var region = roi;
var image2= sentImage.select(bands);
var scale = 1000;
var bandNames = image2.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 = image2.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: region,
scale: scale,
maxPixels: 1e9
});
var means = ee.Image.constant(meanDict.values(bandNames));
var centered = image2.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*0.605-pc2*0.395' ,
{
constant: 1,
pc1:pcImage.select('pc1'),
});
print(rsei_un_unit,'rsei_un_unit')
var rsei = img_normalize(rsei_un_unit);
//Map.addLayer(rsei, {}, 'PCA');
//Map.addLayer(lst, visParams1, "LST");
//Map.addLayer(NDBSI, {}, "NDBSI");
//Map.addLayer(ndvi, visParams, "NDVI");
//Map.addLayer(wet, visParams2, "Wet");
//print(wet)
print(rsei);
Map.addLayer(rsei, {}, "rsei");
Export.image.toDrive({
image: rsei,//分类结果
description: 'rsei',//文件名
folder: 'test',
scale: 30,//分辨率
region: roi,//区域
maxPixels:34e10//此处值设置大一些,防止溢出
});
基于GEE使用Landsat 8和Landsat 5影像计算RSEI相关推荐
- GEE(8):使用MODIS填补由去云后的Landsat影像计算得到的NDVI数据
最近想要在GEE中使用Landsat影像计算一下广州的NDVI值,发现这片区域云覆盖较多,去云以后部分月份的数据很少,就造成NDVI计算结果缺失的问题.经过查阅相关资料,可以使用MODIS的NDVI产 ...
- gee去云处理Landsat、Sentinel和Modis影像
gee去云处理Landsat.Sentinel和Modis影像 1.Landsat 和 MODIS影像的去云函数 function cloudfree_mod09a1(image){var qa = ...
- 【基于GEE计算年度植被覆盖度】
基于GEE计算年度植被覆盖度 1.植被覆盖度计算公式参考地理学报的<2001-2013年华北地区植被覆盖度与干旱条件的相关分析> 2.代码 var roi=ee.FeatureCollec ...
- GEE:基于GEE的单个湖泊的实时水体提取(以武汉东湖为例)
前言 博主不主修遥感方向,不是专业人士,由于毕设需要使用GEE,故临时学习并做了记录,由于是博主自己钻研的,不知道有无其他更便捷的方式.若文中有错误欢迎指正. 在做毕设的时候是对武汉市所有湖泊进行分别 ...
- 【Earth Engine】基于GEE对季节性地物进行分类(多源数据叠图+监督分类)
目录 1 简介与摘要 2 采样点的选取 3 合成多季节多波段影像 4 分类器参数设置 5 监督分类的实现与分类精度计算 6 影像的显示与结果 7 后记 1 简介与摘要 最近导师给了个新方向,其中要用到 ...
- GEEer成长日记十九:使用Landsat 8影像计算水体指数NDWI、MNDWI并下载到本地
目录 一.NDWI和MNDWI计算公式介绍 1.NDWI(归一化差异水体指数) 2.MNDWI(改进的归一化差异水体指数) 二.使用Landsat8影像计算NDWI和MNDWI 1.获取Landsat ...
- 基于GEE与哨兵1号影像数据提取水体
在进行光学数据提取水体时,常会发现部分时间.部分区域内由于云的存在而出现大面积的空白区域,从而使得我们提取水体面积过程中存在不精确的问题.因此,基于微波数据进行云的提取就成为我们较好的一个选择,通过阅 ...
- 如何使用 Landsat 8 卫星影像计算地表温度
首先当然是要下载Landsat 8卫星影像.打开以下网址 https://eos.com/landviewer/或者https://earthexplorer.usgs.gov/或者 访问之前数据资源 ...
- pandas使用cut函数基于分位数进行连续值分箱(手动计算分位数)处理后出现NaN值原因及解决
pandas使用cut函数基于分位数进行连续值分箱(手动计算分位数)处理后出现NaN值原因及解决 目录 pandas使用cut函数基于分位数进行连续值分箱(手动计算分位数)处理后出现NaN值原因及解决 ...
- python 语义similarity_GitHub - samelltiger/word_similarity: 基于《知网》的语义相似度计算 python2.7 API...
基于<知网>的语义相似度计算 python2.7 API 本项目使用python语言实现根据义原树来计算词语之间的语义相似度,并提供对应的 API. 词语距离有两类常见的计算方法,一种是根 ...
最新文章
- Oracle可以处理LOB字段的常用字符函数
- [学习windows/记录篇]安装TMG防火墙(三向外围)
- HDU5863 cjj's string game(DP + 矩阵快速幂)
- 数据结构实验三:Huffman树及Huffman编码的算法实现
- 浅析SparkRPC源码(spark2.11)
- M1 mac 使用docker 安装mysql
- 【BZOJ-1952】城市规划 [坑题] 仙人掌DP + 最大点权独立集(改)
- 做购物车系统时利用到得几个存储过程
- Python读取IRIS数据集并转换为PaddlePaddle中使用的reader
- linux查看文件时显示行号,linux中查看文件时显示行号
- VS Code离线安装C/C++插件cpptools-linux-aarch64.vsix
- fstab损坏无法开机的修复
- 35. Consider alternatives to virtual functions
- Linux内核ncsi驱动源码分析(二)
- 工程力学(1)-公理以及简单的受力分析
- 齐天大圣蟠桃园吃桃子
- 报错 Error from server (InternalError): an error on the server (““) has prevented the request from suc
- 虚拟机VMware安装Ubuntu记录
- 树莓派3B+ 编译Qt源码
- 计算机课图画的变形教案,《义务教育课程标准实验教科书美术(四年级上册)》提示及教学要点...