利用MATLAB进行VFC的逐像元的线性回归(二) 逐像元的线性回归

逐像元的线性回归

tic
[a,R]=geotiffread('1.tif');%先导入某个图像的投影信息，为后续图像输出做准确
info=geotiffinfo('1.tif');
[m,n]=size(a);
years=10; %表示有多少年份需要做回归
data=zeros(m*n,years);
k=1;
for year=1:10 %起始年份
    file=[int2str(year),'.tif'];%注意自己的名字形式，这里使用的名字是年1.tif，根据这个可修改
    bz=importdata(file);
    bz=reshape(bz,m*n,1);
    data(:,k)=bz;
    k=k+1;
end
    xielv=zeros(m,n);
    p=zeros(m,n);
    %b = zeros(m,n);
    F = zeros(m,n);
    %R_2 = zeros(m,n);
for i=1:length(data)
    bz=data(i,:);
    if isnan(bz)    %注意这是进行判断有效值范围，如果有效范围是-1到1，则改成max(bz)>-1即可
        bz = bz';                   %背景不进行回归
        p(i) = NaN;
        xielv(i) = NaN;
        F(i) = NaN;
    elseif numel(unique(bz)) ==2;   %水体像元不进行回归
         bz = bz';
        p(i) = -3;
        xielv(i) = -3;
        F(i) = -3;
    elseif unique(bz) == 0;         %如果是VFC一直为0
        bz = bz';
        p(i) = -4;
        xielv(i) = -4;
        F(i) =-4;
    else
        bz=bz';
        X=[ones(size(bz)) bz];
        X(:,2)=[1:years]';
        [b,bint,r,rint,stats] = regress(bz,X);
        pz=stats(3);
        p(i)=pz;
        xielv(i)=b(2);
        %R_2 = stats(1);
        %B(i) = b(1);     
        F(i) = stats(2);
    end
end
name = 'k_.tif'
name2='P值.tif';
name4 = 'F的值.tif';
geotiffwrite(name,xielv,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(name2,p,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
geotiffwrite(name4,F,R,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);
%一般来说，只有通过显著性检验的趋势值才是可靠的
F;
toc  

VFC其实是植被覆盖度的另一种指数，是利用像元二分模型计算出来的。是NDVI指数的一种进步。
其实逐像元就是在一个像元的基础之后进行了每个像元的操作，但是由于数据量庞大，会导致时间变得非常的长，根据电脑的配置以及数据量的大小决定，像我就得跑10多分钟。需要耐心等待。其次，由于每个人在处理 遥感影像的时候有不同的习惯，比如我会把水体赋值成-3，但是这些操作的本质都是对矩阵赋值，至于程序怎么写还是需要自己多多琢磨的。
这里建议把VFC改为’1-10.tif’不仅管理方便，更是为了后续的读取更加方便。由于在处理VFC的时候我已经做过一些处理了。
主要是借鉴了这篇文章并且自己进行了一系列的修改。

结果

大概要好长时间才能更新了。。之后要认真准备考研了。
基于Matlab的栅格数据一元线性回归及显著性检验（slope趋势分析）

最后，我将在另外一个账号写博客，有兴趣的可以关注。