EOF分析

记录帖：
clear
clc
%% 读取nc
lon=ncread('cru_ts4.03.1901.2018.pre.dat.nc','lon'); %经度
lat=ncread('cru_ts4.03.1901.2018.pre.dat.nc','lat'); %纬度
[X,Y]=meshgrid(lon,lat);
pre=ncread('cru_ts4.03.1901.2018.pre.dat.nc','pre'); %降水量
pp=reshape(pp,[720,360,12,83]);
p1=squeeze(nansum(pp,3));
%% 区域挑选
b=[35 55];
a=[45 90];
aa=find(lon<a(2)&lon>a(1));
bb=find(lat<b(2)&lat>b(1));
pre_region=p1(aa,bb,:);
pdata=reshape(pre_region,[3600 83]);
%% 距平
for i=1:3600
p_r(i,:)=pdata(i,:)-(nanmean(squeeze(pdata(i,:))));
end
%% 开始进行EOF分解
X=p_r;
ll=length(p_r);
R=XX'./83; % 构造空间协方差矩阵
[v,d]=eig(R); %求特征根和特征向量，每个非0的特征根对应一列特征向量值，也称为EOF。
%由于MATLAB分解出来的特征根是从小到大排列，所以，最后一个特征根是最大的，也就对应的是最后一特征向量。
EOF1=v(:,ll); %第一模态
%计算主成分：将EOF投影到原始资料矩阵X上，就得到所有空间特征向量对应的时间系数（即主成分）
PC=v'X; %PC是主成分，PC的每一行数据就是对应每个特征向量的时间系数
%求方差贡献
diagonal=diag(d);
sum_e=sum(diagonal);
G1=diagonal(ll)/sum_e; %方差贡献率
%% 模态转换为二维数组，还原。
u1=reshape(EOF1(1:3600),90,40); %转换为二维数组