eof.m

提取变量场的主要成分

filename=textread('I:\SODA\filename.dat','%21c');%filename是字符串数组，存放SODA资料nc文件的文件名。
[di dj]=size(filename);%di是所有文件的个数，对应于月份的个数。

% Pacific region
nlat=122;
nlon=303;
for ii=1:di 
  f=netcdf(strcat('I:\SODA\',filename(ii,:)),'nowrite');
  temp=f{'temp'}(1:1,91:212,240:542);
  temp(abs(temp)>100)=NaN;
  sst(:,:,ii)=squeeze(temp);
end
clear temp

%  mean monthly sst 
for t=1:12
    avsst(:,:,t)=mean(sst(:,:,t:12:t+12*39),3);
end


% monthly anomaly annaul mean
for t=1:50
    asst(:,:,1+(t-1)*12:12+(t-1)*12)=sst(:,:,1+(t-1)*12:12+(t-1)*12)-avsst(:,:,:);
end
clear sst
clear avsst



% quanzhong
lat=-30.25*pi/180:0.5*pi/180:30.25*pi/180;
quan(1:122)=sqrt(cos(lat));

for la=1:122
    ansst(la,:,:)=asst(la,:,:)*quan(la);
end
clear quan 
clear asst


% convert to 2 dimensional data
t=600;
k=0;kkk=0;
for j=1:nlat
    for i=1:nlon
        
        kk=0;
        for m=1:t
            if isnan(ansst(j,i,m))
                kk=kk+1;
            end
        end
        
        if (kk~=0) 
            land(j,i)=1;
            kkk=kkk+1;
        else
            land(j,i)=0;
            k=k+1;
            for m=1:t
                z(k,m)=ansst(j,i,m);
            end
        end
        
    end
end
clear ansst

%------------------- deviation ----------------------------
mz=mean(z,2);
for i=1:t
    zd(:,i)=z(:,i)-mz;
end
clear mz
%------------------- eof ---------------------------------
if k<=m 
    c=zd*zd';
    [v,d]=eig(c);
    x=diag(d);
    [b,ix]=sort(x);
    b=flipud(b);ix=flipud(ix);
    for i=1:length(ix)
        v2(:,i)=v(:,ix(i));
    end
    per=b/sum(b);
    ts=v2'*zd;
else
    c=zd'*zd;
    [v,d]=eig(c);
    x=diag(d);
    [b,ix]=sort(x);
    b=flipud(b);ix=flipud(ix);
    l=z*v;
    clear zd
    clear z 
    clear v
    for i=1:t
        lr(:,i)=l(:,i)/sqrt(x(i));
    end
    %lr=(z*v)/sqrt(x);
    for i=1:length(ix)
        v2(:,i)=lr(:,ix(i));
    end
    clear l
    per=b/sum(b);
    ts=v2'*zd;
end
disp('end');



sss=nan;kmod=6;

for m=1:kmod
    kk=0;
    for j=1:nlat
        for i=1:nlon
            hcurl(j,i,m)=sss;
            if(land(j,i)==0)
                kk=kk+1;
                hcurl(j,i,m)=v2(kk,m);
            end
        end
    end
end
clear v2
 figure;
 h=hcurl(:,:,1);
contourf(squeeze(h))

axis equal
colorbar
figure
plot(ts(2,:))
save I:\mysodadata\variancecontribiction.dat per -ascii

for nmod=1:6
    time(nmod,:)=(ts(nmod,:)-mean(ts(nmod,:)))/std(ts(nmod,:));
end

for nmod=1:6
    vector(:,:,nmod)=hcurl(:,:,nmod)*std(ts(nmod,:));
end






fid=fopen('I:\mysodadata\time.dat','wb');
 for t=1:600
  fwrite(fid,time(1:1,t),'float32');
  fwrite(fid,time(2:2,t),'float32');
  fwrite(fid,time(3:3,t),'float32');
  fwrite(fid,time(4:4,t),'float32');
  fwrite(fid,time(5:5,t),'float32');
  fwrite(fid,time(6:6,t),'float32');
 end
fclose(fid); clear fid ;


vecter(isnan(vector))=10000;

fid=fopen('I:\mysodadata\vector.dat','wb');
  fwrite(fid,(squeeze(vector(:,:,1:1)))','float32');
  fwrite(fid,(squeeze(vector(:,:,2:2)))','float32');
  fwrite(fid,(squeeze(vector(:,:,3:3)))','float32');
  fwrite(fid,(squeeze(vector(:,:,4:4)))','float32');
  fwrite(fid,(squeeze(vector(:,:,5:5)))','float32');
  fwrite(fid,(squeeze(vector(:,:,6:6)))','float32');
fclose(fid); clear fid ;