对时间序列栅格数据进行趋势分析,最后得到趋势图,可以用来分析栅格数据在每个栅格位置的变化趋势。
数据准备:时间序列栅格数据
ps: 最好按照时间序列规则命名
代码:
tic;
clc;
clear;
s1 = 'sm\2000.tif';% 栅格数据存储在sm文件夹下
[A,R] = geotiffread(s1);
info = geotiffinfo(s1);%获取坐标系信息
x = (1:21); %自变量
[row,col] = size(A); %行列数
M = zeros(21,row,col)+NaN; %建立空数据立方体
%将时间序列数据读入空数据立方体
for year = 2000:2020
s2 = ['sm\',int2str(year),'.tif'];
d = geotiffread(s2);
d(d<-200) = NaN; % 本数据的填充值为负数,所以要去掉填充值,不然会参与运算因此在计算之前要先查看
数据内容,如果有填充值就要根据实际情况去掉
M(year-1999,:,:) = d;% 年份放在第一个位置
end
Slope = zeros(row,col); % 预先建立斜率空数组
Rsq = zeros(row,col); % R2
Pvalue = zeros(row,col)+1; % P值
B = zeros(row,col); % 截距
%一元线性回归
for i = 1:row
for j = 1:col
Y = M(:,i,j);%因变量
X = [ones(length(Y),1),x'];%自变量
[b,bint,r,rint,stats] = regress(Y,X);
Pvalue(i,j) = stats(3); %P值
Slope(i,j) = b(2); %斜率
B(i,j) = b(1); %截距
Rsq(i,j) = stats(1); %R2
end
disp(row);
end
Slope(isnan(Slope)) = 0;
Rsq(isnan(Rsq)) = 0;
Pvalue(isnan(Pvalue)) = 1; %替换空值
%设置存储路径
f1 = 'Trend\Slope.tif';%
f2 = 'Trend\Pvalue.tif';
f3 = 'Trend\Rsq.tif';
f4 = 'Trend\B.tif';
%导出数据
geotiffwrite(f1,Slope,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(f2,Pvalue,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(f3,Rsq,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(f4,B,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
toc;
PS:代码和时间序列数据文件要放在同一个文件夹下,保证路径正确(如下图:sm中存储的是时间序列数据,右边为matlab代码,中间是趋势分析结果存放的位置,在运行程序前确保三个文件都处于同一文件夹下),其次,Trend文件夹需要在运行代码前先建立。