matlab eof分析

神秘的宝宝

what m_map
rehash toolboxcache
ncinfo('HadISST_sst.nc')
sst1=ncread('HadISST_sst.nc','sst');     % 海温
lat1=ncread('HadISST_sst.nc','latitude');     % 纬度
lon1=ncread('HadISST_sst.nc','longitude');     % 经度
time1=ncread('HadISST_sst.nc','time');   % 时间
%% 选区域
% 框选南太平洋地区(Tropical Pacific 30S-90S，130E-60W）
sst_TP1=sst1(1:120,120:180,1:1800);
sst_TP2=sst1(310:360,120:180,1:1800);
sst_TP12=[sst_TP1;sst_TP2]
% 把数据展开为四维（经度、纬度、月、年）
sst_TP_122=reshape(sst_TP12,171,61,12,150);
[lon_TP,lat_TP,month_TP,year_TP]=size(sst_TP_122);
%% 资料预处理
% 求平均(每个格点不同月份的年平均)
average=squeeze(nanmean(sst_TP_122(:,:,:,1:year_TP),4));  % 150年
% 算距平（减去对应月份的平均值）
sst_TP_anomaly_12=sst_TP_122-repmat(average,[1,1,1,year_TP]);
% 降维，排列成m×n的形式（m：站点；n：时间）
sst_TP_anomaly=reshape(sst_TP_anomaly_12,lon_TP,lat_TP,month_TP*year_TP);
X_NaN=reshape(sst_TP_anomaly,lon_TP*lat_TP,month_TP*year_TP);
% 去除NaN数据（NaN为陆地，没有海温）
X_logical=isnan(X_NaN);  % 逻辑判断
i=0;
for j=1:lon_TP*lat_TP
    if sum(X_logical(j,:))==0
       i=i+1;
       X(i,:)=X_NaN(j,:);
    end
end
X=X/sqrt(month_TP*year_TP);
%% EOF运算
[X_x,X_y]=size(X);
if X_x<=X_y             % 不考虑时空转换
   
    S=X*X';             % 计算协方差
    [v,d]=eig(S);       % 使用eig函数求特征向量v和特征值d
   
    % 使用sort函数将特征值、特征向量按升序排序
    [~,ind]=sort(diag(d));
    Ds=d(ind,ind);
    Vs=v(:,ind);
   
    V=fliplr(Vs);  % 排序后的特征向量V
    D=rot90(Ds,2); % 排序后的特征值D
else                     % 考虑时空转换
    S=X'*X;              % 计算协方差
    [v,d]=eig(S);        % 使用eig函数求特征向量v和特征值d
    % 使用sort函数将特征值、特征向量按升序排序
    [~,ind]=sort(diag(d));
    Ds=d(ind,ind);
    Vs=v(:,ind);
   
    VR=fliplr(Vs);  % 排序后的特征向量V
    D=rot90(Ds,2);  % 排序后的特征值D
   
    % 求X*X'的特征值
    VR=X*VR;        
    for i=1:length(D)
        V(:,i)=VR(:,i)/sqrt(D(i,i));
    end
end
% 时间系数
T=V'*X*sqrt(month_TP*year_TP);
% 特征向量的方差贡献
diagonal=diag(D);
total=0;
for i=1:length(diagonal)
    R(i)=diagonal(i)/sum(diagonal);  % 方差贡献率
    total=total+diagonal(i);
    G(i)=total/sum(diagonal);        % 累积方差贡献率
end
% 标准化处理
for i=1:length(diagonal)
    V(:,i)=V(:,i)*sqrt(D(i,i));    % 特征向量乘以特征值的平方根
    T(i,:)=T(i,:)/sqrt(D(i,i));    % 时间系数除以特征值的平方根
end
%% 把删去的NaN数据补回
[a,b]=size(V);
V_NaN=NaN(lon_TP*lat_TP,b);
j=1;
for i=1:lon_TP*lat_TP
    if sum(X_logical(i,:))==0
       V_NaN(i,:)=V(j,:);
       j=j+1;
    end
end
%% 升维
V_NaN=reshape(V_NaN,lon_TP,lat_TP,b);
%% 显著性检验
figure
err=sqrt(2/(month_TP*year_TP))*diagonal/sum(diagonal); % 特征值的取样误差
dot=1:4;
R_1_4=R(1:4)*100;
err=err(1:4)'*100;
plot(dot,R_1_4,'.','Markersize',12);
hold on
errorbar(dot,R_1_4,err,'LineStyle','none');
axis([0 5 0 50]);
set(gca,'YTick',0:10:50,'tickdir','out','yminortick','off');
set(gca,'YTicklabel',{'0','10','20','30','40','50'});
set(gca,'XTick',0:1:5,'tickdir','out','xminortick','off');
set(gca,'XTicklabel',{'0','1','2','3','4','5'});
set(gca,'linewidth',1,'fontsize',15,'fontname','Times');
set(gcf,'color','white','Position',[100,100,1000,500]);
title('Percentage variance (%)','fontsize',15,'color','b','position',[2.5,52]);