- 积分
- 9
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2020-3-25
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
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]);
|
|