- 积分
- 107
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-9-5
- 最后登录
- 1970-1-1
|
MATLAB
| 问题截图: |
- |
| 问题概况: |
在做大西洋sst的EOF分解，算贡献率时出错了，求大神指点！
%清空先前变量，关闭窗口
clear; clc;
close all;
%读入数据
filename='HadISST_sst.nc';
ncdisp(filename);
sst1=ncread(filename,'sst');
lat=ncread(filename,'latitude');
lon=ncread(filename,'longitude');
%数据选取，1991——2010年，20°N——80°N，90°W——40°E
sst=sst1(90:221,10:71,1441:1680);
%将数据转化为二维数组
sst_a=zeros(240,8184);
for i=1,240
sq=squeeze(sst(:,:,i));
sst_a(i,:)=reshape(sq,1,8184);
end
%选取有用(海水)数据
sst_a(sst_a<-100)=NaN;
i=1;
for j=1:8184
if sst_a(i,j)==-1.79999995231628
sst_a(i,j)=NaN; %冰点填充值为-1.79999995231628
i=i+1;
end
end
sst_b=isnan(sst_a);
i=0;
for j=1:8184
if sum(sst_b(:,j))==0;
i=i+1;
sst_c(:,i)=sst_a(:,j);
end
end
%------------------------------------------------
%求距平值
X=zeros(size(sst_c));
for i=1:12;
X(i:12:228+i,:)=sst_c(i:12:end,:)-repmat(mean(sst_c(i:12:end,:),1),size(sst_c(i:12:end,:),1),1);
end
%------------------------------------------------------------
%时空转换
R=X*X';
[v,d]=eig(R);%d为x的特征值组成的对角阵，v为X*X'的特征向量
v=fliplr(v); %矩阵左右翻转
d=rot90(d,2);%矩阵上下翻转后左右翻转
diagonal=zeros(size())
diagonal=diag(d);
spacef=X'*v;
for i=1:240;
spacef(:,i)=spacef(:,i)/sqrt(diagonal(i));%空间本征函数
end
timef=X*spacef; %时间本征函数
sum_d=sum(diagonal);
count=0;
for i=1:10;
count=count+diagonal(i);
G1(i)=count/sum_d; %G1（i）为累积方差贡献率
end
for i=1,10;
G2(i)=diagonal(i)/sum_d;%G2(i)为方差贡献率
end
%--------------------------------------------------
%将去除数据补回
sst_a2=zeros(240,8184);
sst_a2(:,:)=NaN;
i=1;
spacef2=spacef';
for j=1,8184
if sum(sst_b(:,j))==0;
sst_a2(:,j)=spacef2(:,j);
i=i+1;
end
end
sst_a3=sst_a2';
%--------------------------------------------
%画图 |
| 我看过提问的智慧: |
看过 |
| 自己思考时长(天): |
2 |
| 系统平台: |
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
在线等
|
|
窥视卡
雷达卡
发表于 2014-4-21 20:25:57
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
发表于 2014-4-22 12:54:53
楼主