栅格图像线性回归趋势分析

#导入必要的库

import os
import numpy as np
import rasterio as ras
from scipy.stats import linregress

• os 用于文件路径操作。
• numpy 用于数值计算。
• rasterio 用于读取和写入遥感影像数据。
• scipy.stats.linregress 用于进行线性回归分析

#写入栅格影像

# 写影像
def writeImage(image_save_path, height1, width1, para_array, bandDes, transform1):
    with ras.open(
            image_save_path,
            'w',
            driver='GTiff',
            height=height1,
            width=width1,
            count=1,
            dtype=para_array.dtype,
            crs='+proj=latlong',
            transform=transform1,
    ) as dst:
        dst.write_band(1, para_array)
        dst.set_band_description(1, bandDes)
    del dst

 #定义线性回归函数

def linear_trend(image_path, output_path):
    filepaths = [os.path.join(image_path, file) for file in os.listdir(image_path) if file.endswith('.tif')]
    num_images = len(filepaths)

 #读取影像数据

    # 读取影像数据
    img1 = ras.open(filepaths[0])
    transform1 = img1.transform
    height1 = img1.height
    width1 = img1.width
    array1 = img1.read(1)
    img1.close()

• 读取第一张影像，获取其投影变换参数、尺寸和数据。

#创建三维数组存储时间序列影像

    # 创建3D数组
    array3d = np.zeros((num_images, height1, width1))
    array3d[0, :, :] = array1

    for i, filepath in enumerate(filepaths[1:], start=1):
        img = ras.open(filepath)
        array3d[i, :, :] = img.read(1)
        img.close()

 

#创建slop增长率数组 

    slope_array = np.full((height1, width1), np.nan)

#遍历每个像元，提取时间序列数据 ，使用 linregress 进行线性回归分析，计算斜率并存储到 slope_array 中。

    for x in range(height1):
        for y in range(width1):
            pixel_series = array3d[:, x, y]
            if np.isnan(pixel_series).any():
                continue

            # 创建时间序列
            time_series = np.arange(num_images)
            # 计算线性回归
            slope, intercept, r_value, p_value, std_err = linregress(time_series, pixel_series)
            slope_array[x, y] = slope

    # 保存结果
    slope_save_path = os.path.join(output_path, "slope3.tif")
    writeImage(slope_save_path, height1, width1, slope_array, 'slope', transform1)

if __name__ == '__main__':
    image_path = r"D:\Master\气温\sanbei"
    output_path = r"D:\Master\蒙古"
    linear_trend(image_path, output_path)

完整代码：

import os
import numpy as np
import rasterio as ras
from scipy.stats import linregress


# 写影像
def writeImage(image_save_path, height1, width1, para_array, bandDes, transform1):
    with ras.open(
            image_save_path,
            'w',
            driver='GTiff',
            height=height1,
            width=width1,
            count=1,
            dtype=para_array.dtype,
            crs='+proj=latlong',
            transform=transform1,
    ) as dst:
        dst.write_band(1, para_array)
        dst.set_band_description(1, bandDes)
    del dst


def linear_trend(image_path, output_path):
    filepaths = [os.path.join(image_path, file) for file in os.listdir(image_path) if file.endswith('.tif')]
    num_images = len(filepaths)

    # 读取影像数据
    img1 = ras.open(filepaths[0])
    transform1 = img1.transform
    height1 = img1.height
    width1 = img1.width
    array1 = img1.read(1)
    img1.close()

    # 创建3D数组
    array3d = np.zeros((num_images, height1, width1))
    array3d[0, :, :] = array1

    for i, filepath in enumerate(filepaths[1:], start=1):
        img = ras.open(filepath)
        array3d[i, :, :] = img.read(1)
        img.close()

    slope_array = np.full((height1, width1), np.nan)

    for x in range(height1):
        for y in range(width1):
            pixel_series = array3d[:, x, y]
            if np.isnan(pixel_series).any():
                continue

            # 创建时间序列
            time_series = np.arange(num_images)
            # 计算线性回归
            slope, intercept, r_value, p_value, std_err = linregress(time_series, pixel_series)
            slope_array[x, y] = slope

    # 保存结果
    slope_save_path = os.path.join(output_path, "slope4.tif")
    writeImage(slope_save_path, height1, width1, slope_array, 'slope', transform1)


if __name__ == '__main__':
    image_path = r"D:\Master\中国年度降雨量栅格数据--平均要素差值--2000-2020\三北地区"
    output_path = r"D:\Master\蒙古"
    linear_trend(image_path, output_path)