clear all;
close all;

% read file
fid=fopen('v2.mean_adj.txt');
A=textscan(fid,'%3d%5d%3d%1d%4d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d');
fclose(fid);

% convert to matrix
A=cell2mat(A);
A=double(A);

% sort
A(A==-9999)=nan;
A=sortrows(A,[2 5]);

% gather data from stations
stations=unique(A(:,2));
numstations=length(unique(A(:,2)));
data=NaN(numstations, 12*(max(A(:,5))-min(A(:,5))));

Arows=size(A);
yearmin=min(A(:,5));

disp('Parsing temperature data into chronological order...');

for nn=1:Arows
    % count station number
    for ii=1:numstations
        if A(nn,2)==stations(ii)
            break;
        end
    end
    disp(['Processing station ' num2str(stations(ii))...
            '/' num2str(stations(numstations)) ', year '...
                num2str(A(nn,5)) '...' ]);
    % append temperature data
    year=A(nn,5);
    for jj=1:12
        data(ii,12*(year-yearmin)+jj)=A(nn,jj+5)./10;
    end
end

% number of years
n=(max(A(:,5))-min(A(:,5)));

tempmean=zeros(1,n);
tempsem=zeros(1,n);
tempstd=zeros(1,n);

disp('Computing yearly temperature statistics...');

for nn=1:n
    tempmean(nn)=nanmean(nanmean(data(:,12*(nn-1)+[1:12]),2));    
    tempstd(nn)=nanstd(nanmean(data(:,12*(nn-1)+[1:12]),2));        
    tempsem(nn)=tempstd(nn)./(sqrt(numstations));
end

% plot standard deviation
hold on;
yearlow=1850-yearmin;
yearhi=2010-yearmin-1;
alpha=1-0.95;
z=norminv(1-alpha/2,0,1);
val=[yearlow:8:yearhi];
for ii=1:length(val)
    plot([val(ii) val(ii)]+yearmin,[tempmean(val(ii))-tempstd(val(ii))...
                            tempmean(val(ii))+tempstd(val(ii))],'k-');
end
% plot mean with 95% confidence intervals
val=[yearlow:1:yearhi];
errorbar(val+yearmin, tempmean(val), z*tempsem(val));
set(gca,'xlim',[1850 2010],'ylim',[3 20]);
title('global mean temperature (95% confidence intervals, one sigma standard deviation)');
ylabel('temperature (C degrees)');
xlabel('year');
box on;
hold off;

print(gcf,'-dpdf','WORLD_MEAN.pdf');
print(gcf,'-djpeg90','WORLD_MEAN.jpg');