source: trunk/src/nonprecipperiods.m @ 327

function [meantimestepof0, stdtimestepof0] = nonprecipperiods(ij_concat, time_concat, Longitudes, Latitudes, nomvar,  infocompl, minmean, maxmean, minstd, maxstd)

%NONPRECIPPERIODS est une fonction permettant de suivre les periodes de
%pauses dans les precipitations et de tracer la carte des moyennes des
%jours de non precipitation dans une pause ainsi que la carte des
%ecart-types.

%
%+
%
% ==================
% nonprecipperiods.m
% ==================
%
% .. function:: nonprecipperiods(ij_concat, time_concat, Longitudes, Latitudes, nomvar, infocompl, minmean, maxmean, minstd, maxstd)
%
% DESCRIPTION
% ===========
%
% - Donnees d'entrees :
%
%   * ij_concat : Matrice des donnees de precipitations dont les dimensions
%     sont temps, longitudes, latitudes.
%   * time_concat : Variable temporelle des donnees de precipitation en 
%     jour julien.
%   * Longitudes : Variable de la seconde dimension des donnees, les
%     longitudes.
%   * Latitudes : Variable de la seconde dimension des donnees, les
%     latitudes.
%   * nomvar (en caracteres) : Nom de la variable de precipitation, source
%     des donnees par exemple ('cumuls EPSAT (mm)' ou 'cumuls par analyses 
%     des Tb MSG (mm)').
%   * infocompl (en caracteres, factultatif) : Informations complementaires 
%     sur les donnees pour le titre. (La date etant automatiquement indique
%     dans le titre, informations suite a une operation special, 'echelle
%     restreinte entre 0 et 7' par exemple).
%   * minmean (Facultatif) : Borne inferieure de l'echelle des moyennes.
%   * maxmean (Facultatif) : Borne superieure de l'echelle des moyennes.
%   * minstd (Facultatif) : Borne inferieure de l'echelle des ecart-types.
%   * maxstd (Facultatif) : Borne inferieure de l'echelle des ecart-types.
%
% - Donnees de sorties :
%
%   * meantimestepof0 : Matrice comportant les donnees moyennes de jours de
%     non precipitations et dont les dimensions sont longitudes et 
%     latitudes.
%   * stdtimestepof0 : Matrice comportant les donnees des ecart-types des 
%     jours de non precipitations et dont les dimensions sont longitudes et 
%     latitudes.
%
% Cette fonction permet de suivre les periodes de pauses dans les 
% et de tracer la carte des moyennes des jours de non precipitation dans 
% une pause ainsi que la carte des ecart-types.
%
% EXAMPLES
% ========
%
% Voir la 'demo'.
%
% SEE ALSO
% ========
%
% :ref:`data_epsat`
%
% :ref:`data_msg`
%
% :func:`cumulEPSAT_concat`
% :func:`intensityMSG_concat`
%
% TODO
% ====
%
% describe demo
%
% demo with simulated data
%
% cratos octave 3.4 : warning: ft_render: skipping missing glyph for character `
%
% should long text be written on the figure ? now : pb of adjustment between plot size and text length
%
% fix colorbar limits to make comparison between matlab/octave and machine easier
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20110804T132212Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * mimize decoration (but not enougth ?)
%
% - jaclod 2011-07-29
%
%   * Ajout des donnees de sorties.
%   * Ajout d'une restriction sur les echelles de moyennes et ecart-types
%     pour une meilleure observation.
%   * Ajout de la demonstration.
%   * Ajout des cotes.
%
% - jaclod 2011-07-28
%
%   * Ajout de la suppression de la variable 'timestepof0' qui induit des
%     erreurs sinon (il peut rester des valeurs du pixel precedent).
%   * Suppression des autres variables egalement.
%   * Suppression des pas de temps ou les donnees sont des 'NaN'.
%   * Modifications des titres.
%
% - jaclod 2011-07-27
%
%   * Creation.
%
%-

%%%%%%%%%%%%%%%%%
% Initialisation
%%%%%%%%%%%%%%%%%

% Initialisation de la matrice qui contiendra les moyennes du temps en jour
% des periodes de non precipitations.
meantimestepof0 = zeros(size(ij_concat,2),size(ij_concat,3));

% Initialisation de la matrice qui contiendra les ecart-types du temps en
% jour des periodes de non precipitations.
stdtimestepof0 = zeros(size(ij_concat,2),size(ij_concat,3));

% Suppressions des donnees pour les jours ou nous avons que des 'NaN'.
% On cherche s'il y en a, quels sont les pas de temps concernes et leur
% nombre.
indnan = find(isnan(ij_concat(:,1,1)));
nbnan = size(indnan,1);

% S'il y en a, on va les supprimer.
if ~isempty(indnan);
    
    % On boucle pour supprimer chaque pas de temps concernes en commencant
    % par les derniers (si on souhaite par exemple supprimer les pas de
    % temps 2 et 5 et que l'on supprime le 2 d'abord, le 5e pas de temps
    % sera ensuite decale au 4e. En commencant par la fin, cela ne se
    % produit pas).
    for nansuppr = 1 : nbnan;
       
        ij_concat(indnan(nbnan - nansuppr + 1),:,:) = [];
        
        time_concat(indnan(nbnan - nansuppr + 1)) = []; 
        
    end
    
end

clear('indnan','nbnan','nansuppr');

%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%
% Boucle sur les pixels
%%%%%%%%%%%%%%%%%%%%%%%%

% Boucle sur les longitudes.
for ilon = 1 : size(ij_concat,2);
    
    % Boucle sur les latitudes.
    for ilat = 1 : size(ij_concat,3);
        
        % Recuperation des donnees du pixel.
        data = squeeze(ij_concat(:,ilon,ilat));
        
        % Recuperation des indices ou les precipitations sont nulles
        % pour le pixel considere.
        indzero = find(data==0);
        
        % Boucle sur les indices trouves (Elle ne se fait pas si aucun
        % 0 n'a ete trouve).
        indindzero = 1;
        nbperiods = 0;
        while indindzero <= size(indzero,1)
            
            % Recuperation de l'indice de debut de la periode de 0.
            indstart = indzero(indindzero);
            
            % On cherche l'indice de fin de la periode de 0.
            while indindzero ~= size(indzero,1) && indzero(indindzero+1) - indzero(indindzero) == 1;
                
                indindzero = indindzero+1;
                
            end
            
            indstop = indzero(indindzero);
            
            nbperiods = nbperiods + 1;
            
            % On enregistre pour chaque periode, le nombre de jours de non
            % precipitation.
            timestepof0(nbperiods) = 1 + indstop - indstart;
            
            indindzero = indindzero + 1;
            
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%
        % Moyenne et ecart-type
        %%%%%%%%%%%%%%%%%%%%%%%%
        
        % Calcul de la moyenne.
        if ~isempty(indzero);
            ij_meantimestepof0 = sum(timestepof0) / nbperiods;
        else
            ij_meantimestepof0 = 0;
        end
        
        meantimestepof0(ilon,ilat) = ij_meantimestepof0;
        
        % Calcul de l'ecart-type.
        if ~isempty(indzero)
            ij_stdtimestepof0 = sqrt(sum(power(timestepof0,2))/nbperiods - power(ij_meantimestepof0,2));
        else
            ij_stdtimestepof0 = 0;
        end
        stdtimestepof0(ilon,ilat) = ij_stdtimestepof0;
        
        clear('ij_meantimestepof0','ij_stdtimestepof0','indstart','indstop','indzero','indindzero','timestepof0','data');
        
        %%%%%%%%%%%%%%%%%%%%%%%%
        
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Images des moyennes et ecart-types
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Ouverture d'une figure.
figure();

% Creation de variables comportant les titres et noms des axes.
if nargin >= 6;
    Titlegen = (['Suivi des pauses de precipitations d''apres les ' nomvar ' sur ' num2str(size(ij_concat,1)) ' jours du ' datestr(time_concat(1),'dd mmm yyyy') ' au ' datestr(time_concat(size(time_concat,2)),'dd mmm yyyy') char(10) infocompl]);
else
    Titlegen = (['Suivi des pauses de precipitations d''apres les ' nomvar ' sur ' num2str(size(ij_concat,1)) ' jours du ' datestr(time_concat(1),'dd mmm yyyy') ' au ' datestr(time_concat(size(time_concat,2)),'dd mmm yyyy')]);
end

Title1 = ('Moyennes des jours de non precipitation');
Title2 = ('Ecart-types des jours de non precipitation');
Abs = 'Longitudes';
Ord = 'Latitudes';

% Chargement des donnees cotieres.
load('cmgco_traitdecote.mat','Xgco','Ygco');

% Moyennes (la transposee est utilisee pour mettre les longitudes selon
% l'axe x).
subplot(1,2,1);
imagesc(Longitudes,Latitudes,meantimestepof0');
set(gca, 'DataAspectRatio', [1 1 1]);
set(gca,'ydir','normal');
title(Title1);
xlabel(Abs)
ylabel(Ord)
if nargin >= 8;
    caxis([minmean maxmean])
end
colorbar
hold on

hold on
% Ajout des cotes.
plot(Xgco,Ygco,'k.');

% Ecart-types (la transposee est egalement utilise ici).
subplot(1,2,2);
imagesc(Longitudes,Latitudes,stdtimestepof0');
set(gca, 'DataAspectRatio', [1 1 1]);
set(gca,'ydir','normal');
title(Title2);
xlabel(Abs)
ylabel(Ord)
if nargin == 10;
    caxis([minstd maxstd])
end
colorbar

hold on
% Ajout des cotes.
plot(Xgco,Ygco,'k.');

text(-0.2,1.2,Titlegen,...
    'HorizontalAlignment','center','verticalAlignment','middle')

end

%!demo
%! clear all
%! close all
%! varamma_startup
%! more off
%! [ij_cumul_concat, time_concat, Longitudes, Latitudes] = cumulEPSAT_concat('extracted', 0, 0, 0, 0, 0, 0, 2006);
%! [meantimestepof0, stdtimestepof0] = nonprecipperiods(ij_cumul_concat, time_concat, Longitudes, Latitudes, 'cumuls EPSAT (mm)',  'Restriction des echelles entre 0 et 7 pour la demonstration', 0, 7, 0, 7);