source: trunk/src/read_lai_gis.m @ 327

function [lon_value, lat_value, datejul_value, datestr_value, lai_value, gridcode_value, classes_value,msd_value, msdnom_value]=read_lai_gis(yyyy,nbligne)

%READ_LAI_GIS fill LAI 2D array and associated arrays date, lat, lon, vegetation classe, pedologiy classe, etc. with data in :file:`${PROJECT_ID}/LAI/LAI{yyyy}F_veg2010_sol.csv`

%
%+
%
% ==============
% read_lai_gis.m
% ==============
%
% .. function:: read_lai_gis(yyyy)
%
% DESCRIPTION
% ===========
%
%    :param yyyy: year
%    :type yyyy: int16
%                ++todo++ je ne comprend pas: avec un int16, num2str ne fonctionne pas sous octave
%    :type yyyy: int8 (ordinaire)
%    :raise yyyy: required ex 2008
%
%    :param nbligne : nb ligne a lire ex 10000
%    :type nbligne : int
%    :raise nbligne: required
%
% read a file :file:`${PROJECT_ID}/LAI_CSV/LAI{yyyy}F_veg2010_sol.csv`
%
% return associated arrays lat, lon, lai, vegetation class , pedology class
%
% Les dates au format `yyyymmddhhmn` sont stockées dans le tableau
% **datestr_value** et les dates en jour julien sont stockées dans
% **datejul_value**.
% test realise sur des vraies donnees
% corrige pb de colonne mapusub souvent vide de donnees avec
% textedit (remplace ;; par ;" ";)
%lecture reussie sur un nb limite de lignes, mais pb memoire pour
%traitement du fichier complet (sur mon Mac)
%corrige quelques pbs dans le prgm (inversion colonnes, ismember
%remplace par isnan,...)
%
% pb de memoire: impossible d'avoir l'allocation memoire pour les 50000 lignes du fichier
% fonctionne jusqu'a 45000 lignes (sur cratos).
%pour developpements d'analyse LAI par Soukeye, ajout du parametre
%nbligne dans la fonction 
% lelod 20110728
%
% EXAMPLES
% ========
%
%
% TIPS
% ====
%
% To see the first line (header) of this kind of file::
%
%   $  sed -ne "1,1p" lai2008_Clip_Intersect_Inter.csv
%
%
% To see the second line (first line of data) of this kind of file::
%
%   $  sed -ne "2,2p" lai2008_Clip_Intersect_Inter.csv
%
%
% SEE ALSO
% ========
%
% :ref:`data_lai`
% :ref:`zones`
%
% :func:`simul_lai_gis` ++
% :func:`write_lai_gis`
% :func:`plot_lai_gis` ++
%
% :func:`read_lai_2d.m`
%
% TODO
% ====
%
% add zi_code parameter when filename convention will be changed
%
%
% il ne devrait plus y avoir des lignes avec tous le LAI à NaN ici
%
%
% hard coded value of nb_line_data_max; 
% handling eof vs nb_line_data_max
%
% coding rules
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20110708T085037Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * reading data line using technique of read_zi.m
%     because of string values dmlread connaot be used
%   * file written in correct format (. for decimal, ";" as
%     separator, float values for numbers (no more string for
%     numeric values)   - lelod 20110720
%   * add datejul_value output and all GIS info 
%
% - fplod 20110706T154713Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * creation draft from la2008i_Clip_Intersect_Inter.csv produce by
%     openoffice saving la2008i_Clip_Intersect_Inter.dbf in CSV format
%     and read_lai_2d.m
%
%     cautions : some modification on the original test file
%     because too many user no friendly features.
%
%-
%
global PROJECT_ID
%
%usage='[lon_value, lat_value, datejul_value, datestr_value,lai_value, gridcode_value, classes_value,msd_value, msdnom_value]=read_lai_gis(yyyy, nbligne)';
%
%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 filename to be read++
%fullfilename=['/Users/lelod/progX11/varamma_id/senegal/LAI/LAI_CSV/LAI' num2str(yyyy) 'F_veg2010_sol.csv'];
%
%fullfilename=['/usr/temp/lelod/varamma_d/LAI_CSV/LAI' num2str(yyyy) 'F_veg2010_sol.csv'];
%
fullfilename=[ PROJECT_ID 'LAI_CSV/LAI' num2str(yyyy) 'F_veg2010_sol.csv'];
disp(['iii : opening for reading ', fullfilename]);
% opening file
[fid,message]=fopen(fullfilename,'r');
% test if fullfilename exists
if (fid == -1)
    disp(message)
    error([ fullfilename ' not exist'])
end
clear message
%
% first line contains "lat","lon" and
% an unknown quantity of days of years
% so we read this line and loop to find the number of columns
%
% get first 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);
% attention les semaines sont de 7 jours )--> 52 semaines (avec 8j,
% on tombe a 46)
doy=zeros(nb_doy_max,1);
doy(:)=NaN;
while ( ~isempty(header_line) )
  % parse next column label
  [next,header_line] = strtok(header_line,';');
  switch next
     case {'"Y"'}
     case {'"X"'}
    case {'"FID_occups"'}
     case {'"ID"'}
     case {'"GRIDCOD"'}
     case {'"FID_bvferl"'}
    case {'"ID_1"'}
     case {'"GRIDCODE_1"'}
     case {'"classe"'}
     case {'"AREA"'}
     case {'"PERIMETER"'}
     case {'"SGSOLRU_"'}
    case {'"SGSOLRU_ID"'}
     case {'"MAPU"'}
     case {'"MAPUSUB"'}
     case {'"MSD"'}
     case {'"MSDNOM"'}
     otherwise
        % here we expect somesthing like "doy"
        % for example day 9 will be coded "9"
        % ++ todo regexp
        [sdoy, extra_sdoy] = strtok(next,';');
        clear extra_sdoy 
        doy_string = sdoy;
        %(1,3:end); plus utile, puisque header nettoye a la fabrication
        clear sdoy
        doy(index_column - 1) = str2double(doy_string);
        clear doy_string
  end
  clear next
  % increase index of column
  index_column = index_column + 1;
