source: trunk/src/read_lai_2d.m @ 327

function [lon_value, lat_value, datestr_value, lai_value]=read_lai_2d(yyyy)

%READ_LAI_2D fill array LAI 2D array and associated arrays date, lat, lon
% with data in :file:`${PROJECT_ID}/LAI/laisen{yyyy}_float.txt`

%
%+
%
% =============
% read_lai_2d.m
% =============
%
% .. function:: read_lai_2d(yyyy)
%
% DESCRIPTION
% ===========
%
%    :param yyyy: year
%    :type yyyy: int16
%    :raise yyyy: required
%
% read a file :file:`${PROJECT_ID}/LAI/laisen{yyyy}_float.txt`
%
% Les LAI sont stockées dans le tableau 2D (temps,ligne) **lai_value**
%
% Les longitudes sont stockées dans le tableau 1D **lon_value**.
%
% Les latitudes sont stockées dans le tableau 1D **lat_value**.
%
% Les dates au format `yyyymmdd` sont stockées dans le tableau 1D **datestr_value**.
%
% EXAMPLES
% ========
%
% To read :file:`${PROJECT_ID}/LAI/laisen2001_float.txt`:
%
% see demo 1
%
% A file :file:`${PROJECT_ID}/LAI/laisen2001_write_lai_2d_float.txt` can be written
% with data we just read::
%
%  octave> result=write_lai_2d(lon_value, lat_value, datestr_value, lai_value);
%
% Files comparison might help detecting problems (either on files or in programs involved)::
%
%  $ diff ${PROJECT_ID}/LAI/laisen2001_float.txt ${PROJECT_ID}/LAI/laisen2001_write_lai_2d_float.txt
%
% SEE ALSO
% ========
%
% :ref:`guide data LAI <data_lai>`
%
% :func:`simul_lai_2d`
% :func:`write_lai_2d`
% :func:`showgrid`
%
% TODO
% ====
%
% status de retour
%
% description du résultat en ReST
%
% hard coded value of nb_data_line_max;
% should be calculated according to latitude and longitude range and resolution
%
% EVOLUTIONS
% ==========
%
% $Id: read_lai_2d.m 312 2011-07-07 17:13:54Z pinsard $
%
% $URL$
%
% - fplod 20110228T142213Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * move example to demo
%   * use int2str instead of num2str
%
% - fplod 20101208T102711Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * explicit yyyy type int16 (side effect : modification of datenum call)
%   * check argument
%
% - fplod 20101206T152338Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * check argument
%
% - fplod 20101206T103915Z adonis.locean-ipsl.upmc.fr (Linux)
%
%   * cleanning thanks to mlint
%
% - fplod 20101203T171418Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * comma also on the first line (for qgis)
%
% - fplod 20101201T153337Z zeus.locean-ipsl.upmc.fr (Linux)
%
%   * new terminology and new format
%   * add parameter information in header
%
% - fplod 20101022T140456Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add test on octave version because of dlmread problem running
%     octave 3.0.2
%
% - fplod 20101022T104521Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * bug fix : replace '' by ' ' in dlmread instruction
%
% - fplod 20100913T100256Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * creation from read_lai.m because impossible to define
%     lai_values(46,839041,839041) - too big
%   * faster using dlmread
%
%-
%
global PROJECT_ID
%
usage='lon_value, lat_value, datestr_value, lai_value]=read_lai_2d(yyyy)';
%
if nargin~=1
   disp(['Incorrect number of arguments']);
   error(usage);
end
%
arg_info=whos('yyyy');
if ~strcmp(arg_info.class,'int16')
   disp(['Incorrect type of arg yyyy']);
   whos yyyy
   error(usage);
end
clear arg_info
%
% build LAI ASCII filename to be read
fullfilename=[PROJECT_ID 'LAI/laisen' int2str(yyyy) '_float.txt'];
disp(['iii : opening for reading ', fullfilename]);
% opening file
fid=fopen(fullfilename,'r');
% test if fullfilename exists
if (fid == -1)
    disp(message)
    error([ fullfilename ' not exist'])
end
clear fullfilename
%
% read file
%
% first line contains "Y,X" and an unknown quantity of days of years
% so we read this line and loop to find the number of columns
%
% get next line as a string
header_line = fgetl(fid);
%
% loop to find the number of columns
% and initialize the "day of year" array
index_column = 0;
% the maximum of doy is the number of days in year yyyy
% divided by 8 (according to time resolution of this dataset)
yyyy_d=double(yyyy);
nb_doy_max=ceil((datenum(yyyy_d,12,eomday(yyyy_d,12)) - datenum(yyyy_d,1,1))/8);
doy=zeros(nb_doy_max,1);
doy(:)=NaN;
datestr_value=doy;
while ( ~isempty(header_line) )
  % parse next column label
  [next,header_line] = strtok(header_line,',');
  switch next
     case {'X'}
     case {'Y'}
     otherwise
        doy(index_column - 1) = str2double(next);
  end
  clear next
  % increase index of column
  index_column = index_column + 1;
end
% close file
fclose(fid);
%
% we know now how many columns are written on header
nb_column = index_column;
nb_doy = nb_column - 2;
clear index_column
%
% remove extra doy NaN values due to preallocation
doy=doy(~isnan(doy));
%
% initialisation d'un tableau string de nb_doy elements à 9999999
datestr_value=repmat('99999999',nb_doy,1);
for index_doy=1:nb_doy
    datestr_value(index_doy,:)=datestr((datenum(yyyy_d,1,1) + doy(index_doy) - 1),'yyyymmdd');
end
%
% following lines contains lat, lon and LAI values
% we first preallocate arrays (boite,dx,dy formula or
% equivalent of IDL function file_lines ++)
% and initialize values to NaN
nb_data_line_max = 231002 - 2 ; %++
lat_value=zeros(1,nb_data_line_max);
lat_value(:)=NaN;
lon_value=zeros(1,nb_data_line_max);
lon_value(:)=NaN;
lai_value=zeros(nb_doy,nb_data_line_max);
lai_value(:,:)=NaN;
%
% reading data from file (from second line till end of file)
data = dlmread(fullfilename,',', 1, 0);
%
% put values into arrays to be used for plot or computation
lat_value=data(:,1);
lon_value=data(:,2);
lai_value=data(:,3:size(data,2));
clear data
%
end

%!demo
%! % To read :file:`${PROJECT_ID}/LAI/laisen2001_float.txt`::
%! varamma_startup
%! more off
%! yyyy=int16(2001);
%! [lon_value, lat_value, datestr_value, lai_value]=read_lai_2d(yyyy);