Python+GEE遥感开发之计算遥感生态指数RSEI
前言：遥感生态指数RSEI， 是一种使用卫星遥感影像数据通过反演计算得到的数据，可以用来对城市的生态状况进行快速监测和评价，该指数主要利用主成分分析方法，对植被指数、湿度、地表温度以及建筑指数四项指数指标进行集成。这个指数是由徐涵秋发明的，更多详细信息还请参考阅读：徐涵秋.城市遥感生态指数的创建及其应用[J].生态学报,2013,33(24):7853-7862.
0 使用的遥感数据
本篇博客主要介绍MODIS下计算RSEI的方法和Landsat计算RSEI的方式，两者的差别会在下面的指数计算中给出不同点。
其中MODIS计算RSEI主要用到了MOD09A1波段（1~7波段分别表示为：红光、近红外、蓝光、绿光、近红外2、短波红外1、短波红外2）的数据和MOD11A1下的LST数据。
Landsat计算RSEI主要用到了Landsat8 OLI波段数据。
1 植被指数的计算（NDVI）
植被指数主要推荐使用NDVI，NDVI的计算主要是用到了红外和近红外波段，可以使用MODIS或者Landsat的NDVI产品，也可以使用波段完成计算，这里给出一个使用GEE平台和MODIS下的MOD09A1产品计算年合成的NDVI，具体代码如下所示。当然也可以参考我的其它GEE计算NDVI的博客。
var geometry = ee.FeatureCollection('users/www1573979951/hkh_boundary');
Map.centerObject(geometry,6);
var mod09a1 = ee.ImageCollection('MODIS/061/MOD09A1').filterBounds(geometry);
print(mod09a1)
mod09a1 = mod09a1.map(function(img)
{ var ndvi = img.expression('((nir-red)/(nir+red))',
{
'red': img.select(['sur_refl_b01']).multiply(0.0001),
'nir': img.select(['sur_refl_b02']).multiply(0.0001)
})
// Clip NDBSI values to the range [-1, 1] and mask out values outside this range
var clipped_ndvi = ndvi.clip(geometry).updateMask(ndvi.gte(-1).and(ndvi.lte(1)));
return img = img.addBands(clipped_ndvi.rename('NDVI'));
});
for(var i=2001;i<=2023;i++){
var data_collection = mod09a1.filterDate(i+'-01-01',i+'-12-31').select('NDVI');
var YR_collection = data_collection.mean().clip(geometry);
Export.image.toDrive({
image: YR_collection,
description: i,
fileNamePrefix: i,
scale: 1000,
crs: "EPSG:4326",
region: geometry,
maxPixels: 1e13,
folder: 'Modis'
})
}
2 湿度指数的计算（WET）
2.1 MOD09A1计算WET（GEE代码）
湿度的计算主要用到了MOD09A1产品的7个波段，主要参考CSDN博主“空中旋转篮球”给出的WET的计算方式。
MOD09A1计算WET的GEE代码如下：
var geometry = ee.FeatureCollection('users/www1573979951/hkh_boundary');
Map.centerObject(geometry,6);
var mod09a1 = ee.ImageCollection('MODIS/061/MOD09A1').filterBounds(geometry);
print(mod09a1)
mod09a1 = mod09a1.map(function(img)
{ // 湿度函数:Wet
var Wet = img.expression('((0.1147 * red) + (0.2489 * infrared) + (0.2408 * blue) + (0.3132 * green) - (0.3122 * short_infrared) - (0.6416 * middle_infrared1) - (0.5087 * middle_infrared2))',{
'red': img.select(['sur_refl_b01']).multiply(0.0001),
'infrared': img.select(['sur_refl_b02']).multiply(0.0001),
'blue': img.select(['sur_refl_b03']).multiply(0.0001),
'green': img.select(['sur_refl_b04']).multiply(0.0001),
'short_infrared': img.select(['sur_refl_b05']).multiply(0.0001),
'middle_infrared1': img.select(['sur_refl_b06']).multiply(0.0001),
'middle_infrared2': img.select(['sur_refl_b07']).multiply(0.0001)
})
var clipped_wet = Wet.clip(geometry).updateMask(Wet.gte(-1).and(Wet.lte(1)));
return img.addBands(clipped_wet.rename('WET'))
});
for(var i=2001;i<=2023;i++){
var data_collection = mod09a1.filterDate(i+'-01-01',i+'-12-31').select('WET');
var YR_collection = data_collection.mean().clip(geometry);
Export.image.toDrive({
image: YR_collection,
description: 'wet'+i,
fileNamePrefix: 'wet'+i,
scale: 1000,
crs: "EPSG:4326",
region: geometry,
maxPixels: 1e13,
folder: 'Modis'
})
}
2.2 Landsat8计算Wet（Python代码）
计算方法来源于《基于Python计算Landsat8OLI遥感生态指数RSEI》
def wet(blue, green, red, nir, swir1, swir2):
"""
计算湿度 wet
:param blue: landsat8的B2波段,蓝色波段blue
:param green: landsat8的B3波段,绿色波段green
:param red: landsat8的B4波段,红外波段red
:param nir: landsat8的B5波段,近红外波段nir
:param swir1: landsat8的B6波段,短波红外1 swir1
:param swir2: landsat8的B7波段,短波红外2 swir2
:return: wet 计算结果
"""
blue = blue * 0.0001
green = green * 0.0001
red = red * 0.0001
nir = nir * 0.0001
swir1 = swir1 * 0.0001
swir2 = swir2 * 0.0001
wet = blue * 0.1511 + green * 0.1972 + red * 0.3283 + nir * 0.3407 + swir1 * (-0.7117) + swir2 * (-0.4559)
return wet
3 温度指数的计算（LST）
3.1 MOD11A1获取LST（GEE代码）
这里使用的是MOD11A1产品下的LST产品，取了Night和Day的平均值。
var geometry = ee.FeatureCollection('users/www1573979951/hkh_boundary');
Map.centerObject(geometry,6);
var lst = ee.ImageCollection('MODIS/061/MOD11A1').filterBounds(geometry);
print(lst)
lst = lst.map(function(img)
{ var lst_mean = img.expression('((Day+Night)/2)',
{
'Day': img.select(['LST_Day_1km']),
'Night': img.select(['LST_Night_1km'])
})
return img = img.addBands(lst_mean.rename('LST'));
});
for(var i=2001;i<=2023;i++){
var data_collection = lst.filterDate(i+'-01-01',i+'-12-31').select('LST');
data_collection = data_collection.map(function(img){
var date = img.get('system:time_start');
return img.multiply(0.02).subtract(273.15).set('system:time_start', date);
});
var YR_collection = data_collection.mean().clip(geometry);
Export.image.toDrive({
image: YR_collection,
description: 'lst'+i,
fileNamePrefix: 'lst'+i,
scale: 1000,
crs: "EPSG:4326",
region: geometry,
maxPixels: 1e13,
folder: 'Modis'
})
}
3.2 Landsat获取LST（Python代码）
计算方法来源于《基于Python计算Landsat8OLI遥感生态指数RSEI》
landsat8的第10波段是热红外波段，可以用其辐射亮度结合植被覆盖度和比辐射率来计算热度LST。
def lst(tirs1, ndvi):
"""
计算热度指数 lst
:param tirs1: landsat8的B10波段,热红外波段1 tirs1
:param ndvi: ndvi
:return:
"""
# 选择第Band10辐射亮度波段
tirs1_fsld = tirs1 * 0.1
# 计算植被覆盖度pv，假设最低值为0.05，最大值为0.95
pv = np.power((ndvi-0.05) / (0.95-0.05), 2)
# 计算比辐射率c
c = 0.004 * pv + 0.986
# 计算lst
lst = (tirs1_fsld / (1 + (0.00109 * (tirs1_fsld / 1.438)) * np.log(c))) - 273.15
return lst
4 干度指数的计算（NDBISI）
NDBSI的计算需要先算出裸土指数SI和建筑指数IBI，进而再求出NDBSI。
4.1 MOD09A1计算NDBISI（GEE代码）
NDBISI的计算方式同样来源于CSDN博主“空中旋转篮球”给出的方法。
var geometry = ee.FeatureCollection('users/www1573979951/hkh_boundary');
Map.centerObject(geometry,6);
var mod09a1 = ee.ImageCollection('MODIS/061/MOD09A1').filterBounds(geometry);
mod09a1 = mod09a1.map(function(img) {
// 干度函数:ndbsi = ( ibi + si ) / 2
var ibi = img.expression('(((float(2*b5))/(float(b5+b4))-((float(b4))/(float(b4+b3))+(float(b2))/(float(b2+b5))))/((float(2*b5))/(float(b5+b4))+((float(b4))/(float(b4+b3))+(float(b2))/(float(b2+b5)))))', {
'b2': img.select('sur_refl_b02').multiply(0.0001),
'b3': img.select('sur_refl_b03').multiply(0.0001),
'b4': img.select('sur_refl_b04').multiply(0.0001),
'b5': img.select('sur_refl_b05').multiply(0.0001)
});
var si = img.expression('((float((b3+b5)-(b1+b4)))/((b3+b5)+(b1+b4)))', {
'b1': img.select('sur_refl_b01').multiply(0.0001),
'b2': img.select('sur_refl_b02').multiply(0.0001),
'b3': img.select('sur_refl_b03').multiply(0.0001),
'b4': img.select('sur_refl_b04').multiply(0.0001),
'b5': img.select('sur_refl_b05').multiply(0.0001)
});
var ndbsi = (ibi.add(si)).divide(2);
// Clip NDBSI values to the range [-1, 1] and mask out values outside this range
var clipped_ndbsi = ndbsi.clip(geometry).updateMask(ndbsi.gte(-1).and(ndbsi.lte(1)));
return img.addBands(clipped_ndbsi.rename('NDBISI'));
});
for(var i=2023; i<=2023; i++){
var data_collection = mod09a1.filterDate(i+'-01-01', i+'-12-31').select('NDBISI');
var YR_collection = data_collection.mean().clip(geometry);
Export.image.toDrive({
image: YR_collection,
description: 'ndbisi'+i,
fileNamePrefix: 'ndbisi'+i,
scale: 1000,
crs: "EPSG:4326",
region: geometry,
maxPixels: 1e13,
folder: 'Modis'
})
}
4.2 Landsat计算NDBISI（python代码）
计算方法来源于《基于Python计算Landsat8OLI遥感生态指数RSEI》
def ndbsi(blue, green, red, nir, swir1):
"""
计算干度 ndbsi
:param blue: landsat8的B2波段,蓝色波段blue
:param green: landsat8的B3波段,绿色波段green
:param red: landsat8的B4波段,红外波段red
:param nir: landsat8的B5波段,近红外波段nir
:param swir1: landsat8的B6波段,短波红外1 swir1
:return: ndbsi 计算结果
"""
blue = blue * 0.0001
green = green * 0.0001
red = red * 0.0001
nir = nir * 0.0001
swir1 = swir1 * 0.0001
# 计算裸土指数 si
si = ((swir1 + red) - (nir + blue)) / ((swir1 + red) + (nir + blue))
# 计算建筑指数 ibi
ibi = (2 * swir1 / (swir1 + nir) - (nir / (nir + red) + green / (green + swir1))) / \
(2 * swir1 / (swir1 + nir) + (nir / (nir + red) + green / (green + swir1)))
# 计算干度
ndbsi = (si + ibi) / 2
return ndbsi
5 RSEI的计算（Python代码）
首先对NDVI、LST、WET、NDBISI数据的归一化，然后使用PCA进行主成分分析计算出RSEI指数。
import os
import numpy as np
from osgeo import gdal, gdalnumeric
from sklearn.decomposition import PCA
def normalize(array):
"""
线性归一化
:param array: 数组
:return: 处理后的数组
"""
# 线性归一化
max = np.nanmax(array)
min = np.nanmin(array)
normalize_value = (array-min) / (max-min)
return normalize_value
def rsei(ndvi, wet, ndbsi, lst):
"""
PCA主成分分析并计算RSEI
:param ndvi: ndvi
:param wet: wet
:param ndbsi: ndbsi
:param lst: lst
:return: rsei, pc_ratio(保留成分的方差百分比)
"""
# 将四个指标合成一个多维数组，此时data的shape是(4, m, n)，是一个三维数据
data = np.array([wet, ndvi, lst, ndbsi])
# 对数组进行掩膜，不让nan参与PCA
panel_data = np.hstack([data[i, :, :].flatten().reshape(-1, 1) for i in range(data.shape[0])])
mask = np.isnan(panel_data).sum(axis=1) > 0
train_x = panel_data[~mask, :]
# 主成分分析
pca = PCA(n_components=1)
pca_data = pca.fit_transform(train_x)
# 保留成分的方差百分比(如果用不到可以不用返回，可以去掉)
pc_ratio = pca.explained_variance_ratio_
final_result = np.ones_like(mask) * 1.0
final_result[~mask] = pca_data.flatten()
final_result[mask] = np.NAN
result_data = final_result.reshape(data.shape[1:])
# 计算rsei
rsei = 1 - result_data
# 对rsei进行归一化
rsei = normalize(rsei)
return rsei, pc_ratio
def save_tif(data, file, output):
ds = gdal.Open(file)
shape = data.shape
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(output, shape[1], shape[0], 1, gdal.GDT_Float32) # 以float类型进行存储
dataset.SetGeoTransform(ds.GetGeoTransform())
dataset.SetProjection(ds.GetProjection())
dataset.GetRasterBand(1).WriteArray(data)
def read_tif02(filepath):
data = gdalnumeric.LoadFile(filepath)
data = data.astype(np.float32)
a = data[0][0]
data[data == a] = np.nan
return data
def get_data_list(file_path, out = ""):
list1 = [] # 文件的完整路径
if os.path.isdir(file_path):
fileList = os.listdir(file_path)
if out != "":
for f in fileList:
out_data = out + "\\" + f
list1.append(out_data)
else:
for f in fileList:
pre_data = file_path + '\\' + f # 文件的完整路径
list1.append(pre_data)
return list1
if __name__ == '__main__':
file_ndvi_list = get_data_list(r"E:\AAWORK\work\研究方向\rsei\data\月\ndvi01")
file_lst_list = get_data_list(r"E:\AAWORK\work\研究方向\rsei\data\月\lst01")
file_ndbisi_list = get_data_list(r"E:\AAWORK\work\研究方向\rsei\data\月\ndbisi01")
file_wet_list = get_data_list(r"E:\AAWORK\work\研究方向\rsei\data\月\wet01")
file_rsei_list = get_data_list(r"E:\AAWORK\work\研究方向\rsei\data\月\ndvi01",r"E:\AAWORK\work\研究方向\rsei\data\月\rsei01")
for file_ndvi,file_lst,file_ndbisi,file_wet,file_rsei in zip(file_ndvi_list,file_lst_list,file_ndbisi_list,file_wet_list,file_rsei_list):
ndvi_data = normalize(read_tif02(file_ndvi))
lst_data = normalize(read_tif02(file_lst))
ndbisi_data = normalize(read_tif02(file_ndbisi))
wet_data = normalize(read_tif02(file_wet))
print(file_ndvi)
print(file_lst)
print(file_ndbisi)
print(file_wet)
print(file_rsei)
print("---------")
rsei_data, pc_ratio = rsei(ndvi_data,wet_data,ndbisi_data,lst_data)
save_tif(rsei_data,file_ndvi,file_rsei)
3178
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
