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
感谢大家花时间来阅读本文 作者水平有限 有失误之处请大家斧正!