【Python】基于Python计算长时间遥感栅格图像的像元值变化度（斜率）和变异系数

1 简介与技术流程

之前看到一篇文章，用变异系数（CV）计算年际初级净生产力（NPP）的变异系数，研究一个地区的产量稳定性。下面是这篇文章（Hong C, Jin X, Ren J, et al. 2019. Satellite data indicates multidimensional variation of agricultural production in land consolidation area. Science of the Environment, 653: 735-747.）关于这部分的截图，供大家参考。上面一张图是他材料与方法中的介绍，下面一张图是他的结果。


总的来说，就是把多幅影像叠加在一起，逐像元构建一组时间序列，然后计算它的变异系数。比方说，你研究区覆盖的像元是100*100的，时间序列是从2011到2015每年一幅影像，那么你就需要构建10000组时间序列，每组里面包含5个数据。我的思想是，把它放在一个坐标系中，横轴是时间，纵轴是每个时间下的像元值，一组时间序列构建一条折线，比较直观的可以参考我论文里的图，如下（我这篇论文还在撰写，所以不能完全放出来，作了一些修改）。

基于这个思想，我们可以计算得出遥感栅格图像的变异系数。然后我再发散一下，发现还可以通过这个计算变化度（斜率）。我结合的实例也是初级净生产力（NPP），用斜率来探究一个地区的生产力增长情况（Production Improvement）；用变异系数来探究生产稳定性（Production Stability）。

2 数据及其预处理

首先，我通过CASA模型计算出4个年份下的初级净生产力（NPP），栅格数据存成geotif格式（后缀是.tif）。

然后把tif文件放在一个找得到的位置，等下要用。

3 代码

3.1 要调用的包

from osgeo import gdal
import numpy as np
from sklearn import linear_model
import copy

osgeo里的gdal是用来处理tif数据的。
numpy是用来做数据分析统计什么的。
sklearn里的linear_model是用来用直线拟合时间序列求斜率的。
copy用里面的深拷贝以防破坏影像数据的原矩阵。

3.2 代码函数：计算斜率和变异系数

首先，我们先写一个普通的斜率和变异系数的函数做个铺垫。这里写的函数先不考虑处理影像，单独就处理一个数组的斜率和变异系数。

# 计算斜率
def calculate_slope(data):
    reg = linear_model.LinearRegression()
    reg.fit(np.array(range(len(data))).reshape(-1, 1), np.array(data).reshape(-1, 1))
    slope = reg.coef_ # 斜率
    intercept = reg.intercept_ # 截距
    slope = slope[0][0]
    intercept = intercept[0]
    return slope


# 计算变异系数
def coefficient_of_variation(data): # 变异系数
    mean = np.mean(data) # 计算平均值
    std = np.std(data, ddof=0) # 计算标准差
    cv = std/mean
    return cv

3.3 代码：计算栅格图像的斜率

然后开始涉及到栅格图像的斜率计算。这里需要调用上面写的函数。
这里输入的参数images是tif数据路径的数组形式，通过for in逐个读取（要保证输入的的多幅影像的空间范围、像元数和坐标系一致）。outpath是计算后输出的斜率图像。

# 栅格图像组计算斜率
def slope(images, outpath):
    images_pixels = [] # 存放多个图像像元矩阵的空数组
    for image in images:
        print(f'正在读取栅格图像： {image}')
        tif = str(image)
        open_tif = gdal.Open(tif) # 打开栅格图像
        band = open_tif.GetRasterBand(1).ReadAsArray() # 获取波段的矩阵
        images_pixels.append(band) # 把图像的矩阵加入到数组中
    SLP = copy.deepcopy(images_pixels[0]) # 获取一个矩阵作为要写入的模板
    # 读取像元，并计算变异系数
    print('正在计算')
    for i in range(len(SLP)):
        for j in range(len(SLP[1])):
            if SLP[i][j] >= 0:
                SLP_data = [] # 存放ij坐标下像元值的数组，以计算变异系数
                for px in range(len(images)): # 遍历多个图像下ij坐标的像元值
                    SLP_data.append(images_pixels[px][i][j]) # 同一坐标的多点加入数组，以计算变异系数
                SLP_value = calculate_slope(SLP_data) # 变异系数计算
                SLP[i][j] = SLP_value # 写入该坐标下的变异系数
    print('计算完成')
    # 保存栅格图像
    gtiff_driver = gdal.GetDriverByName('GTiff')
    out_tif = gtiff_driver.Create(outpath, SLP.shape[1], SLP.shape[0], 1, gdal.GDT_Float32)
    # 将数据坐标投影设置为原始坐标投影
    out_tif.SetProjection(open_tif.GetProjection())
    out_tif.SetGeoTransform(open_tif.GetGeoTransform())
    out_band = out_tif.GetRasterBand(1)
    out_band.WriteArray(SLP)
    out_band.FlushCache()
    print('栅格图像组斜率计算完成')

