EOF

liuxiuqing

clc;clear;close all
% 参考书籍
% 黄嘉佑《气象统计分析与预报方法》（第四版）
% 魏凤英《现代气候统计诊断与预测技术》（第2版）
% matlab eig 函数：https://ww2.mathworks.cn/help/matlab/ref/eig.html
% matlab M_Map 工具箱：https://www.eoas.ubc.ca/~rich/map.html
% matlab M_Map 工具箱工具箱的安装：http://blog.sina.com.cn/s/blog_8fc890a20102v6pm.html
%% 数据读取
% ncdisp('sst-195001-201712.nc');             % NOAA Extended Reconstructed SST V4
sst=ncread('sst-195001-201712.nc','sst');     % 海温
lon=ncread('sst-195001-201712.nc','lon');     % 经度，分辨率：2°
lat=ncread('sst-195001-201712.nc','lat');     % 纬度，分辨率：2°
time=ncread('sst-195001-201712.nc','time');   % 时间，月平均，68年，12个月
%% 选区域
% 框选热带太平洋地区(Tropical Pacific：120°E-70°W；30°S-30°N）
sst_TP=sst(61:146,30:60,:);
% 把数据展开为四维（经度、纬度、月、年）
sst_TP_12=reshape(sst_TP,86,31,12,68);
[lon_TP,lat_TP,month_TP,year_TP]=size(sst_TP_12);
%% 资料预处理
% 求平均(每个格点不同月份的年平均)
average=squeeze(nanmean(sst_TP_12(:,:,:,1:year_TP),4));  % 68年
% 算距平（减去对应月份的平均值）
sst_TP_anomaly_12=sst_TP_12-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);  % 《气象统计分析与预报方法》142页
%% 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'的特征值（《现代气候统计诊断与预测技术》109页）
    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); % 特征值的取样误差（North et al.,1982）
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]);
%% 第一模态画图
load colormap.mat
figure ('colormap',colormap)
% 空间分布
subplot(2,1,1);
V_1=V_NaN(:,:,1);  
m_proj('miller','lon',[120 290],'lat',[-30 30]);
[X1,Y1]=meshgrid(120:2:290,-30:2:30);
m_contourf(X1,Y1,rot90(V_1),'linestyle','none');
hold on
m_contour(X1,Y1,rot90(V_1),[-0.1 0.1 0.4 0.8],'linecolor','k','ShowText','on');
m_coast('patch',[.85 .85 .85 ],'edgecolor',[0.6 0.6 0.6]);
m_grid('linestyle','none','xtick',[140 180 220 260],'ytick',[30 15 0 -15 -30],'fontsize',14,'tickdir','out');
caxis([-1,1])
ax=colorbar;
set(ax,'linewidth',0.5,'fontsize',12,'edgecolor','k','YTick',[-1,-0.5,0,0.5,1],'YTickLabel',{'-1','-0.5','0','0.5','1'})
set(gcf,'color','white','Position',[100,100,1000,500]);
title(['(a) EOF1 = ' num2str(R(1:1)*100,4) '%'],'fontsize',15,'color','b');
% 时间系数
subplot(2,1,2);
T_1=T(1,:);
t=[1:month_TP*year_TP];
plot(t,T_1,'k','linewidth',1);
hold on;
T_0=zeros(1,month_TP*year_TP);
plot(t,T_0,'k--','linewidth',1);
axis([0 840 -4 4]);
set(gca,'YTick',-4:1:4,'tickdir','out','yminortick','on');
set(gca,'YTicklabel',{'-4','-3','-2','-1','0','1','2','3','4'});
set(gca,'XTick',0:60:840,'tickdir','out','xminortick','off');
set(gca,'XTicklabel',{'1950','1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010','2015','2020'});
set(gca,'linewidth',1,'fontsize',15,'fontname','Times');
set(gcf,'color','white','Position',[100,100,1000,500]);
title('(b) PC1','fontsize',15,'color','b');
%% 第二模态画图
load colormap.mat
figure ('colormap',colormap)
% 空间分布
subplot(2,1,1);
V_2=V_NaN(:,:,2);  
m_proj('miller','lon',[120 290],'lat',[-30 30]);
[X1,Y1]=meshgrid(120:2:290,-30:2:30);
m_contourf(X1,Y1,rot90(V_2),'linestyle','none');
hold on
m_contour(X1,Y1,rot90(V_2),[-0.1 0.1 0.4 0.8],'linecolor','k','ShowText','on');
m_coast('patch',[.85 .85 .85 ],'edgecolor',[0.6 0.6 0.6]);
m_grid('linestyle','none','xtick',[140 180 220 260],'ytick',[30 15 0 -15 -30],'fontsize',14,'tickdir','out');
caxis([-1,1])
ax=colorbar;
set(ax,'linewidth',0.5,'fontsize',12,'edgecolor','k','YTick',[-1,-0.5,0,0.5,1],'YTickLabel',{'-1','-0.5','0','0.5','1'})
set(gcf,'color','white','Position',[100,100,1000,500]);
title(['(a) EOF2 = ' num2str(R(2:2)*100,4) '%'],'fontsize',15,'color','b');
% 时间系数
subplot(2,1,2);
T_2=T(2,:);
t=[1:month_TP*year_TP];
plot(t,T_2,'k','linewidth',1);
hold on;
T_0=zeros(1,month_TP*year_TP);
plot(t,T_0,'k--','linewidth',1);
axis([0 840 -4 4]);
set(gca,'YTick',-4:1:4,'tickdir','out','yminortick','on');
set(gca,'YTicklabel',{'-4','-3','-2','-1','0','1','2','3','4'});
set(gca,'XTick',0:60:840,'tickdir','out','xminortick','off');
set(gca,'XTicklabel',{'1950','1955','1960','1965','1970','1975','1980','1985','1990','1995','2000','2005','2010','2015','2020'});
set(gca,'linewidth',1,'fontsize',15,'fontname','Times');
set(gcf,'color','white','Position',[100,100,1000,500]);
title('(b) PC2','fontsize',15,'color','b');

dengbin313

厉害，学习了

wsjqxjy

{:eb502:}{:eb502:}{:eb502:}