source: trunk/src/intensityMSG_concat.m @ 327

function [ij_intensity_concat, time_concat, Longitudes, Latitudes] = intensityMSG_concat(ncfilesset, lonmin, lonmax, latmin, latmax, longresol, latresol, Annee, Mois_start, Jour_start, Mois_stop, Jour_stop)

%INTENSITYMSG_CONCAT permet d'analyser les donnees de temperature de
%brillance par MSG et de concatener les intensites resultantes en sortie.
%L'analyse s'effectue sur une seule saison, sur tous les fichiers ou entre
%deux dates predefinies.

%
%+
%
% ======================
% intensityMSG_concat.m
% ======================
%
% .. function:: intensityMSG_concat(ncfilesset, lonmin, lonmax, latmin, latmax, longresol, latresol, Annee, Mois_start, Jour_start, Mois_stop, Jour_stop)
%
% DESCRIPTION
% ===========
%
% - Donnees d'entrees :
%
%   * ncfilesset (en caracteres) : Type de fichier NetCDF que l'on cherche,
%     'normal' ou 'extracted'.
%   * lonmin : Borne inferieure des longitudes du domaine etudie.
%   * lonmax : Borne superieure des longitudes du domaine etudie.
%   * latmin : Borne inferieure des latitudes du domaine etudie.
%   * latmax : Borne superieure des latitudes du domaine etudie.
%   * longresol : Nouvelle resolution souhaitee (nombre de points) pour les
%     longitudes (0 pour garder la resolution d'origine).
%   * latresol : Nouvelle resolution souhaitee (nombre de points) pour les
%     latitudes (0 pour garder la resolution d'origine).
%   * Annee : Annee de la saison d'etude.
%   * Mois_start (Facultatif) : Mois du premier jour d'etude pour la
%     selection des fichiers.
%   * Jour_start (Facultatif) : Jour du premier jour d'etude pour la
%     selection des fichiers.
%   * Mois_stop (Facultatif) : Mois du dernier jour d'etude pour la
%     selection des fichiers.
%   * Jour_stop (Facultatif) : Jour du dernier jour d'etude pour la
%     selection des fichiers.
%
% - Donnees de sorties :
%
%   * ij_intensity_concat : Matrice de donnees (temps, longitudes,
%     latitudes)contenant les intensites concatenees des differents jours
%     d'etude.
%   * time_concat : Variable de temps, en jour julien, pour identifier les
%     jours etudies.
%   * Longitudes : Variable de longitudes dans le domaine selectionne.
%   * Latitudes : Variable de latitudes dans le domaine selectionne.
%
% Cette fonction permet d'analyser les donnees de temperature de brillance
% par MSG et de concatener les intensites resultantes en sortie. L'analyse
% s'effectue sur une seule saison, sur tous les fichiers ou entre deux
% dates predefinies.
%
% EXAMPLES
% ========
%
% Voir la 'demo'.
%
% SEE ALSO
% ========
%
% :ref:`data_msg`
%
% :ref:`varamma_startup.m`
%
% :func:`MSGbuildfullfilename`
% :func:`extractedmatrixdata`
% :func:`spatialresolution`
% :func:`MSGread`
% :func:`intensityMSG_day`
%
% :func:`int2precip`
% :func:`nonprecipperiods`
% :func:`MSGanalysiscrossvalidation`
% :func:`missingdaytonan`
%
% TODO
% ====
%
% coding rule and description of demo
%
% demo normal
%
% demo simul
%
% use :func:`simul_msg`
%
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20110805T093216Z cratos (Linux)
%
%   * remove varamma_startup call

% - fplod 20110804T150823Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * octave and matlab 7.4 (zeus) do not like [argout1,~, argout3]=myfunc() so
%     replace ~ by temporary variable
%
% - jaclod 2011-07-27
%
%   * Revision de la documentation et des commentaires.
%
% - jaclod 2011-07-22
%
%   * Compatibilite avec des fichiers 'extracted'.
%
% - jaclod 2011-07-20
%
%   * Passage sur Cratos, fonctionne!
%   * Modifications de la documentation et des commentaires.
%
% - jaclod 2011-07-19
%
%   * Ajout de boucles pour la lecture des fichiers avec ou sans selection.
%
% - jaclod 2011-07-13
%
%   * Creation.
%
%-

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

% Appel de la variable global 'PROJECT_ID'.
global PROJECT_ID

% Sauvegarde de la variable 'Annee' sous forme de caracteres.
annee = num2str(Annee);

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Boucles sur les fichiers
%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Les boucles se realisent de differentes facons selon la selection ou non
% de jours.

% Si le nombre d'arguments est superieur a 8 c'est qu'une selection est
% demandee. Sinon, on analyse tous les fichiers present dans le repertoire
% MSG.
if nargin > 8;

    % Creation d'un vecteur comportant le nombre maximum de jours que l'on
    % a dans chaque mois.
    maxjourmois = [31 29 31 30 31 30 31 31 30 31 30 31];

    %%%%%% Boucle sur les mois %%%%%%
    for Mois = Mois_start : Mois_stop;

        % Recuperation de la variable 'Mois' sous forme de caracteres.
        if Mois < 10;
            mois = ['0' num2str(Mois)];
        else
            mois = num2str(Mois);
        end

        % S'il s'agit du premier mois, le jour de depart sera celui entre
        % en argument, sinon ce sera le premier du mois.
        if Mois == Mois_start;
            Jour_ini = Jour_start;
        else
            Jour_ini = 1;
        end

        % S'il s'agit du dernier mois, le jour de fin sera celui entre
        % en argument, sinon ce sera le dernier jour possible du mois.
        if Mois == Mois_stop;
            Jour_fin = Jour_stop;
        else
            Jour_fin = maxjourmois(Mois);
        end

        %%%%%% Boucle sur les jours %%%%%%
        for Jour = Jour_ini : Jour_fin;

            % Recuperation de la variable 'Jour' sous forme de caracteres.
            if Jour < 10;
                jour = ['0' num2str(Jour)];
            else
                jour = num2str(Jour);
            end

            % Creation du chemin d'acces complet au fichier correspondant.
            [fullfilename] = MSGbuildfullfilename(ncfilesset, annee, mois, jour);

            %%%%%% Analyse du fichier %%%%%%

            % Verification que le fichier du jour existe bien. Si ce n'est
            % pas le cas un message d'avertissement est donne. Sinon, on
            % peut faire l'analyse.
            % Cette verification differe selon le type de fichier en raison
            % d'un changement de repertoire et de nom.

            switch ncfilesset

                case 'normal'
                    repertoire = [PROJECT_ID 'MSG/' annee '/' mois '/'];

                case 'extracted'
                    repertoire = [PROJECT_ID 'MSG/extracted/' annee '/' mois '/'];

                otherwise
                    error('La variable ncfilesset doit etre ''normal'' ou ''extracted''.');

            end

            path(path,repertoire);

            switch ncfilesset

                case 'normal'
                    if exist(['msg-tb108_' annee '-' mois '-' jour '_15min.nc'],'file') == 2;

                        % Message confirmant l'analyse du fichier.
                        disp(['iii : Analyse du fichier msg-tb108_' annee '-' mois '-' jour '_15min.nc.']);

                        % Lecture du fichier NetCDF, recuperation des
                        % donnees et des variables.
                        [TbMSG, Temps, Longitudes, Latitudes] = MSGread(ncfilesset, fullfilename);

                        % Extraction des donnees.
                        [TbMSG, Temps, Longitudes, Latitudes] = extractedmatrixdata(TbMSG, Temps, Longitudes, Latitudes, Temps(1), Temps(size(Temps,1)), lonmin, lonmax, latmin, latmax);

                        %%%%%% Changement de resolution %%%%%%
                        % Si longresol et/ou latresol ont pour valeur 0, alors on
                        % recupere la valeur de la resolution initiale pour ne pas
                        % faire de changement.
                        if longresol == 0;
                            longresol = size(TbMSG,2);
                        end

                        if latresol == 0;
                            latresol = size(TbMSG,3);
                        end

                        % Si la resolution demandee est differente de celle
                        % initiale, alors on interpole les donnees avec la
                        % fonction 'spatial resolution'.
                        if longresol ~= size(TbMSG,2) || latresol ~= size(TbMSG,3);
                            [Longitudes, Latitudes, TbMSG] = spatialresolution(Longitudes, Latitudes, TbMSG, longresol, latresol);
                        end

                        %%%%%% Calcul des intensites %%%%%%
                        [ij_intensity, unused_Tps_SeuilInf, unused_Tb_SeuilSup, unused_Tb_LimNoyau] = intensityMSG_day(Temps, TbMSG, 3, 235, 190);
                        clear unused_Tps_SeuilInf
                        clear unused_Tb_SeuilSup
                        clear unused_Tb_LimNoyau

                        %%%%%% Concatenation %%%%%%
                        % Sauvegarde des intensites (la dimension 3 est ici
                        % le temps).
                        if Mois == Mois_start && Jour == Jour_start;
                            ij_intensity_concat = ij_intensity;
                        else
                            ij_intensity_concat = cat(3,ij_intensity_concat,ij_intensity);
                        end

                        % Sauvegarde de l'information sur la date. On la
                        % sauve en jour julien, il faut donc convertir par
                        % rapport au temps fourni, en heure depuis le 1er
                        % janvier 1960.
                        Tempsjourjulien = Temps(1) / 24 + datenum(1960,01,01);
                        if Mois == Mois_start && Jour == Jour_start;
                            time_concat = Tempsjourjulien;
                        else
                            time_concat = cat(2,time_concat,Tempsjourjulien);
                        end

                    else
                        % Message d'avertissement en cas de non presence du
                        % fichier.
                        disp(['www : le fichier msg-tb108_' annee '-' mois '-' jour '_15min.nc n''est pas present.']);
                    end

                case 'extracted'
                    if exist(['extracted-msg-tb108_' annee '-' mois '-' jour '_15min.nc'],'file') == 2;

                        % Message confirmant l'analyse du fichier.
                        disp(['iii : Analyse du fichier extracted-msg-tb108_' annee '-' mois '-' jour '_15min.nc.']);

                        % Lecture du fichier NetCDF, recuperation des
                        % donnees et des variables.
                        [TbMSG, Temps, Longitudes, Latitudes] = MSGread(ncfilesset, fullfilename);

                        % Pas d'extraction des donnees.

                        %%%%%% Changement de resolution %%%%%%
                        % Si longresol et/ou latresol ont pour valeur 0,
                        % alors on recupere la valeur de la resolution
                        % initiale pour ne pas faire de changement.
                        if longresol == 0;
                            longresol = size(TbMSG,2);
                        end

                        if latresol == 0;
                            latresol = size(TbMSG,3);
                        end

                        % Si la resolution demandee est differente de celle
                        % initiale, alors on interpole les donnees avec la
                        % fonction 'spatial resolution'.
                        if longresol ~= size(TbMSG,2) || latresol ~= size(TbMSG,3);
                            [Longitudes, Latitudes, TbMSG] = spatialresolution(Longitudes, Latitudes, TbMSG, longresol, latresol);
                        end

                        %%%%%% Calcul des intensites %%%%%%
                        [ij_intensity, unused_Tps_SeuilInf, unused_Tb_SeuilSup, unused_Tb_LimNoyau] = intensityMSG_day(Temps, TbMSG, 3, 235, 190);
                        clear unused_Tps_SeuilInf
                        clear unused_Tb_SeuilSup
                        clear unused_Tb_LimNoyau

                        %%%%%% Concatenation %%%%%%
                        % Sauvegarde des intensites (la dimension 3 est ici
                        % le temps).
                        if Mois == Mois_start && Jour == Jour_start;
                            ij_intensity_concat = ij_intensity;
                        else
                            ij_intensity_concat = cat(3,ij_intensity_concat,ij_intensity);
                        end

                        % Sauvegarde de l'information sur la date. On la
                        % sauve en jour julien, il faut donc convertir par
                        % rapport au temps fourni, en heure depuis le 1er
                        % janvier 1960.
                        Tempsjourjulien = Temps(1) / 24 + datenum(1960,01,01);
                        if Mois == Mois_start && Jour == Jour_start;
                            time_concat = Tempsjourjulien;
                        else
                            time_concat = cat(2,time_concat,Tempsjourjulien);
                        end

                    else
                        % Message d'avertissement en cas de non presence du
                        % fichier.
                        disp(['www : le fichier extracted-msg-tb108_' annee '-' mois '-' jour '_15min.nc n''est pas present.']);
                    end

                otherwise
                    % Cas impossible. Erreur deja donnee.
            end

        end

    end

else

    %%%%% Recherche des mois de donnees %%%%%
    % Repertoire correspondant aux fichiers MSG de l'annee en question.
    % differents selon le type de fichier.
    switch ncfilesset

        case 'normal'
            repertoireMSGsaison = [PROJECT_ID 'MSG/' annee];

        case 'extracted'
            repertoireMSGsaison = [PROJECT_ID 'MSG/extracted/' annee];

        otherwise
    end

    path(path,repertoireMSGsaison);

    % On cherche les dossier dans ce dernier repertoire correspondant aux
    % mois dont nous disposons des donnees.
    listdir = dir(repertoireMSGsaison);

    % Nombre de dossier. Attention les deux premiers resultats sont '.' et
    % '..' et ne compte donc pas comme dossier de fichiers.
    nbdir = numel(listdir)-2;

    %%%%% Boucle sur les mois %%%%%
    for numdir = 1 : nbdir;

        % Recuperation du numero du mois sous forme de caracteres.
        mois = listdir(numdir+2).name;

        % Repertoire du mois correspondant.
        repertoireMSGmois = [repertoireMSGsaison '/' mois];

        path(path,repertoireMSGmois);

        %%%%% Recherche des fichiers netCDF %%%%%
        ext = '*.nc';
        chemin = fullfile(repertoireMSGmois,ext);
        listncfiles = dir(chemin);

        nbfiles = numel(listncfiles);

        %%%%% Boucle sur les fichiers, les jours %%%%%
        for numfiles = 1 : nbfiles;

            % Recuperation du nom du fichier et de son chemin d'acces
            % complet.
            fullfilename = [repertoireMSGmois '/' listncfiles(numfiles).name];

            % Message confirmant l'analyse du fichier.
            disp(['iii : Analyse du fichier ' listncfiles(numfiles).name '.']);

            % Lecture du fichier NetCDF, recuperation des donnees et des
            % variables.
            [TbMSG, Temps, Longitudes, Latitudes] = MSGread(ncfilesset, fullfilename);

            switch ncfilesset

                case 'normal'
                    % Extraction des donnees.
                    [TbMSG, Temps, Longitudes, Latitudes] = extractedmatrixdata(TbMSG, Temps, Longitudes, Latitudes, Temps(1), Temps(size(Temps,1)), lonmin, lonmax, latmin, latmax);

                case 'extracted'
                    % Pas d'extraction.

                otherwise
            end

            %%%%%% Changement de resolution %%%%%%
            % Si longresol et/ou latresol ont pour valeur 0, alors on
            % recupere la valeur de la resolution initiale pour ne pas
            % faire de changement.
            if longresol == 0;
                longresol = size(TbMSG,2);
            end

            if latresol == 0;
                latresol = size(TbMSG,3);
            end

            % Si la resolution demandee est differente de celle
            % initiale, alors on interpole les donnees avec la fonction
            % spatial resolution.
            if longresol ~= size(TbMSG,2) || latresol ~= size(TbMSG,3);
                [Longitudes, Latitudes, TbMSG] = spatialresolution(Longitudes, Latitudes, TbMSG, longresol, latresol);
            end

            %%%%%% Calcul des intensites %%%%%%
            [ij_intensity, unused_Tps_SeuilInf, unused_Tb_SeuilSup, unused_Tb_LimNoyau] = intensityMSG_day(Temps, TbMSG, 3, 235, 190);
            clear unused_Tps_SeuilInf
            clear unused_Tb_SeuilSup
            clear unused_Tb_LimNoyau

            %%%%%% Concatenation %%%%%%

            % Sauvegarde des intensites (la dimension 3 est ici le temps).
            if numdir==1 && numfiles==1;
                ij_intensity_concat = ij_intensity;
            else
                ij_intensity_concat = cat(3,ij_intensity_concat,ij_intensity);
            end

            % Sauvegarde de l'information sur la date. On la sauve en jour
            % julien, il faut donc convertir par rapport au temps fourni,
            % en heure depuis le 1er janvier 1960.
            Tempsjourjulien = Temps(1) / 24 + datenum(1960,01,01);
            if numdir==1 && numfiles==1;
                time_concat = Tempsjourjulien;
            else
                time_concat = cat(2,time_concat,Tempsjourjulien);
            end

        end

    end

end

% On passe le temps en premiere dimension pour plus de concordance avec les
% autres variables.
ij_intensity_concat = permute(ij_intensity_concat,[3 1 2]);

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

end

%!demo
%! clear all
%! close all
%! varamma_startup;
%! more off
%! [ij_intensity_concat, time_concat, Longitudes, Latitudes] = intensityMSG_concat('normal',-20, 0, 0, 20, 200, 200, 2006, 07, 26, 08, 09);