代码分享
方法介绍:Sen 斜率估计用于计算趋势值,通常与MK非参数检验结合使用。即首先计算Sen趋势值,然后使用MK方法判断趋势显著性
示例:1984-NDVI年最大值趋势分析
注意:在对NDVI进行趋势分析时,绝对值0.1以下的NDVI值需要去除
代码1:MKTrend(代码2会用)
function MKResult = MKTrend(X,Alpha)%tic% 时间序列数据的Mann-Kendall趋势分析%% 原假设:H0 Beta == 0% 当 |Z| <= pNorm ,表明在当前显著性水平Alpha下接受原假设(无显著变化)% 当 Z > pNorm , 表明在当前显著性水平Alpha下认为有上升趋势:此时Beta > 0% 当 Z < -pNorm ,表明在当前显著性水平Alpha下认为有下降趋势:此时Beta < 0%% MKResult = [Z,pNorm]; (统计量S的标准化)Z统计量,上侧分位数%% Input% X - 时间序列数据,向量% Alpha- 显著性水平%% Output% Beta- 衡量趋势大小的倾斜度或斜率(整体变化趋势的速率)% MKResult - MK趋势检验的Z统计量及Z分位数,[Z,pNorm]% S_stdS - 统计量S及标准差,[S,sqrt(varS)];% ConInterval - 置信区间% -----------------------------------------------------if (nargin < 2) || (isempty(Alpha)), Alpha = 0.05; endN = length(X);if N < 8, error('Error:MannKendallTrend: 要求数据样本数目不小于8!'); end% if N < 41, error('Error:MannKendallTrend: 要求数据样本数目不小于40!'); endS = 0; % 时间序列X的MK趋势检验的统计量Sfor ii = 1:(N - 1)for jj = ii:NS = S + sign(X(jj) - X(ii));endend% 计算方差varS = ( N * ( N - 1 ) * ( 2 * N + 5 ) ) / 18;% 计算 Z 统计量,该统计量服从标准正态分布if S == 0Z = 0;elseif S > 0Z = (S - 1) / sqrt(varS);elseZ = (S + 1) / sqrt(varS);end% 原假设:H0 Beta == 0% 当 |Z| <= pNorm ,表明在当前显著性水平Alpha下接受原假设% 当 Z > pNorm , 表明在当前显著性水平Alpha下认为有上升趋势% 当 Z < -pNorm ,表明在当前显著性水平Alpha下认为有下降趋势pNorm = norminv(1 - Alpha / 2);MKResult = [Z,pNorm]; % MK趋势检验结果%toc% end of function MannKendallTrend%%end
代码2:调用代码1。MKTrendAnalysis(需要修改的地方已在代码中注释)
clc; clear;
imgdir1='D:\homework2\data\data\NDVIyear\';%%修改为所要处理的数据路径
addpath(genpath(imgdir1));
%% MK趋势分析
filenames = dir([imgdir1 '*.tif']);
for i = 1:numel(filenames)
data(:,:,i) = single(imread(filenames(i).name)); %% 原始数据
end
%%
[row,col, N]=size(data);
timeslice = N;
beg=1984;%%数据起始年份
last=;%%数据结束年份
NA = data(1,1,1);
%MK_para=zeros(row,col,2);
K=zeros(row,col)*NaN;
Z=zeros(row,col)*NaN;
X=zeros(1,timeslice)*NaN;
t=beg:1:last;
Alpha=0.05;%%置信区间
for i=1:row
i
for j=1:col
if ismember(data(1,1,1),data(i,j,:)) % 当某位置的时间序列里有无效的数据时, assign NaN to Z and K
Z(i,j)=-9999;
K(i,j)=-9999;
else
MKResult=MKTrend(data(i,j,:),Alpha);
X=squeeze(data(i,j,:));
p=polyfit(t',X,1);
K(i,j)=p(1); %% 变化量
Z(i,j)=MKResult(1);%%显著性
end
end
end
%%
[~,R0]=readgeoraster('D:\homework2\data\data\NDVIyear\ndvi_1984.tif');%%输入一幅标准的栅格数据来获取属性信息
info = geotiffinfo('D:\homework2\data\data\NDVI year\ndvi_1984.tif'); %%输入一幅标准的栅格数据来获取属性信息
geoTags = info.GeoTIFFTags.GeoKeyDirectoryTag;
outPath = 'D:\homework2\data\data\NDVI year\'; %%输出路径
outName1=[outPath,'Z.tif'];%%输出数据名称
geotiffwrite(outName1,Z,R0,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);%%若输出Slope值,将本行中的Z改为K即可
ARCGIS制图
首先将更改栅格图像Nodata值
在arcgis目录处找到需要处理的栅格数据 - 选中数据 - 单击右键 - 打开属性 - 点击编辑 - 将Nodata值改为-9999 - 计算 - 确定 - 应用
然后对Z值和Slope值分别进行重分类,重分类的结果用栅格计算器相乘,得到最终趋势变化图
其中分类标准如下:
对于Slope:-0.0005 - 0.0005为稳定区域,大于或等于0.0005为植被改善区域,小于-0.0005为植被退化区域
对于Z值:置信水平为0.05,-1.96 - 1.96之间为不显著,绝对值大于1.96为显著。
最终可以分为五类:严重退化、轻微退化、稳定不变、轻微改善、明显改善
结果图如下所示:
END
编辑|南波婉帆
审核| 南波婉琳