year=1981:2015;

files=dir('/Users/yangyang/Extreme_Aerosol/VIS/dat/1980_2014/*.txt');

lon=zeros(length(files),1);
lat=zeros(length(files),1);

vis_mean=zeros(length(files),length(year));
vis_mean_average=zeros(length(year),1);
vis_mean_trend=zeros(length(files),1);

vis_ndays_temp=zeros(length(files),length(year));
vis_ndays=zeros(length(files),length(year));

vis_ndays_average=zeros(length(year),1);
vis_ndays_trend=zeros(length(files),1);


for FileID=1:length(files)
    
    temp=importdata(['/Users/yangyang/Extreme_Aerosol/VIS/dat/1980_2014/',files(FileID).name],' ');
    lat(FileID)=temp(2);
    lon(FileID)=temp(3);
    temp=importdata(['/Users/yangyang/Extreme_Aerosol/VIS/dat/1980_2014/',files(FileID).name],' ',2);
    ncount=length(temp.data);
    month=temp.data(:,3);
    vis=temp.data(:,5);
    prcp=temp.data(:,6);
    rh=temp.data(:,7);
    t=temp.data(:,8);
    
    YearID=1;
    month(ncount+1)=12;
    days=0;
    
    for DayID=1:ncount
        
        days=days+1;
        vis_mean(FileID,YearID)=vis_mean(FileID,YearID)+vis(DayID);
        
        if vis(DayID)<10
            vis_ndays_temp(FileID,YearID)=vis_ndays_temp(FileID,YearID)+1;
        end
        
        if month(DayID)==2 && month(DayID+1)==12
            vis_mean(FileID,YearID)=vis_mean(FileID,YearID)/days;
%            vis_ndays(FileID,YearID)=vis_ndays_temp(FileID,YearID);
            vis_ndays(FileID,YearID)=vis_ndays_temp(FileID,YearID)/days*90;
            days=0;
            YearID=YearID+1;
        end
    end
    
end


temp=zeros(length(year),1);
temp2=zeros(length(year),1);
n=0;
for FileID=1:length(files)
    temp(:)=vis_ndays(FileID,:);
    temp2(:)=vis_mean(FileID,:);
    if lon(FileID)>=110 && lon(FileID)<=122.5 && lat(FileID)>=30 && lat(FileID)<=45
        n=n+1;
        vis_ndays_average=vis_ndays_average+temp;
        vis_mean_average=vis_mean_average+temp2;
    end
end
vis_ndays_average=vis_ndays_average/n;
vis_mean_average=vis_mean_average/n;


y=zeros(length(year),1);
for FileID=1:length(files)
    y(:)=vis_ndays(FileID,:);
    X=[ones(length(year),1),year'];
    [b,bint,r,rint,stats]=regress(y,X);
    vis_ndays_trend(FileID)=b(2)*10;
    if stats(3)>0.05
       vis_ndays_trend(FileID)=-999;
    end
end

y=zeros(length(year),1);
for FileID=1:length(files)
    y(:)=vis_mean(FileID,:);
    X=[ones(length(year),1),year'];
    [b,bint,r,rint,stats]=regress(y,X);
    vis_mean_trend(FileID)=b(2)*10;
    if stats(3)>0.05
       vis_mean_trend(FileID)=-999;
    end
end

save('VIS_1981_2015_NC.mat','lat','lon','vis_mean','vis_mean_trend','vis_ndays','vis_ndays_trend')

nccreate('VIS_1981_2015_NC.nc','Year','Dimensions',{'year',length(year)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','Year', year)
ncwriteatt('VIS_1981_2015_NC.nc','Year','long_name','Year 1973-2014')
nccreate('VIS_1981_2015_NC.nc','lat','Dimensions',{'nstation',length(files)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','lat', lat)
ncwriteatt('VIS_1981_2015_NC.nc','lat','units','degrees_north')
ncwriteatt('VIS_1981_2015_NC.nc','lat','long_name','Latitude')
nccreate('VIS_1981_2015_NC.nc','lon','Dimensions',{'nstation',length(files)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','lon', lon)
ncwriteatt('VIS_1981_2015_NC.nc','lon','units','degrees_east')
ncwriteatt('VIS_1981_2015_NC.nc','lon','long_name','Longitude')
nccreate('VIS_1981_2015_NC.nc','vis_mean','Dimensions',{'nstation',length(files),'year',length(year)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','vis_mean', vis_mean)
nccreate('VIS_1981_2015_NC.nc','vis_mean_average','Dimensions',{'year',length(year)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','vis_mean_average', vis_mean_average)
nccreate('VIS_1981_2015_NC.nc','vis_mean_trend','Dimensions',{'nstation',length(files)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','vis_mean_trend', vis_mean_trend)
nccreate('VIS_1981_2015_NC.nc','vis_ndays','Dimensions',{'nstation',length(files),'year',length(year)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','vis_ndays', vis_ndays)
nccreate('VIS_1981_2015_NC.nc','vis_ndays_average','Dimensions',{'year',length(year)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','vis_ndays_average', vis_ndays_average)
nccreate('VIS_1981_2015_NC.nc','vis_ndays_trend','Dimensions',{'nstation',length(files)},'Format','classic')
ncwrite('VIS_1981_2015_NC.nc','vis_ndays_trend', vis_ndays_trend)
ncdisp('VIS_1981_2015_NC.nc')
