matlab EOF程序

% eof第一模态图用变量eof的第一列，以此类推；相应的时间系数用pc第一行
clear
clc

x = [2 6 1 5 2;9 4 0 5 4;12 2 55 9 10;4 55 78 2 13];        %原始数据，列为站点或格点，行为时间序列上的值
x(1,:) = x(1,:) - mean(x(1,:));
x(2,:) = x(2,:) - mean(x(2,:));       %距平处理

c = x*x'/5;    %协方差矩阵

[eof,eigenvalue] = eig(c);% eof为特征向量，eigenvalue为特征值
pc = eof' * x;% pc为主成分，即时间系数

eigenvalue = fliplr(flipud(eigenvalue));% flipud矩阵上下翻转；fliplr矩阵左右翻转
lambda = diag(eigenvalue);%lambda(λ)
eof = fliplr(eof);
pc = flipud(pc);

% check
I2 = eof * eof';    % = I
lambda2 = pc * pc'/5;     % = lambda
x2 = eof * pc;     % = x

%计算方差贡献
r = lambda(1,1)/sum(lambda(:,1));

clear;
clc;
in_data = importdata('Precipitation_mo(20-20).txt');
stations = xlsread('数据库对应站名.xlsx');
for i = 2 : length(in_data)
    for j = 2 : 88
        if in_data(i,j) == 9999
            in_data(i,j) = 0;
        end
    end
end

%补上经纬度
n = 1;
for i = 2 : 88
    for j = 1 : 91
        if in_data(1,i) == stations(j,1)
            jwd(n,1) = in_data(1,i);
            jwd(n,2) = stations(j,2);
            jwd(n,3) = stations(j,3);
            n = n + 1;
        end
    end
end
%% 提取冬季数据
time_str = num2str(in_data(:,1));
r = 1;
for i = 2 : length(in_data)
    if str2num(time_str(i,5:6)) == 12 | str2num(time_str(i,5:6)) == 1 | str2num(time_str(i,5:6)) == 2
        data(r,:) = in_data(i,:);
        r = r + 1;
    end
end
        
%% 冬季降水距平
for k = 1 : 41
    pre(k,1) = 1977 + k;
    pre(k,2:88) = sum(data((1 + 3*(k-1)):(3 + 3*(k-1)),2:end),'omitnan');%各年9-11月（3个月）总降水量
end
pre(1,:)=[]; % 删除1978年数据
pre_m = mean(pre(:,2:88),'omitnan');
for i = 1 : 40
    for j = 2 : 88
        jupin(i,1) = pre(i,1);
        jupin(i,j) = pre(i,j)-pre_m(1,j-1);
    end
end


%eof第一模态图用变量eof的第一列，以此类推；相应的时间系数用pc第一行
x = jupin(:,2:end)';                  %原始数据，列为站点或格点，行为时间序列上的值
c = x*x'/87; %协方差矩阵
[eof,eigenvalue] = eig(c);% eof为特征向量，eigenvalue为特征值
pc = eof' * x;% pc为主成分，即时间系数
eigenvalue = rot90(eigenvalue,2); % flipud矩阵上下翻转；fliplr矩阵左右翻转
lambda = diag(eigenvalue); %lambda(λ)
eof = fliplr(eof);
pc = flipud(pc);
draw_eof = [jwd,eof]; %绘图用

% check
I2 = eof * eof';    % = I 
lambda2 = pc * pc';   % = lambda
x2 = eof * pc;     % = x

%计算方差贡献
r(1,1) = lambda(1,1)/sum(lambda(:,1));
r(1,2) = lambda(2,1)/sum(lambda(:,1));
r(1,3) = lambda(3,1)/sum(lambda(:,1));
r(1,4) = lambda(4,1)/sum(lambda(:,1));
r(1,5) = lambda(5,1)/sum(lambda(:,1));
disp(['总方差贡献率',num2str(sum(r))])

%时间系数图
for i = 1:40
    time_xishu(1,i) = mean(pc(1,i));
end
plot(1979:2018,pc(1,1:40),'r')
hold on
plot(1979:2018,pc(2,1:40),'b')
plot(1979:2018,pc(3,1:40),'g')
legend('第一模态','第二模态','第三模态')

有问题联系我，备注CSDN