end
%
% we know now how many columns are written on header
nb_column = index_column;
nb_doy = nb_column - 2 - 15
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
    datejul_value(index_doy)=datenum(yyyy_d,1,1) + doy(index_doy) - 1;
end
datestr_value=datestr(datejul_value,'yyyymmdd');
%
% preallocation of output arrays to missing value: attention:
% entre en parametre
nb_line_data_max = nbligne;
%45000 valeur max;
%warning!!! 20 lignes pour le test - a ajuster sur le plus gros
%fichier (50247 lignes?)
lat_value=zeros(1,nb_line_data_max);
lat_value(:)=NaN;
lon_value=zeros(1,nb_line_data_max);
lon_value(:)=NaN;
lai_value=zeros(nb_doy,nb_line_data_max);
lai_value(:,:)=NaN;
gridcode_value=zeros(1,nb_line_data_max);
gridcode_value(:)=NaN;
%gridcode_1_value=zeros(nb_line_data_max);
%gridcode_1_value(:)=NaN;
classes_value=repmat({'NaN'},nb_line_data_max,1);
%area_value=zeros(nb_line_data_max);
%area_value(:)=NaN;
%perimeter_value=zeros(nb_line_data_max);
%perimeter_value(:)=NaN;
%sgsolru_value=zeros(nb_line_data_max);
%sgsolru_value(:)=NaN;
%mapu_value=repmat({'NaN'},nb_line_data_max,1);
%mapusub_value=repmat({'NaN'},nb_line_data_max,1);
msd_value=zeros(1,nb_line_data_max);
msd_value(:)=NaN;
msdnom_value=repmat({'NaN'},nb_line_data_max,1);
%
% reading data from file (from second line till end of file)
%lat_value=dlmread(fid,';',0,1)

