source: trunk/src/essailai.m @ 325

%+
% .. _essailai.m:
%
% ==========
% essailai.m
% ==========
%
% DESCRIPTION
% ===========
%
% read LAI file :file:`${PROJECT_ID}/LAI/laisen{yyyy}_float.txt`.
%
% plot maps of LAI for every dates on figures (4 days/page) if matlab is
% running (sorry for octave users : griddata is not ok either in 3.0.2 or
% 3.3.52)
%
% compute sum, std, moy
%
% plot sum, std, moy
% Figures are saved in :file:`${PROJECT_OD}/essailai_sum.ps`,
% :file:`${PROJECT_OD}/essailai_mean.ps` and
% :file:`${PROJECT_OD}/essailai_stddev.ps`.
%
% save sum, std, moy for each date in file :file:`${PROJECT_OD}/essailai.txt`
%
% EXAMPLES
% ========
%
% ::
%
%   $ octave
%   octave> varamma_startup
%   octave> more off
%   octave> essailai
%
%
% CAUTIONS
% ========
%
% If matlab is running, Statistical Toolbox is needed because of
% usage of :func:`nanmean`, :func:`nansum` and :func:`nanstd`.
%
% SEE ALSO
% ========
%
% :ref:`guide data LAI <data_lai>`
%
% :func:`varamma_startup`
%
% :func:`read_lai_2d`
% :func:`showgrid`
%
% :func:`nanmean`
% :func:`nanstd`
% :func:`nansum`
%
% TODO
% ====
%
% transform into a function with yyyy as parameter
%
% fixer les bornes lat, long (axes) et lai (z) pour pouvoir comparer
% les images de lai=f(lat,long) dans le temps
%
% utilisation de datetick avec octave > 3.0.2 et matlab
%
% compute sum, mean and std within box
%
% evaluation of interpolation quality of griddata
% may be usage of other method than linear
%
% IO error
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20101213T114841Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * save data (sum, mean, stddev) in ASCII format
%   * save figures
%   * add year in titles
%
% - fplod 20101208T095815Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * restruction because of griddata problem running octave
%
%   * 4 plots/page
%
%   * explicit types
%
%   * usage of :func:`griddata`
%     thanks to "Using MESHGRID and GRIDDATA to Fit Vector Data and Plot
%     Unevenly Spaced Data"
%     http://www.mathworks.com/support/tech-notes/1200/1212.html
%
%   * add info of Matlab Statiscal tool box (see
%
% - fplod 20101206T112730Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * usage of read_lai_2d
%
% - fplod 20100923T144224Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * replace simul_lai by simul_lai_2d
%   * replace read_lai by read_lai_2d
%   * add call to showgrid_lai
%
% - fplod 20100830T160653Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add argument to simul_lai call
%   * replace nb_date with nb_datestr
%
% - fplod 20100824T143715Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add header
%   * correction for reading LAI ASCII files
%     Aim is to have arrays doy(nb_doys), lat(nb_lat), lon(nb_lon) lai(nb_doy,nb_lat,nb_lon)
%     thanks to http://web.cecs.pdx.edu/~gerry/MATLAB/plotting/examples/readColData.m
%     and http://enacit1.epfl.ch/cours_matlab/mfiles.html#entrees_sorties
%
%   * make work running matlab : change ! operator to ~, replace " by ', endwhile by while
%
%   * improve indentation (including comments)
%
% - soclod 20100822
%
%   * draft
%-
%
global application;
global application_version;
global PROJECT_OD;
%
% default yyyy
yyyy=int16(2001);
yyyy_d=double(yyyy);
%
[lon_value, lat_value, datestr_value, lai_value]=read_lai_2d(yyyy);
%
% uncommment following lines to use simulated data
%index_simulation=int8(1);
%[lon_value, lat_value, datestr_value, lai_value]=simul_lai_2d(index_simulation);
%
% initialise les longitudes et les latitudes
lonmax=15.929;
lonmin=15.921;
latmax=-13.402;
latmin=-10.455;
% ++ forcage au domaine initial
lonmax=max(lon_value);
lonmin=min(lon_value);
latmax=max(lat_value);
latmin=min(lat_value);
%
% initialisation du nb de pas en latitude et longitude
lonres=10;
latres=10;
%
% Define the range and spacing of the x- and y-coordinates,
% and then fit them into X and Y
xv = linspace(lonmin, lonmax, lonres);
yv = linspace(latmin, latmax, lonres);
[xinterp,yinterp] = meshgrid(xv,yv);
clear xv
clear yv
%
% initialise le tableau LAI et des tableau associés
% tracé de la grille lat-lon
% mode='line'
% ++ ok result=showgrid_lai(lon_value, lat_value, mode);
%
nb_datestr=size(datestr_value,1);
%
ifigure=0;
%
switch application
   case {'matlab'}
      % plot lai_value=f(lat,long) for each day of year
      % definition du nb de subplots par page
      nb_column_page=2;
      nb_line_page=2;
      nb_subplot=nb_column_page * nb_line_page;
      %
      % loop on every date
      index_subplot=1;
      for index_datestr=1:nb_datestr
          z=squeeze(lai_value(:,index_datestr));
          % test if not only NaN for this day of year
          if ( isnan(z) )
             disp(['no data for ' datestr_value(index_datestr,:) ]);
          else
             disp(['data for '    datestr_value(index_datestr,:) ]);
             % début de figure si index_subplot = 1
             if (index_subplot == 1)
                ifigure=ifigure+1;
                figure(ifigure);
             end
             subplot(nb_line_page,nb_column_page,index_subplot);
             % Calculate Z in the X-Y interpolation space, which is an
             % evenly spaced grid:
             method='linear';
             zinterp = griddata(lon_value, lat_value, z, xinterp, yinterp, method);
             %mesh(xinterp,yinterp,zinterp)
             contourf(xinterp,yinterp,zinterp)
             xlabel('longitude');
             ylabel('latitude');
             title(['LAI ' datestr_value(index_datestr,:)]);
             %++colormap(cool(8))
             colorbar('WestOutside');
             % préparation subplot suivant : si index_subplot = nb_subplot,
             % remise à 1 de compteur de subplots
             if (index_subplot == nb_subplot )
                index_subplot=1;
             else
                index_subplot=index_subplot+1;
             end
          end
          % clear z
          % clear zinterp
      end
   case {'octave'}
      disp([mfilename ' pb with griddata with octave < 3.3.52']);
   otherwise
      error('unknown application running');
end
%
% calcul de la somme, de la moyenne et de l'écart-type des LAI
% pour chaque date sur toute la zone géographique
som=zeros(nb_datestr,1);
som(:)=NaN;
moy=zeros(nb_datestr,1);
moy(:)=NaN;
ecartyp=zeros(nb_datestr,1);
ecartyp(:)=NaN;
%
for index_datestr=1:nb_datestr
   som(index_datestr)=nansum(nansum(lai_value(:,index_datestr)));
   moy(index_datestr)=nanmean(nanmean(lai_value(:,index_datestr)));
   ecartyp(index_datestr)=nanstd(nanstd(lai_value(:,index_datestr)));
   disp(['index_datestr som moy std ', num2str(index_datestr), ' ', num2str(som(index_datestr)), ' ',  num2str(moy(index_datestr)), ' ', num2str(ecartyp(index_datestr))]);
end

% initialise les jours (les valeurs de x, les valeurs de y étant les sommes, les moyennes,
% les ecart type calculées
startdate=datenum(yyyy_d,1,1);
enddate=datenum(yyyy_d,12,eomday(yyyy_d,12));
xdata=zeros(nb_datestr,1);
xdata(:)=NaN;
for index_datestr=1:nb_datestr
    xdata(index_datestr)=datenum(datestr_value(index_datestr,:),'yyyymmdd');
end
xticklabel=datestr(xdata,'yyyymmdd');

% dessin de l'évolution de la somme des LAI
ifigure=ifigure+1;
figure(ifigure);
plot(xdata,som,'o');
set(gca,'XTickLabel',xticklabel)
xlabel('date')
ylabel('sum LAI 8d');
title(['sum LAI 8d of the year ', num2str(yyyy)]);
%
% save figure
printer='ps';
print_printer=['-d', printer];
fullfilename=[PROJECT_OD, mfilename, '_sum.', printer];
clear printer;
print(print_printer,fullfilename);
clear print_printer;
clear fullfilename;
%
% dessin de l'évolution de la moyenne des LAI
ifigure=ifigure+1;
figure(ifigure);
plot(xdata,moy,'o');
set(gca,'XTick',xdata)
set(gca,'XTickLabel',xticklabel)
xlabel('date')
ylabel('mean LAI 8d');
title(['mean LAI 8d of the year ', num2str(yyyy)]);
%
% save figure
printer='ps';
print_printer=['-d', printer];
fullfilename=[PROJECT_OD, mfilename, '_mean.', printer];
clear printer;
print(print_printer,fullfilename);
clear print_printer;
clear fullfilename;
%
% dessin de l'évolution des écart-types des LAI
ifigure=ifigure+1;
figure(ifigure);
plot(xdata,ecartyp,'o');
set(gca,'XTick',xdata)
set(gca,'XTickLabel',xticklabel)
xlabel('date')
ylabel('std LAI 8d');
title(['std LAI 8d of the year ', num2str(yyyy)]);
% save figure
printer='ps';
print_printer=['-d', printer];
fullfilename=[PROJECT_OD, mfilename, '_stddev.', printer];
clear printer;
print(print_printer,fullfilename);
clear print_printer;
clear fullfilename;
%

% sauvegarde des infos en ASCII
fullfilename=[PROJECT_OD, mfilename, '.txt'];
disp(['iii : ', mfilename ' : ','creation of ', fullfilename]);
fid=fopen(fullfilename,'w');
% écriture en tête
fprintf(fid,'# %s\n',fullfilename);
clear fullfilename;
% écriture data
fprintf(fid,'# yyyymmdd mean sum stddev\n');
for index_datestr=1:nb_datestr
    fprintf(fid,'%s %f %f %f\n', datestr_value(index_datestr,:), moy(index_datestr), som(index_datestr), ecartyp(index_datestr));
end
fclose(fid);
whos