MATLAB 建立多元线性回归模型——一以栅格数据气候因子与植被NDVI为例

基于先前博主的代码，我应用gpt修改了下，生成一份更方便理解学习的代码供大家参考！

原博主代码链接：基于matlab的长时间栅格序列逐像元多元回归_逐像元多元线性回归代码-CSDN博客

% @author yinlichang3064@163.com
[aa,R]=geotiffread('I:\裁剪降水\2001pre.tif');%先导入纬度数据
info=geotiffinfo('I:\裁剪降水\2001pre.tif');
[m,n]=size(aa);
temsum=zeros(m*n,20);%20表示时间序列长度
presum=zeros(m*n,20);%20表示时间序列长度
NDVIsum=zeros(m*n,20);%20表示时间序列长度
% 导入每年的数据
for year=2001:2020
     tem=importdata(['I:\裁剪气温\',int2str(year),'tm.tif']) ; %根据自己名称修改，本数据名称为2001-tem.tif
     pre=importdata(['I:\裁剪降水\',int2str(year),'pre.tif']) ; %根据自己名称修改，本数据名称为2001-pre.tif
     NDVI=importdata(['I:\裁剪NDVI\',int2str(year),'yndvi.tif']) ; %根据自己名称修改，本数据名称为'2001-NDVI.tif'
     %注意数据的有效范围
     tem(tem<-1000)=NaN;%温度有效范围
     pre(pre<0)=NaN;%有效范围大于0
     NDVI(NDVI<-1)=NaN; %有效范围是-1到1
     %累加时间序列数据
     temsum(:,year-2001+1)=reshape(tem,m*n,1);%在累加时间序列数据时，我们使用 year - 2001 + 1 来计算列索引，这样可以确保索引从1开始，而不是从0开始。
     presum(:,year-2001+1)=reshape(pre,m*n,1);% 降水数据重塑为列向量
     NDVIsum(:,year-2001+1)=reshape(NDVI,m*n,1);
  
end
%多元回归，ndvi=a+b*pre+c*tem
% 初始化存储回归结果的矩阵
intercept = zeros(m, n) + NaN;
pre_slope = zeros(m, n) + NaN;
tem_slope = zeros(m, n) + NaN;
pz = zeros(m, n) + NaN;
r2 = zeros(m, n) + NaN; % 初始化 R^2 矩阵
for i = 1:m * n
    % 提取特定位置的时间序列数据
    pre = presum(i, :);
    tem = temsum(i, :);
    NDVI = NDVIsum(i, :);
    
    % 检查降水数据是否有效
    if min(pre) >= 0
        % 将列向量转换为行向量
        pre = pre(:);
        tem = tem(:);
        NDVI = NDVI(:);
        % 构建设计矩阵 X
        X = [ones(size(NDVI, 1), 1), tem, pre];
        
        % 进行多元线性回归
        [b, ~, r, rint, stats] = regress(NDVI, X);
        
        % 存储回归结果
        intercept(i) = b(1);
        pre_slope(i) = b(3);
        tem_slope(i) = b(2);
        pz(i) = stats(3);
        r2(i) = stats(1);
    end
end

% 定义输出文件夹路径
outputFolder = 'I:\多元线性回归\';
if ~exist(outputFolder, 'dir')
    mkdir(outputFolder);
end

% 保存回归系数、显著性和截距到 GeoTIFF 文件
filenames = {
    'Intercept.tif';
    'Precipitation_Slope.tif';
    'Temperature_Slope.tif';
    'Equation_Significance.tif';
    'Model_R2.tif'
};
for i = 1:length(filenames)
    filename = fullfile(outputFolder, filenames{i});
    
    % 根据文件名提取相应的变量
    switch filenames{i}
        case 'Intercept.tif'
            variable = intercept;
        case 'Precipitation_Slope.tif'
            variable = pre_slope;
        case 'Temperature_Slope.tif'
            variable = tem_slope;
        case 'Equation_Significance.tif'
            variable = pz;
        case 'Model_R2.tif'
            variable = r2;
    end
    % 保存变量到 GeoTIFF 文件
    geotiffwrite(filename, variable, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
end
% 应用显著性检验并保存显著的结果
pz(pz > 0.05) = NaN;
intercept(pz > 0.05) = NaN;
pre_slope(pz > 0.05) = NaN;
tem_slope(pz > 0.05) = NaN;
r2(pz > 0.05) = NaN;
% 再次保存显著性检验后的回归结果
for i = 1:length(filenames)
    filename = fullfile(outputFolder, filenames{i});
    
    % 根据文件名提取相应的变量
    switch filenames{i}
        case 'Intercept.tif'
            variable = intercept;
        case 'Precipitation_Slope.tif'
            variable = pre_slope;
        case 'Temperature_Slope.tif'
            variable = tem_slope;
        case 'Equation_Significance.tif'
            variable = pz;
        case 'Model_R2.tif'
            variable = r2;
    end
    
    % 保存变量到 GeoTIFF 文件
    geotiffwrite(filename, variable, R, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
end

 

gpt改的代码 顺利运行啦！