% data filling of Colombia IDEAM data using stochastic daily weather
% generator

% March 2009
% May 2009: fill missings in updated dataset (station names have changed,
% and there are now six stations)

% required functions:
% stochastic_filling.m:
% -- here http://portal.iri.columbia.edu/portal/server.pt?open=512&objID=697&PageID=7575&cached=true&mode=2&userID=2

% fill_missing_temperature_days.m
% -- here:
% http://iri.columbia.edu/~awr/TP/Colombia/May2009/fill_missing_temperature_days.m

clear
cd /Users/andy/data/Americas/Colombia/May09
% read data from IRI Data Library into Matlab (here via NetCdf files)
% precip:
% http://iridl.ldeo.columbia.edu/home/.ciph/.Colombia/.pilot_sites/.climate
% /.station/.daily/.PRCP/
% ISTA
%     grid: /ISTA (ids) ordered [ (2501001) (2502516) (3101501) (3210507) (5304503) (5311501)] :grid 
% Time
%     grid: /T (julian_day) ordered (1 Jan 1942) to (28 Feb 2009) by 1. N= 24531 pts :grid 
% % (67 years)
ncload('prcp_IDEAM2.nc')
ncload('Tmin_IDEAM2.nc')
ncload('Tmax_IDEAM2.nc')
ncload('Tmean_IDEAM2.nc')

nsta=6; % number of stations
nyr=67; % number of years: 1942-08

% station names: NB, these have CHANGED (as per 5/19/09)
sname{1}='PTO LIBERTADOR';
sname{2}='CUBA HDA';
sname{3}='TRUENO EL';
sname{4}='SAN JOSE DEL GUAVIARE';
sname{5}='GUAPI';
sname{6}='APTO BUENAVENTURA';

% Files extend to Feb 09. Select only 1942-08 complete years.
ndat=24472;
PRCP=PRCP(:,1:ndat);
TMIN=TMIN(:,1:ndat);
TMAX=TMAX(:,1:ndat);
TMEAN=TMEAN(:,1:ndat);

% change missing data flag to NaN
PRCP(find((PRCP<-99)))=NaN;
TMIN(find((TMIN<-99)))=NaN;
TMAX(find((TMAX<-99)))=NaN;
TMEAN(find((TMEAN<-99)))=NaN;

