[Matlab] 如何实现Matlab主成分分析?

[复制链接]

本文作者:中国海洋大学,华睿

联系邮箱:958660470@qq.com


本文数据来源:

位势高度场,风场等月平均数据:

https://psl.noaa.gov/data/gridded/data.ncep.reanalysis.pressure.html

海表面温度(SST)月平均数据:

https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html

降水数据(CMAP月平均,Enhanced Monthly Mean):

https://psl.noaa.gov/data/gridded/data.cmap.html

6 b- p  m  l" _+ c, T) [

[C++] 纯文本查看 复制代码
clc,clear,close all;fclose all;
figure('color',[1,1,1])
lon=ncread('F:/RMatlab/data20211219/sst.mnmean.nc','lon');
lat=ncread('F:/RMatlab/data20211219/sst.mnmean.nc','lat');
time=ncread('F:/RMatlab/data20211219/sst.mnmean.nc','time');
ntimes = length(time);
nlats=length(lon);nlons=length(lat);
sst_data=[nlats, nlons, ntimes];
lon1=-30;lon2=30;
[anod] = region(sst_data, lat, lon, [lon1 lon2], [-250, -60]+360 );
m1=size(anod,1);m2=size(anod,2);n=size(anod,3);
%----
for i=1:15   
    i1=(16-i)*2;%---需要调整一下个数
    anod(:,i,:)=cosd(i1)*anod(:,i,:);
    anod(:,15+i,:)=cosd(i1)*anod(:,15+i,:);
end
x=reshape(anod,[m1*m2,n]);
r=2;
%---处理缺省值
a=find(isnan(x(:,:))==0);
c=find(~isnan(x(:,:))==0);
c0=c(1:393);
x(c0,:)=0;
[e,p,lmd,vf,evf]=EOF_2D(x,r);
%--NO.1
x1=reshape(e(:,r),[m1,m2]);%---修改主成分的空间系数
x2=reshape(e(:,r-1),[m1,m2]);%---修改主成分的空间系数
for i=1:15   
    i1=(16-i)*2;%---需要调整回来
    x1(:,i)=x1(:,i)/cosd(i1);
    x1(:,15+i)=x1(:,15+i)/cosd(i1);
end
p1=p(r,:);%---修改主成分的时间系数
p2=p(r-1,:);%---修改主成分的时间系数
pp1(1:40)=0;
plot(1980:2019,pp1,'linewidth',1,'color','r');
hold on
plot(1980:2019,p2,'linewidth',1.5,'color','b');
xlabel('Year','FontSize',10,'FontWeight','bold','Color','k');
set(gca,'YLim',[-3 3]);
set(gca,'FontSize',10,'FontName','Times New Roman','FontWeight','Bold') 
title('PC-2');
figure('color',[1,1,1])
boundary=[110 300 -30 30];
lon_scope = find(lon >= boundary(1) & lon<=boundary(2));
lat_scope = find(lat >= boundary(3) & lat<=boundary(4));
lon_number = length(lon_scope);     
lat_number = length(lat_scope);
lat_1=linspace(boundary(3),boundary(4),lat_number);
lon_1=linspace(boundary(1),boundary(2),lon_number); 
[plon,plat]=meshgrid(lon_1,lat_1);
guojie=shaperead('bou2_4l.shp');
m_proj('Equidistant','long', [110 180],'lat',[-30,30]);
bou1_4lx=[guojie(:).X];
bou1_4ly=[guojie(:).Y];
m_plot(bou1_4lx,bou1_4ly,'linewidth',2,'color',[0.35 0.35 0.35]);
hold on
eof1_plot=imrotate(x1,90);
eof2_plot=imrotate(x2,90);
% [cs,h]=m_contourf(plon,plat,eof2_plot,'linewidth',0.5);%---估计EOF
[cs,h]=m_contourf(plon,plat,p2(40)*eof2_plot+p1(40)*eof1_plot,'linewidth',0.5);%%----(3)利用EOF的前2模态对2019年12月的海温距平场做出估计(绘图)
D2019=anod(:,:,40);D2019_plot=imrotate(D2019,90);
% [cs,h]=m_contourf(plon,plat,D2019_plot,'linewidth',0.5);%%----(3)实际测量的2019年12月ssta
set(h,'LineColor','none')%取消等值线
set(h,'LevelStep',get(h,'LevelStep')*0.2);
m_coast('color',[0 0 0],'linewidth',1);
hold on
str=sprintf( 'EOF-2');
% str=sprintf( 'SST anomaly for December 2019');
str=sprintf( 'SST anomaly estimates for December 2019(EOF-1,2)');
title(str,'fontsize',10,'color','k')
set(gca,'FontSize',10,'FontName','Times New Roman','FontWeight','Bold') 
% m_grid('xtick',6,'ytick',[-30 -15 0 15 30],'fontsize',18);
m_grid('ytick',[-30:10:30],'xtick',[110:10:180],'tickdir','out','linest','none','fontname','Times','fontsize',12,'linewidth',1.5);
% ---设置经纬度,m_grid.m第400和500多行
c=colorbar('eastoutside','ticklength',0);
caxis([-0.8,0.8])
hold on
% 调用ncl色带
load('colorbar-mat/BkBlAqGrYeOrReViWh200_r.mat');
CMap1=flipud(BkBlAqGrYeOrReViWh200_r); 
colormap(CMap1);
ax = gca;
axpos = ax.Position;
c.Position(3) = 0.5*c.Position(3);
ax.Position = axpos;
cbarrow;
e1=vf(2)*sqrt(2/40);
e2=vf(1)*sqrt(2/40);

! ]" Z& }* j' K! R" u