index_line_data=1;
%
%Attention: fichier trop gros. On ne lit que les premieres nbligne lignes
%while (feof(fid)~= 1)
while (index_line_data <nbligne+1)
%
%
  data_line = fgetl(fid);
  % loop to handle eahc item on the data line (";" separated value)
  % attention, il arrive qu'une cellule soit vide (mapusub) : corrige en principe
  index_item=1;
  while ( ~isempty(data_line))
     % parse next column label
     [item,data_line] = strtok(data_line,';');
     if (index_item == 1 )
       lat_value(index_line_data)=str2double(item);
     elseif (index_item == 2 )
       lon_value(index_line_data)=str2double(item);
     elseif (index_item > 2 && index_item <= 2 + nb_doy)
       lai_value(index_item - 2, index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 1)
   %    fid_occups_value(index_line_data)=str2double(item);
    elseif (index_item ==  2 + nb_doy + 2)
    %   id_value(index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 3)
       gridcode_value(index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 4)
    %   fid_bvferl_value(index_line_data)=str2double(item);
    elseif (index_item ==  2 + nb_doy + 5)
   %   id_1_value(index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 6)
    %   gridcode_1_value(index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 7)
       classes_value(index_line_data)=cellstr(item);
     elseif (index_item ==  2 + nb_doy + 8)
    %   area_value(index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 9)
    %   perimeter_value(index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 10)
    %   sgsolru_value(index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 11)
    %   sgsolru_id_value(index_line_data)=str2double(item);
     elseif (index_item ==  2 + nb_doy + 12)
    %   mapu_value(index_line_data)=cellstr(item);
     elseif (index_item ==  2 + nb_doy + 13)
    %   mapusub_value(index_line_data)=cellstr(item);
    elseif (index_item ==  2 + nb_doy + 14)
       msd_value(index_line_data)=str2double(item);
    elseif (index_item ==  2 + nb_doy + 15)
       msdnom_value(index_line_data)=cellstr(item);
    else
       disp([' eee : data_line ', data_line ]);
       disp([' eee : item ', item]);
       disp([' eee : index_line_data ', num2str(index_line_data) ]);
       disp([' eee : index_item ', num2str(index_item) ]);
       error([ 'invalid data line (too many items) ']);
    end
    clear item
    % increase index of item
    index_item = index_item + 1;
  end
  % check for number of item in current data line
  if (index_item < nb_doy + 2 + 15)
     disp([' eee : data_line ', data_line ]);
     disp([' eee : index_line_data ', num2str(index_line_data) ]);
     disp([' eee : index_item ', num2str(index_item) ]);
     error([ 'invalid data line (missing items) :'])
  end
  index_line_data=index_line_data+1;
  clear data_line
  clear index_item
end
if (index_line_data > nb_line_data_max+1)
  error([' need to increase nb_line_data_max up to ' num2str(index_line_data) ]);
end
%
clear index_line_data
%
fclose(fid);
clear fid
%
% remove extra data NaN values due to preallocation
% criteria = when lon_value is not NaN
missing_data_indices=find(isnan(lon_value));
if (~isempty(missing_data_indices))
    nb_data=min(missing_data_indices)-1
    lat_value=lat_value(1:nb_data);
    lon_value=lon_value(1:nb_data);
    lai_value=lai_value(:,1:nb_data);
    %fid_occups_value=fid_occups_value(1:nb_data);
    %id_value=id_value(1:nb_data);
    gridcode_value=gridcode_value(1:nb_data);
    %fid_bvferl_value=fid_bvferl_value(1:nb_data);
    %id_1_value=id_1_value(1:nb_data);
    %gridcode_1_value=gridcode_1_value(1:nb_data);
    classes_value=classes_value(1:nb_data);
    %area_value=area_value(1:nb_data);
   % perimeter_value=perimeter_value(1:nb_data);
    %sgsolru_value=sgsolru_value(1:nb_data);
    %sgsolru_id_value=sgsolru_id_value(1:nb_data);
    %mapu_value=mapu_value(1:nb_data);
   % mapusub_value=mapusub_value(1:nb_data);
    msd_value=msd_value(1:nb_data);
    msdnom_value=msdnom_value(1:nb_data);
end
%
clear nb_line_data_max
%
%whos
%
result=0;
end

%!demo
%! % To read :file:`${PROJECT_ID}/LAI/la2008i_Clip_Intersect_Inter.csv`:
%! clear all
%! close all
%! varamma_startup
%! more off
%! yyyy=int16(2008);
%! [lon_value, lat_value, datejul_value, datestr_value, lai_value, fid_occups_value, id_value, gridcode_value, fid_bvferl_value, id_1_value, gridcode_1_value, classes_value, area_value, perimeter_value, sgsolru_value, sgsolru_id_value, mapu_value, mapusub, msd_value, msdnom_value]=read_lai_gis(yyyy);