Python实现遥感生态指数计算

最近在做一些遥感相关的图像处理项目，涉及到遥感生态指数的计算。由于项目要求Python实现，搜索互联网关于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.img_width, self.img_height, self.img = self.read_img(img_path)
        self.deno_bias = 0.00001  # 分母偏置，防止除0
        self.res_save_dir = res_save_dir

    def read_img(self, img_path):
        """读取遥感数据信息"""
        dataset = gdal.Open(img_path, gdalconst.GA_ReadOnly)

        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))

程序可以按照下面的方式进行使用（注意输入图像是Landsat七波段的图像）：

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))

因为对遥感领域不很了解，这也是我的初次尝试。如果程序中有什么错误及不足，欢迎读者朋友们在评论区批评指正。

转载自：https://blog.csdn.net/hfutdog/article/details/87878787