记录帖：
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); %转换为二维数组
EOF分析
©著作权归作者所有,转载或内容合作请联系作者
- 文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
- 文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
- 文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
推荐阅读更多精彩内容
- 写了个程序，对Numpy的绝大部分函数及其说明进行了中文翻译。原网址：https://docs.scipy.org...
