# Python实现遥感生态指数计算

``````"""

"""
import os
import time
import numpy as np
from scipy import misc
from osgeo import gdal, gdalconst
from sklearn.decomposition import PCA

class CalIntegratedEcoIndex:

def __init__(self, img_path, res_save_dir):
self.deno_bias = 0.00001  # 分母偏置，防止除0
self.res_save_dir = res_save_dir

"""读取遥感数据信息"""

img_width = dataset.RasterXSize
img_height = dataset.RasterYSize
img_data = np.array(dataset.ReadAsArray(0, 0, img_width, img_height), dtype=float)  # 将数据写成数组，对应栅格矩阵

del dataset

return img_width, img_height, img_data

def get_wet_degree(self):
"""获取湿度指标"""
return 0.2626 * self.img[0] + 0.2141 * self.img[1] + 0.0926 * self.img[2] + \
0.0656 * self.img[3] - 0.7629 * self.img[4] - 0.5388 * self.img[6]

def get_green_degree(self):
"""获取绿度指标"""
return (self.img[3] - self.img[2]) / (self.img[3] + self.img[2] + self.deno_bias)

def get_temperature(self):
"""获取热度指标"""
gain = 3.20107  # landsat5 第6波段的增益值
bias = 0.25994  # 第6波段的偏置值

K1 = 606.09
K2 = 1282.71
_lambda = 11.45
_rho = 1.438e10-2
_epsilon = 0.96  # 比辐射率

DN = self.img[4] * 299/1000 + self.img[3] * 587/1000 + self.img[2] * 114/1000  # 获取象元灰度值

L6 = gain * DN + bias
T = K2 / np.log(K1 / L6 + 1)
LST = T / (1 + (_lambda * T / _rho) * np.log(_epsilon))

return LST

def get_dryness_degree(self):
"""获取干度指标"""
band5_plus_band3 = self.img[4] + self.img[2]
band4_plus_band1 = self.img[3] + self.img[0]

SI = (band5_plus_band3 - band4_plus_band1) / (band5_plus_band3 + band4_plus_band1 + self.deno_bias)

left_expr = 2 * self.img[4] / (self.img[4] + self.img[3] + self.deno_bias)
right_expr = self.img[3] / (self.img[3] + self.img[2] + self.deno_bias) + \
self.img[1] / (self.img[1] + self.img[4] + self.deno_bias)
IBI = (left_expr - right_expr) / (left_expr + right_expr + self.deno_bias)

NDBSI = (IBI + SI) / 2

return NDBSI

def arr2img(self, save_path, arr):
misc.imsave(save_path, arr)

def normlize(self, img_arr):
arr = np.array(img_arr)
return (arr - arr.min()) / (arr.max() - arr.min())

def get_rsei(self):
"""获取遥感生态指数"""
wet = self.get_wet_degree()
ndvi = self.get_green_degree()
lst = self.get_temperature()
ndbsi = self.get_dryness_degree()

data = np.array([self.normlize(wet), self.normlize(ndvi), self.normlize(lst), self.normlize(ndbsi)])
data = data.reshape(data.shape[0], -1).T

pca = PCA(n_components=1)
rsei = self.normlize(1 - np.reshape(pca.fit_transform(data), newshape=(self.img_height, self.img_width)))

return wet, ndvi, lst, ndbsi, 1 - rsei

def save_result(self, wet, ndvi, lst, ndbsi, rsei):
"""将各指数结果保存为图片"""
res_dict = {"wet": wet, "green": ndvi, "temperature": lst, "dry": ndbsi, "rsei": rsei}

for k, v in res_dict.items():
self.arr2img(os.path.join(self.res_save_dir, k + ".jpg"), v)

def get_average_index(self, wet, ndvi, lst, ndbsi, rsei):
"""获取归一化后各指标的平均值"""
return np.mean(self.normlize(wet)), np.mean(self.normlize(ndvi)), \
np.mean(self.normlize(lst)), np.mean(self.normlize(ndbsi)), np.mean(self.normlize(rsei))
``````

``````if __name__ == '__main__':
start = time.time()
cal = CalIntegratedEcoIndex("/your/path/to/Landsat/image", "./result")
wet, ndvi, lst, ndbsi, rsei = cal.get_rsei()
aver_wet, aver_ndvi, aver_lst, aver_ndbsi, aver_rsei = cal.get_average_index(wet, ndvi, lst, ndbsi, rsei)
cal.save_result(wet, ndvi, lst, ndbsi, rsei)
end = time.time()

print("归一化后各指数的平均值\n湿度", aver_wet, "绿度", aver_ndvi, "热度", aver_lst, "干度", aver_ndbsi,
"遥感生态指数", aver_rsei, "Total cost time %.2f s" % (end - start))
``````