收藏本站网站公告 | 吾爱海洋论坛交流QQ群：835383472
资迅
资迅
升学
升学
工作
工作
吹水
吹水
吾爱海洋»论坛 › 科研交流区 › 数据处理软件 › 如何实现Matlab主成分分析？

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

[复制链接]
關錕嶺 楼主 发表于 数据处理软件 2022-3-1 21:40:11 | 显示全部楼层 阅读 818

本文作者：中国海洋大学，华睿

联系邮箱：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


8 x- o1 m7 z, H8 D( l

[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);

! G" n8 S5 k' S0 M$ a" ]: r

17b48b0a5032e7f8f17a9c4a34186d9d.png

4fd704fcd3691a27555795bce3affdb4.png

) n/ t+ t7 j9 E

[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);


& y+ Y- h+ v6 J- ?0 c1 R$ K

7e050f16a5013b8f2b763a9828b24fa7.png

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

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

如何 实现 matlab 主成分分析 成分分析
回复

举报 使用道具

上一篇: 科研向——海洋科学（物理海洋）专业可视化软件综述 ...

下一篇: 用4行代码在matlab地图上绘制省界&国界

相关帖子

全部回帖
梨花白 | 显示全部楼层
感谢楼主的分享，楼主身体好
发表于 2022-8-8 15:30:13

举报 回复 使用道具

返回列表 发新帖
懒得打字？点击右侧快捷回复 【吾爱海洋论坛发文有奖】
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册
關錕嶺
活跃在2024-12-1

手机版|吾爱海洋 ( 粤ICP备2020130938号-2 )|友链申请|网站地图

Copyright © 2001-2024, 吾爱海洋.　Powered by 吾爱海洋团队 X3.4

快速回复 返回顶部 返回列表