3.4 代码：计算栅格图像的变异系数

与计算斜率类似。

# 栅格图像组计算变异系数
def CV(images, outpath):
    images_pixels = [] # 存放多个图像像元矩阵的空数组
    for image in images:
        tif = str(image)
        open_tif = gdal.Open(tif) # 打开栅格图像
        band = open_tif.GetRasterBand(1).ReadAsArray() # 获取波段的矩阵
        images_pixels.append(band) # 把图像的矩阵加入到数组中
    CV = copy.deepcopy(images_pixels[-1]) # 获取一个矩阵作为要写入的模板
    # 读取像元，并计算变异系数
    for i in range(len(CV)):
        for j in range(len(CV[1])):
            CV_data = [] # 存放ij坐标下像元值的数组，以计算变异系数
            for px in range(len(images)): # 遍历多个图像下ij坐标的像元值
                CV_data.append(images_pixels[px][i][j]) # 同一坐标的多点加入数组，以计算变异系数
            CV_value = coefficient_of_variation(CV_data) # 变异系数计算
            CV[i][j] = CV_value # 写入该坐标下的变异系数
    # 保存栅格图像
    gtiff_driver = gdal.GetDriverByName('GTiff')
    out_tif = gtiff_driver.Create(outpath, CV.shape[1], CV.shape[0], 1, gdal.GDT_Float32)
    # 将数据坐标投影设置为原始坐标投影
    out_tif.SetProjection(open_tif.GetProjection())
    out_tif.SetGeoTransform(open_tif.GetGeoTransform())
    out_band = out_tif.GetRasterBand(1)
    out_band.WriteArray(CV)
    out_band.FlushCache()
    print('栅格图像组变异系数计算完成')

3.5 代码：函数的调用与结果

我这里用的是2017年到2020年四年的四幅影像。影像长什么样在“数据及其预处理”部分已经展示过了。
下面是执行处理的代码。

images = [
    r'E:\2017.tif',
    r'E:\2018.tif',
    r'E:\2019.tif',
    r'E:\2020.tif',
]

slope(images, r'E:\斜率2017_2020.tif')
CV(images, r'E:\变异系数2017_2020.tif')

这里我把images单独拿出来存成数组作为输入，来保证代码美观（记得一定要按时间顺序排）。
然后后面直接调用上面写的slope函数和CV函数计算斜率和变异系数。images是输入栅格时间序列图像路径，后面那个参数是输出存放的路径。
结果如下面，这里只展示斜率的。

我这里仅结合我自己的实验做个展示，一个时间序列只有4幅影像，如果要做更长的几十幅几百幅的都不是问题，如果太多的话可以自己另写个批量处理什么的。

4 后记

之前还是人文地理学专业的我，这两个月才在学习遥感，可能有些地方不太专业不太科学，希望诸位同行前辈不吝赐教或者可以和我共同探讨一下。可以在csdn私信我或者联系我邮箱（chinshuuichi@qq.com），不过还是希望大家邮箱联系我，我最近不经常上csdn。

如果觉得有帮助，希望大家能打赏支持一下（我的乞讨码在下面）~
-----------------------分割线（以下是乞讨内容）-----------------------