% calendar dates (Don't know what "T" variable contains.)
dateoffset=datenum('1 Jan 1942');
tcal=datevec((1:24472)+dateoffset-1);

%% data screening
min(PRCP')
max(PRCP')
% ans =
%      0     0     0     0     0     0
% ans =
%   170.0000  145.0000  146.9000  140.0000  166.3000  274.0000

min(TMIN')
max(TMIN')
% ans =
%        NaN   18.0000   13.4000   12.2000   13.2000   16.9000
% ans =
%        NaN   26.5000   28.0000   25.2000   26.0000   26.0000

min(TMAX')
max(TMAX')
% ans =
%        NaN   24.7000   21.8000   22.0000   24.6000   24.0000
% ans =
%        NaN   39.6000   39.0000   38.7000   37.2000   35.0000

min(TMEAN')
max(TMEAN')
% ans =
%        NaN   22.1000   19.0000   20.9000   22.8000   22.7000
% ans =
%        NaN   32.5000   31.6000   30.6000   30.0000   29.3000

% further censoring: from Luis, 5/20/09:
% Precipitation (PRCP) of 274.0 mm (it is a doubtful data)..
% Minimal temperature (TMIN) below 15 celsius degree....
% Maximal temperature (TMAX) below 24 celsius degree...

for s=1:nsta
  knt(s)=length(find(PRCP(s,:)>250));
end

% knt = 0     0     0     0     0     1

for s=1:nsta
  knt(s)=length(find(TMIN(s,:)<15));
end
% knt = 0     0     4     7    17     0

for s=1:nsta
  knt(s)=length(find(TMAX(s,:)<24));
end
% knt = 0     0    15    44     0     0

PRCP(find((PRCP>250)))=NaN;
TMIN(find((TMIN<15)))=NaN;
TMAX(find((TMAX<24)))=NaN;

%% for plotting, remove leap days to make 365-day years
notFeb29=find(~(tcal(:,2)==2 & tcal(:,3)==29));
prcp2=PRCP(:,notFeb29); tmin2=TMIN(:,notFeb29); 
tmax2=TMAX(:,notFeb29); tmean2=TMEAN(:,notFeb29);
prcp2=reshape(prcp2,nsta,365,nyr);
tmin2=reshape(tmin2,nsta,365,nyr);
tmax2=reshape(tmax2,nsta,365,nyr);
tmean2=reshape(tmean2,nsta,365,nyr);

%% plot number of non-missing days per year at each station
% all years, precip
lseas=365;
yrs=1942:2008;
for s=1:nsta
  for yr=yrs
    fpres_prcp_all(s,yr-1941)=length(prcp2(isfinite(prcp2(s,:,yr-1941))))/lseas;
  end
end
figure(1); clf
  plot(yrs,fpres_prcp_all'); grid on
  legend(sname,'Location','NorthEastOutside')
  xlabel('Year'); ylabel('Fraction of days present')
  title('Station precipitation')
print -f1 -dpdf fractdata_precip_4208.pdf

%% plot number of non-missing days per year at each station: 2001-08
lseas=365;
yrs=2001:2008;
for s=1:nsta
  for yr=yrs
    fpres_prcp(s,yr-2000)=length(prcp2(isfinite(prcp2(s,:,yr-1941))))/lseas;
    fpres_tmin(s,yr-2000)=length(tmin2(isfinite(tmin2(s,:,yr-1941))))/lseas;
    fpres_tmax(s,yr-2000)=length(tmax2(isfinite(tmax2(s,:,yr-1941))))/lseas;
    fpres_tmean(s,yr-2000)=length(tmean2(isfinite(tmean2(s,:,yr-1941))))/lseas;
  end
end
figure(2); clf
subplot(2,2,1)
  plot(yrs,fpres_prcp'); grid on
  xlabel('Year'); ylabel('Fraction of days present')
  title('Station precipitation')
subplot(2,2,2)
  plot(yrs,fpres_tmin'); grid on
  xlabel('Year'); ylabel('Fraction of days present')
  title('Station Tmin')
subplot(2,2,3)
  plot(yrs,fpres_tmax'); grid on
  xlabel('Year'); ylabel('Fraction of days present')
  title('Station Tmax')
subplot(2,2,4)
  plot(yrs,fpres_tmean'); grid on
  legend(sname)
  xlabel('Year'); ylabel('Fraction of days present')
  title('Station Tmean')
print -f2 -dpdf fractdata_2001-08.pdf

%% fill missing PRCP data using a stochastic weather generator
% first insert missing data records for Feb 29s
nrec=24472;
i=0;
clear PRCP2
for t=1:nrec
  i=i+1;
  PRCP2(:,i)=PRCP(:,t);
  tcal2(i,:)=tcal(t,:);
  idum(i)=0;
  if(tcal(t,2)==2 && tcal(t,3)==28 && tcal(t+1,3)==1)
    i=i+1;
    PRCP2(:,i)=NaN;
    tcal2(i,:)=tcal(t,:); tcal2(i,3)=29;
    idum(i)=1; % tag
  end
end
% fill missings 
swg_months=[ones(31,1);ones(29,1)*2;ones(31,1)*3;ones(30,1)*4;...
  ones(31,1)*5;ones(30,1)*6;ones(31,1)*7;ones(31,1)*8;...
  ones(30,1)*9;ones(31,1)*10;ones(30,1)*11;ones(31,1)*12];

X2=zeros(366*nyr,nsta);
for j=1:nsta
  [X2(:,j)]=stochastic_filling(PRCP2(j,:),NaN,swg_months,0.0001,'gamma',1);
end

% remove Feb 29s in non-leap years
PRCPf=X2(find(idum==0),:)';

%% write out precip
fid = fopen('Precip_IDEAM2_daily_filled.txt', 'wt');
date=tcal(:,1)*10000 + tcal(:,2)*100 + tcal(:,3);
fmt='%8i'; for j=1:5; fmt=[fmt '%8.2f']; end; fmt=[fmt '\n'];
mat=[date PRCPf'];
%fprintf(fid,'%8i%8.2f%8.2f%8.2f%8.2f%8.2f%8.2f\n',mat');
fprintf(fid,fmt,mat');
fclose(fid)
clear mat fid

%% fill missing temperature data
TMINf = fill_missing_temperature_days(TMIN', PRCP', PRCPf', tcal);
TMAXf = fill_missing_temperature_days(TMAX', PRCP', PRCPf', tcal);
TMEANf = fill_missing_temperature_days(TMEAN', PRCP', PRCPf', tcal);
TMINf=TMINf'; TMAXf=TMAXf'; TMEANf=TMEANf';

%% write out
fid1 = fopen('Tmin_IDEAM2_daily_filled.txt', 'wt');
fid2 = fopen('Tmax_IDEAM2_daily_filled.txt', 'wt');
fid3 = fopen('Tmean_IDEAM2_daily_filled.txt', 'wt');
fmt='%8i'; for j=1:6; fmt=[fmt '%8.2f']; end; fmt=[fmt '\n'];
mat1=[date TMINf']; mat2=[date TMAXf']; mat3=[date TMEANf'];
fprintf(fid1,fmt,mat1'); fprintf(fid2,fmt,mat2'); fprintf(fid3,fmt,mat3');
fclose(fid1); fclose(fid2); fclose(fid3); 

%% check filled against unfilled data using epid. weekly averages
% first make the epid. weekly avgs
% organize weekly according to epidemiological calendar
% The first epidemiological week of the year ends on the first Saturday 
% of January, provided that it falls at least four or more days into the month.

[n, s]=weekday('1-Jan-1942')
% = thursday
% th f sa su m tu w
% 1  2 3  4  5
% so 1st epid weeks starts Jan 4, 1942

startday=4;
nweeks=24472/7 - 1; % one less week than length of time series
tcale=reshape(tcal(startday:startday-1+nweeks*7,:),7,nweeks,6);

temp=PRCP(:,startday:startday-1+nweeks*7);
PRCPe=reshape(temp,nsta,7,nweeks);
PRCPfe=reshape(PRCPf(:,startday:startday-1+nweeks*7),nsta,7,nweeks);

TMINe=reshape(TMIN(:,startday:startday-1+nweeks*7),nsta,7,nweeks);
TMINfe=reshape(TMINf(:,startday:startday-1+nweeks*7),nsta,7,nweeks);

TMAXe=reshape(TMAX(:,startday:startday-1+nweeks*7),nsta,7,nweeks);
TMAXfe=reshape(TMAXf(:,startday:startday-1+nweeks*7),nsta,7,nweeks);

TMEANe=reshape(TMEAN(:,startday:startday-1+nweeks*7),nsta,7,nweeks);
TMEANfe=reshape(TMEANf(:,startday:startday-1+nweeks*7),nsta,7,nweeks);

%% average weekly according to epid cal
PRCPe_m=squeeze(nanmean(PRCPe,2));
PRCPfe_m=squeeze(nanmean(PRCPfe,2));
TMINe_m=squeeze(nanmean(TMINe,2));
TMINfe_m=squeeze(nanmean(TMINfe,2));

TMAXe_m=squeeze(nanmean(TMAXe,2));
TMAXfe_m=squeeze(nanmean(TMAXfe,2));

TMEANe_m=squeeze(nanmean(TMEANe,2));
TMEANfe_m=squeeze(nanmean(TMEANfe,2));
%% check filled against unfilled: PRCP
%  precip_week: (31 Dec 2000 - 6 Jan 2001) to (23-29 Dec 2007)
%  PRCPe_m: 3079 - 3443

figure(1); clf
x=(1:365)/52+2001;
for s=1:6
  subplot(6,1,s)
  plot(x,PRCPfe_m(s,3079:3443)'*7,'b')
  hold on
  plot(x,PRCPe_m(s,3079:3443)'*7,'r')
  ylabel('Precip [mm/week]')
  axis([2001 2008 0 400])
  title(sname{s})
    if(s==1),legend('Filled','Unfilled'); end
end
print -f1 -dpdf check_precip2_filling.pdf

%% check filled against unfilled: T

figure(1); clf
x=(1:365)/52+2001;
for s=1:6
  subplot(6,1,s)
  plot(x,TMAXfe_m(s,3079:3443)','b')
  hold on
  plot(x,TMAXe_m(s,3079:3443)','r')
  ylabel('TMAX [degC]')
  title(sname{s})
    if(s==1),legend('Filled','Unfilled'); end
end
print -f1 -dpdf check_TMAX2_filling.pdf
