Matlab--EOF分析

%% EOF分析
clc;clear;close all

%% 读取降水相关数据
datadir = 'E:\prep\';
% 循环读取降水数据
filelist=dir([datadir,'*.nc']);     
for  i = 1:129          
    ncid =[datadir,filelist(i).name];    
    time(:,i)=ncread(ncid,'time');    
    precip(:,:,i)=ncread(ncid,'precip');
end
path = [datadir,filelist(1).name];
ncdisp(path)
mlat = double(ncread(path,'latitude'));
mlon = double(ncread(path,'longitude'));
[nlat,nlon]=meshgrid(mlat,mlon);

%% 读取shp数据
cdL=shaperead("E:\世界地图shp\世界地图2017\2017世界地图.shp");
cx=[cdL(:).X]; cy=[cdL(:).Y];

%% 计算
precip=(reshape(precip,20*14,129))';     % 将三维数据转化为二维数据

F=detrend(precip,0);                     % 移除均值（0阶多项式趋势）
C=F'*F;                                   % 计算协方差矩阵C

[EOFs,D]=eig(C);                         % EOF特征向量，D特征值
PCs=EOFs'*F';                            % 时间权重PC

D=diag(D);                                % 提取对角线元素得到特征值序列
D=D./sum(D);                             % 特征值标准化

%% 第一模态
EOF1=EOFs(:,D==max(D));                  % 选取特征值最大的作为第一模态
PC1=PCs(D==max(D),:);

EOF1=reshape(EOF1,20,14);                % 重塑为20*14的空间数据

EOF1=EOF1.*std(PC1);                     % 标准化
PC1=PC1./std(PC1);
disp(max(D)*100)                         % 显示方差解释度

%% 第二模态
EOF2=EOFs(:,279);                        % 选取特征值第二大的作为第二模态
PC2=PCs(279,:);

EOF2=reshape(EOF2,20,14);

EOF2=EOF2.*std(PC2);
PC2=PC2./std(PC2);
disp(D(279,1)*100)

%% 画图
figure('position',[300,10,1200,800])

% 第一模态空间分布
subplot(2,2,1);
m_proj('miller','lon',[min(mlon) max(mlon)],'lat',[min(mlat) max(mlat)]);
m_contourf(nlon,nlat,EOF1,'linestyle','none');
m_grid('linestyle','none','xtick',([108 120 132 144]));
hold on
m_plot(cx,cy,'k','linewidth',2); 
colormap(m_colmap('diverging'));
s=colorbar('fontname','Times New Roman','fontsize',12);
title(s,'mm/day','fontname','Times New Roman');
set(gca,'fontsize',12)
title('EOF1: 29.88%','fontsize',16,'fontname','Times New Roman');

% 第一模态时间序列
subplot(2,2,3);
plot(1:129,PC1,'r','linewidth',2);
hold on
plot(xlim,[0,0],'k--','LineWidth',1.5)
set(gca,'xtick',[2:30:129],'xticklabels',1979:10:2021, ...
    'fontname','Times New Roman','fontsize',12);
xlabel('Year','fontname','Times New Roman');
ylabel('PC1','fontname','Times New Roman');
xlim([1 129]);
% ylim([-2.5 2.5]);
set(gca,'fontsize',12,'linewidth',2)
title('PC1: 29.88%','fontsize',16,'fontname','Times New Roman');

% 第二模态空间分布
subplot(2,2,2);
m_proj('miller','lon',[min(mlon) max(mlon)],'lat',[min(mlat) max(mlat)]);
m_contourf(nlon,nlat,EOF2,'linestyle','none');
m_grid('linestyle','none','xtick',([108 120 132 144]));
hold on
m_plot(cx,cy,'k','linewidth',2); 
colormap(m_colmap('diverging'));
s=colorbar('fontname','Times New Roman','fontsize',12);
title(s,'mm/day','fontname','Times New Roman');
set(gca,'fontsize',12)
title('EOF2: 12.72%','fontsize',16,'fontname','Times New Roman');

% 第二模态时间序列
subplot(2,2,4);
plot(1:129,PC2,'r','linewidth',2);
hold on
plot(xlim,[0,0],'k--','LineWidth',1.5)
set(gca,'xtick',[2:30:129],'xticklabels',1979:10:2021, ...
    'fontname','Times New Roman','fontsize',12);
xlabel('Year','fontname','Times New Roman');
ylabel('PC2','fontname','Times New Roman');
xlim([1 129]);
set(gca,'fontsize',12,'linewidth',2)
title('PC2: 12.72%','fontsize',16,'fontname','Times New Roman');