[C] 纯文本查看 复制代码
clc,clear,close all;fclose all;
%% 读入原始的GPCP数据
ncdisp('F:/RMatlab/data20211219/precip.mon.mean.nc');
x=ncread('F:/RMatlab/data20211219/precip.mon.mean.nc','lon');
x=double(x);
y=ncread('F:/RMatlab/data20211219/precip.mon.mean.nc','lat');
y=double(y);
gpcp=ncread('F:/RMatlab/data20211219/precip.mon.mean.nc','precip');
%%
gpcp(:,:,517:end)=[]; % % 把末尾数据(2021年)给删掉
gpcp=reshape(gpcp,length(x),length(y),12,[]); % % 经度*纬度*月份*年份
gpcp=gpcp(x>=107 & x<=122 , y>=20 & y<=43, : , :); % % 107-122E, 20-43N
gpcp=squeeze(mean(mean(gpcp(:,:,6:8,:),3),1)); % % JJA平均
for ii=1:9
gpcp(ii,:)=zscore(squeeze(gpcp(ii,:))); % % 标准化变量
end
runEOF1=nan(9,42); % % 最后画图的数据
for ii=8:35
XX=gpcp(:,ii-7:ii+7); % % 纬度*15年滑动
[v,d]=eig(XX*XX');
lam=real(diag(d));
[B,idx]=sort(lam,'descend');
VV=v(:,idx(1:5));
PC=VV'*XX;
PC=PC'; % % EOF的计算过程
PC1=PC(:,1); % % PC1
EOF1=VV(:,1)*std(PC1); % % EOF1
runEOF1(:,ii)=EOF1; % % 滑动EOF1
end
%% 但是这里有个问题,滑动EOF的符号可能突变,画图就不好看了
for ii=8:34
if dot(runEOF1(:,ii),runEOF1(:,ii+1))<0
runEOF1(:,ii+1)=-1*runEOF1(:,ii+1);
% % 如果后一年符号有突变,就乘以个-1,换回来
end
end
%% 画图看看
figure,hold on;
contourf(1979:2020,21.25:2.5:41.25,runEOF1,...
-1.5:0.01:1.5,'linestyle','none');
colormap('jet');colorbar;caxis([-1 1]);
axis([1986 2013 21.25 41.25]);
c=colorbar('eastoutside','ticklength',0);
caxis([-1,1])
% 调用ncl色带
load('colorbar-mat/BkBlAqGrYeOrReViWh200_r.mat');
colormap(BkBlAqGrYeOrReViWh200_r);
ax = gca;
axpos = ax.Position;
c.Position(3) = 0.5*c.Position(3);
ax.Position = axpos;
cbarrow;