- 积分
- 13
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2016-3-16
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
%ncdump
%regresstion analysis for enso
close all
clear
%read sst from hadisst
nc=netcdf('global sst.nc','r');
timer=nc{'TIME'}(:);
%base=datenum(1,1,1,0,0,0);
%time=timer+base;
%datestr(time/24);
lon=nc{'LON'}(:);
lat=nc{'LAT'}(:);
miss=nc{'SST'}.missing_value(:);
fill=nc{'SST'}.FillValue_(:);
sst=nc{'SST'}(:,:,:);%ncfloat('TIME', 'LAT', 'LON');time=360
sst(sst==fill)=nan;
sst(sst==miss)=nan;
sst(sst<-10)=nan; %min sst是 -1000,去除太小值
%%
%通常下载数据直接绘图是从-180 to 180 ;将地图经度换为 0 to 360
inn=find(lon<0);
inp=find(lon>=0);
sstn=sst(:,:,inn);
sstp=sst(:,:,inp);
sstr=cat(3,sstp,sstn);
lon1=lon(inn)+360;
lonr=[lon(inp)
lon1];
[long,latt]=meshgrid(lonr,lat);
%%
%扣除季节性变化
tmp=reshape(sstr,[12,30,180,360]);
season=squeeze(nanmean(tmp,2));
clim=repmat(season,[30,1,1]);
ssta=sstr-clim;
%load ninoindex
nino=load('nino index.txt');
enso=nino(1:360,5);%time 1982-2011
%plot(enso);
%regress
[tt,ii,jj]=size(ssta);
cof=nan*zeros(ii,jj);
cofs=nan*zeros(ii,jj);
reg=[ones(tt,1),enso];
for i=1:ii
for j=1:jj
tmp=squeeze(ssta(:,i,j));
%if~isnan(tmp)
[b,~,~,stats]=regress(tmp,reg);
%end
cof(i,j)=b(2);
cofs(i,j)=stats(3);
end
end
cof(cof==0)=nan;
%%
%画图
m_proj('Robinson','long',[min(lonr),max(lonr)],...
'lat',[min(lat)+5,max(lat)-5]);
%[X,Y]=m_112xy(lon,lat,'clip','on');
m_grid('box','off','tickdir','on',...
'fontsize',20);
hold on
m_contourf(long,latt,cof,'linestyle','none');
m_coast('patch',[.6 .6 .6]);
hold on
m_contour(long,latt,cofs,[0.05,0.05]);
colorbar;
myc1
caxis([-0.5,0.5]);
title('regresstion enso')
%问题,怎么把地图上的lable弄好点,或没有
|
|