17b48b0a5032e7f8f17a9c4a34186d9d.png

4fd704fcd3691a27555795bce3affdb4.png


/ j, x" S& a+ N( C

[C] 纯文本查看 复制代码
%---用主成分时间序列进行回归分析
clc,clear,close all;fclose all;
data='F:/RMatlab/data20211219/sst.mnmean.nc';
ncdisp(data);
ncinfo(data);
m_proj('Equidistant','long', [110 300],'lat',[-30,30]);
boundary=[0 360 -90 90];
lon=ncread('F:/RMatlab/data20211219/vwnd.mon.mean.nc','lon');
lat=ncread('F:/RMatlab/data20211219/vwnd.mon.mean.nc','lat');
lat1=ncread('F:/RMatlab/data20211219/precip.mon.mean.nc','lat');
lon_scope = find(lon >= boundary(1) & lon<=boundary(2));
lat_scope = find(lat >= boundary(3) & lat<=boundary(4));
lon_number = length(lon_scope);     
lat_number = length(lat_scope);
lat_scope1 = find(lat1 >= boundary(3) & lat1<=boundary(4));
lat_number1 = length(lat_scope1);
u=ncread('F:/RMatlab/data20211219/uwnd.mon.mean.nc','uwnd');
v=ncread('F:/RMatlab/data20211219/vwnd.mon.mean.nc','vwnd');
pre=ncread('F:/RMatlab/data20211219/precip.mon.mean.nc','precip');
level=ncread('F:/RMatlab/data20211219/uwnd.mon.mean.nc','level');
ud=u(:,:,3,396:12:865);%---1979年1月为第一维度
vd=v(:,:,3,396:12:865);
pred=pre(:,:,24:12:493);
ud=ud(:,:,:);vd=vd(:,:,:);
% pred0=reshape(pred,[144*40,72]);
% b=ones(144*40,1);
% pred0 = cat(2,b,pred0);
% pred=reshape(pred0,[144,73,40]);
% ----一元回归
for i=1:lon_number
    for j=1:lat_number1
        a=pred(i,j,:);a=a(:);
        [b] = regress(a,p2') ;
        reg(i,j)=b;
    end
end 
for i=1:lon_number
    for j=1:lat_number
        u=ud(i,j,:);u=u(:);
        v=vd(i,j,:);v=v(:);
        [b1] = regress(u,p2') ;
        [b2] = regress(v,p2') ;
        reg1(i,j)=30*b1;
        reg2(i,j)=30*b2;      
    end
end 
lat_1=linspace(boundary(3),boundary(4),lat_number);
lat1_1=linspace(boundary(3),boundary(4),lat_number1);
lon_1=linspace(boundary(1),boundary(2),lon_number); 
[plon,plat]=meshgrid(lon_1,lat_1);
[plon1,plat1]=meshgrid(lon_1,lat1_1);
m_proj('Equidistant','long', [110 180],'lat',[-20,30]);
guojie=shaperead('bou2_4l.shp');
bou1_4lx=[guojie(:).X];
bou1_4ly=[guojie(:).Y];
m_plot(bou1_4lx,bou1_4ly,'linewidth',2,'color',[0.35 0.35 0.35]);
%-------更改这里(印度洋和太平洋)
reg_plot=imrotate(reg,90);
reg1_plot=imrotate(reg1,90);
reg2_plot=imrotate(reg2,90);
figure('color',[1,1,1])
[cs,h]=m_contourf(plon1,plat1,reg_plot,'linewidth',0.5);
hold on
m_quiver(plon(1:2:end),plat(1:2:end),reg1_plot(1:2:end),reg2_plot(1:2:end),2,'color','k')%%风场矢量图,'2'表示扩大多少为多少倍
set(h,'LineColor','none')%取消等值线
set(h,'LevelStep',get(h,'LevelStep')*0.2);
m_coast('color',[0 0 0],'linewidth',1);
hold on
str=sprintf( 'Regression of Precip and 850hPa wind on PC-2 from 1980 to 2019(Dec)\n');
title(str,'fontsize',10,'color','k')
set(gca,'FontSize',10,'FontName','Times New Roman','FontWeight','Bold') 
% m_grid('xtick',6,'ytick',[-30 -15 0 15 30],'fontsize',10);
m_grid('ytick',[-20:10:30],'xtick',[110:10:180],'tickdir','out','linest','none','fontname','Times','fontsize',12,'linewidth',1.5);
%---设置经纬度,m_grid.m第400和500多行
c=colorbar('eastoutside','ticklength',0);
caxis([-2,2])
hold on
% 调用ncl色带
load('colorbar-mat/BkBlAqGrYeOrReViWh200_r.mat');
CMap1=flipud(BkBlAqGrYeOrReViWh200_r); 
colormap(CMap1);
ax = gca;
axpos = ax.Position;
c.Position(3) = 0.5*c.Position(3);
ax.Position = axpos;
cbarrow;
[h1,h2,h3]=legend_vc('northeast',plon,plat,reg1_plot,reg2_plot,2,2);


' x' L! G- {& m  J; c3 U  E

7e050f16a5013b8f2b763a9828b24fa7.png


Matlab主成分分析数据和代码免费下载链接:

游客,如果您要查看本帖隐藏内容请回复

回复

举报 使用道具

相关帖子

全部回帖
感谢楼主的分享,楼主身体好
发表于 2022-8-8 15:30:13

举报 回复 使用道具

懒得打字?点击右侧快捷回复 【吾爱海洋论坛发文有奖】
您需要登录后才可以回帖 登录 | 立即注册
關錕嶺
活跃在2024-1-31
快速回复 返回顶部 返回列表