%%%%% READING GHCN DATA
%  1. Read the data using 'columnar table'
%  2. 

clear;
close('all');
loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCDC/.GHCN/.v2beta/lon/%28-20%29%2830%29masknotrange/SELECT/lat/%2812.5%29%2820%29masknotrange/SELECT/a:/.lon/:a:/.lat/:a:/.prcp/:a/T/%28Jul%201900%29%28Sep%202008%29RANGEEDGES/T/3/1./boxAverage/T/12/STEP/-999./setmissing_value/dods');
data=prcp.prcp;
stat=prcp.IWMO;
[ns,ny]=size(data);
loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCDC/.GHCN/.v2beta/lon/%28-18%29%2830%29masknotrange/SELECT/lat/%2812.5%29%2820%29masknotrange/SELECT/.lon/dods');
lon=lon.lon;
loaddap('http://iridl.ldeo.columbia.edu/expert/SOURCES/.NOAA/.NCDC/.GHCN/.v2beta/lon/%28-18%29%2830%29masknotrange/SELECT/lat/%2812.5%29%2820%29masknotrange/SELECT/.lat/dods');
lat=lat.lat;

% number of existing station per year and the running-mean of the count
for i=1:ny
    c=find(data(:,i)==-999);
    count(i)=ns-length(c);
end
% running mean to find the period with the highest # of available stations
for i=16:ny-15
aa(i-15)=mean(count(i-15:i+15));
end
% in case of tie choose the most recent period
a=find(aa==max(aa));
d=max(a);
c=d(end);
cr=c+15; % actual order
cy=cr-1+1900; % in year unit
% any station having less than 25 out of 31 years of data in the normal period
i=0;
for is=1:ns
    c=find(data(is,cr-15:cr+15)==-999);
    cc=length(c);
    if (cc<6)  % 25 out 0f 31 = 6 missing only
        i=i+1;
        ccount(i)=is;
    end
end
figure(1)
plot(count,'-k','LineWidth',2)
%'MarkerSize',3)
hold on
plot([16:ny-15],aa,'-r','LineWidth',2)
plot([cr-15:cr+15],[100:100],'-ob','LineWidth',3)
xlim([1 109]);
set(gca,'XTick',1:10:109)
set(gca,'XTickLabel',1900:10:2008)
title('GHCN available stations over 12.5N-20N and 20W-30E')
grid on
eval(sprintf('print SahStat0008GHNC.tiff -dtiff'))


% take only those stations
data1=data(ccount,:);
lon1=lon(ccount);
lat1=lat(ccount);
stat1=stat(ccount);
[ns,ny]=size(data1);
%from here just go to excel lazy ...

% writting in a file to compute the anomalies
fid=fopen('GHCN008JAS12p5N20N.txt','wt');
fprintf(fid,'Order ');
for i=1:ns fprintf(fid,'%d ',i); end
fprintf(fid,'\n');

fprintf(fid,'Stations ');
for i=1:ns fprintf(fid,'%d ',stat1(i)); end
fprintf(fid,'\n');

fprintf(fid,'Long ');
for i=1:ns fprintf(fid,'%.2f ',lon1(i)); end
fprintf(fid,'\n');
fprintf(fid,'Lat ');
for i=1:ns fprintf(fid,'%.2f ',lat1(i)); end
fprintf(fid,'\n');

for iy=1:ny
fprintf(fid,'%d ',iy+1899);
	for is=1:ns 	fprintf(fid,'%.2f ',data1(is,iy)); 	end
	fprintf(fid,'%d \n',count(iy));
end
fclose(fid);




