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

10 篇文章 22 订阅

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。

如果觉得有帮助,希望大家能打赏支持一下(我的乞讨码在下面)~
-----------------------分割线(以下是乞讨内容)-----------------------
在这里插入图片描述

评论 37
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值