matlab进行theil_sen斜率估计，

将时间序列数据的变化速率赋值到栅格上，最终结果得到一幅变化速率栅格数据

matlab的矩阵运算速度比python快！
注意代码中有一个判断语句，用来确定有效像元，代码中只将大于0的像元参与计算，若需要计算所有像元则删掉下面代码块，或者将阈值改成一个最小值。
if min(data) > 0 %
end

clear
[a, R] = geotiffread('F:\UHI\四种类型城市_5.23\LST_4\shizong3\2001.tif'); % 读取初始年份的遥感图像数据
[m, n] = size(a); 
cd = 2021 - 2001 + 1;  % 计算数据的年份数量
datasum = zeros(m * n, cd) + 0; % 创建用于存储数据的矩阵，并初始化为全零
p = 1; % 初始化计数器
for year = 2001:2021 % 遍历每一年数据
    % 构建当前年份的文件路径
    filename = ['F:\UHI\四种类型城市_5.23\LST_4\shizong3\', int2str(year), '.tif'];
    data = importdata(filename);
    data = reshape(data, m * n, 1); % 将数据从二维矩阵形式转换为一维列向量
    datasum(:, p) = data;
    p = p + 1;
end
result = zeros(m, n); % 创建用于存储结果的矩阵，并初始化为全零
for i = 1:size(datasum, 1) % 遍历每个像素点
    data = datasum(i, :); % 获取当前像素点在不同年份的数据
    if min(data) > 0 % 有效格点判定，我这里有效值在0以上，如果需要将全部像素参与计算则删除该判断代码块，或者改成一个比较小的数值
        valuesum = [];
        % 计算不同年份之间的斜率
        for k1 = 2:cd
            for k2 = 1:(k1 - 1)
                cz = data(k1) - data(k2);
                jl = k1 - k2;
                value = cz / jl;  % 计算斜率
                valuesum = [valuesum; value];
            end
        end 
        value = median(valuesum);% 计算斜率的中位数
        result(i) = value;   % 将斜率赋值给结果矩阵中的对应位置
    end
end
result_matrix = reshape(result, m, n);% 将计算得到的结果矩阵转换为原始图像的二维矩阵形式

% 设置输出结果的文件路径
output_filename = 'F:\UHI\四种类型城市_5.23\LST_4\result\tif\tif_8.9\shizong.tif';

% 将斜率栅格存为 GeoTIFF 格式的文件
geotiffwrite(output_filename, result_matrix, R);

结果

代码改自：

https://blog.csdn.net/weixin_45024495/article/details/121043230

感谢大家花时间来阅读本文 作者水平有限 有失误之处请大家斧正！