Python Implementation for Remote Sensing Ecological Index Calculation

The author needed a Python solution for remote‑sensing ecological index calculations and, finding few existing resources, implemented a complete script that computes wetness, greenness, temperature, dryness, and a combined ecological index from Landsat multi‑band images.

Below is the full Python code. It reads a Landsat image with GDAL, extracts the required bands, computes each index, normalizes the results, applies PCA to derive the final ecological index, and saves each index as an image.

<code>"""
利用湿度、绿度、热度、干度计算遥感生态指数
"""
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))
</code>

To use the script, provide the path to a Landsat image containing seven bands and a directory for results, then call the class methods as shown in the example below.

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

The script assumes the input image is a Landsat multi‑band raster (seven bands). The original article can be found at the provided CSDN link.