year=1986:2006;
lon=-180:2.5:177.5;
lat=-90:2:90;

WIND=ncread('/Users/yangyang/Extreme_Aerosol/Meteorology/WIND850SPEED/ncfile/WIND850SPEED_daily.nc','wind850speed');
PM25=ncread('/Users/yangyang/Extreme_Aerosol/MET/ncfile/Aerosol_daily.nc','pm25');

WIND_DJF=zeros(length(lon),length(lat),length(year));
WIND_DJF(:,:,:)=mean(WIND,3);
WIND_DJF(:,:,13)=(WIND_DJF(:,:,11)+WIND_DJF(:,:,12))/2;

PM25_DJF=zeros(length(lon),length(lat),length(year));
PM25_DJF(:,:,:)=mean(PM25,3);
PM25_DJF(:,:,13)=(PM25_DJF(:,:,11)+PM25_DJF(:,:,12))/2;

WIND_ECHINA=zeros(21,1);
WIND_ECHINA(:)=mean(mean(WIND_DJF(117:122,57:68,:),2),1);

y=zeros(length(year),1);
x=zeros(length(year),1);

reg=zeros(length(lon),length(lat));
sig=zeros(length(lon),length(lat));
for i=1:length(lon)
    for j=1:length(lat)
        y(:)=PM25_DJF(i,j,:);
        x(:)=WIND_ECHINA;
        X=[ones(length(year),1),x];
        [b,bint,r,rint,stats]=regress(y,X);
        reg(i,j)=b(2);
        sig(i,j)=double(stats(3)<0.05);        
    end
end


reg(isnan(reg))=-999;


save('reg_PM25_WIND_EC.mat','reg')

nccreate('reg_PM25_WIND_EC.nc','lat','Dimensions',{'lat',length(lat)},'Format','classic')
ncwrite('reg_PM25_WIND_EC.nc','lat', lat)
ncwriteatt('reg_PM25_WIND_EC.nc','lat','units','degrees_north')
ncwriteatt('reg_PM25_WIND_EC.nc','lat','long_name','Latitude')
nccreate('reg_PM25_WIND_EC.nc','lon','Dimensions',{'lon',length(lon)},'Format','classic')
ncwrite('reg_PM25_WIND_EC.nc','lon', lon)
ncwriteatt('reg_PM25_WIND_EC.nc','lon','units','degrees_east')
ncwriteatt('reg_PM25_WIND_EC.nc','lon','long_name','Longitude')
nccreate('reg_PM25_WIND_EC.nc','reg','Dimensions',{'lon',length(lon),'lat',length(lat)},'Format','classic')
ncwrite('reg_PM25_WIND_EC.nc','reg', reg)
nccreate('reg_PM25_WIND_EC.nc','sig','Dimensions',{'lon',length(lon),'lat',length(lat)},'Format','classic')
ncwrite('reg_PM25_WIND_EC.nc','sig', sig)
ncdisp('reg_PM25_WIND_EC.nc')
