read_prevah_input(file,nodata)

function [R, header] = read_prevah_input(file,nodata)
% PUPROSE: read input tabular data for PREVAH hydrological model
% -------------------------------------------------------------------
% USAGE: [data, header] = read_prevah_input(file)
% Input:  file - filename and location
%         nodata - nodata value (optional, if given value will be replaced
%         with NaN)
% Output: data   - 4*n column vector: year - month - day - runoff 
%                  where n is the number of timesteps
%         header - 4*1 cell: elevation, x and y coordinate, name
%
% EXAMPLE: [data, header] = read_prevah_input('rain2003.tab')
% 
% See also: write_prevah_grid, read_prevat_grid, read_prevah_result
%
% Felix Hebeler, Geography Dept., University Zurich, Feb 2009

% 
data = nan(1,4);

%open file
fid = fopen(file,'r');

% read header
% read first line and ignore
title = fgetl(fid);
% read second line
remain = fgetl(fid);
col=1;
while true  % loop once over elevation to find number of columns 
   [str, remain] = strtok(remain);
    if isempty(str),  break;  end
   header{1,col}=str;
   col=col+1;
end
col=col-1;
remain = fgetl(fid);
for i=1:col %loop again for xcoords in 3rd row
   [str, remain] = strtok(remain);
    if isempty(str),  break;  end
   header{2,i}=str;
end
remain = fgetl(fid);
for i=1:col %loop 3rd time for ycoords in 3rd row
   [str, remain] = strtok(remain);
    if isempty(str),  break;  end
   header{3,i}=str;
end
remain = fgetl(fid);
for i=1:col %loop 4th time for name in 4th row
   [str, remain] = strtok(remain);
    if isempty(str),  break;  end
   header{4,i}=str;
end

% construct format string for data
format = '%d %d %d %d'; % first 4 columns are read in as integers (YY MM DD HH)
for i=5:col
    format = [format,' %n']; % add a %n to read in double for each data column
end
data = textscan(fid,format);

%intialise result array
l=length(cell2mat(data(1,col))); % get length of last column, this is used because result files can contain gibberish from some point on,
% which will not be read, but the column length (no of timesteps) will not
% be correct anymore for the reformating
R=nan(l,col);

% reformat to column vector
for i= 1:col
    R(:,i)=data{1,i}(1:l); % l (normally, l equals end)
end
clear data;
% check if nodata values should be replaced
if exist('nodata','var')
    R(R==nodata)=NaN;
end
%close file
fclose(fid);