source: trunk/PROGRAMMES/eof_NCEP_sais.m @ 26

%EOF_NCEP_SAIS calcule les EOFs et composantes principales du champs de temperature de surface en moyennes mensuelles climatologiques
% (12 points temporels)
%
% le programme enregistre automatiquement les resultats du calcul
% (attention: penser a changer le nom du fichier de sortie en cas de
% modification de parametres du programme: domaine spatial, moyenne
% temporelle effectuee: annee / mois, etc. parametre "nomfic" ci dessous)
%
% le programme permet egalement de visualiser les resultats. Voir le doc
% .pdf joint pour l'interpretation de cette visualisation.
%
% Date: Avril-Mai 2008
% Auteur: juliette.mignot@locean-ipsl.upmc.fr
% Contexte Ecole d'ete ENVI-STAT 2008

%+
%
% module
% ======
%
% ``eof_NCEP_sais``
%
% DESCRIPTION
% ===========
%
% ``eof_NCEP_sais`` lit les fichiers
% ``../DONNEES/lsmask_tropatl_PacandMedblancs.mat``
% et ``../DONNEES/skt.mon.tropatl.mat``.
%
% ``eof_NCEP_sais`` calcule ++blabla++
%
% ``eof_NCEP_sais`` sauvegarde dans 
% ``../DONNEES/eof_SSTmens_tropAtl30N-20S.mat``.
%
% Ces fichiers sont lus par ``reg_serietemporelles.m_``\ .
%
% ``eof_NCEP_sais`` affiche trois figures :
%  - ++
%  - ++
%  - valeurs propres de la matrice de covariance
%
% ``eof_NCEP_sais`` sauve ces images dans 
% ``./eof_SSTmens_tropAtl30N-20S_[123].ps``.
%
% ++
%
% EXAMPLES
% ========
%
% ::
%
% >> tpacpandreg_startup
% >> eof_NCEP_sais
%
% SEE ALSO
% ========
%
% tpacpandreg_startup.m_
%
% .. _tpacpandreg_startup.m : tpacpandreg_startup.m.html
%
% initfig.m_
%
% .. _initfig.m : initfig.m.html
%
% reg_seriestemporelles.m_
%
% .. _reg_seriestemporelles.m : reg_seriestemporelles.m.html
%
% whorldmap.m_
%
% .. _whorldmap.m : whorldmap.m.html
%
% colorbartype.m_
%
% .. _colorbartype.m : colorbartype.m.html
%
% contlab.m_
%
% .. _contlab.m : contlab.m.html
%
% TODO
% ====
%
% improve description
%
% tester la branche NetCDF
%
% EVOLUTIONS
% ==========
%
% $Id$
%
%  * remplacement de ``eps`` par ``ps` pour pouvoir imprimer les figures 
%    produites par ``octave``
%    Dans le fichier .eps produit par ``octave`` on pouvait lire 
%    ``!PS-Adobe-2.0 EPSF-2.0``.
%    Aucun message d'erreur n'apparaît suite à la demande d'impression 
%    ``zeus$ lpr -o job-sheets='none' eof_SSTmens_tropAtl30N-20S_1.octave.ps``, 
%    mais aucune feuille ne sort.
%    Pour info dans le fichier .eps produit par matlab, on pouvait lire
%    ``%!PS-Adobe-2.0 EPSF-1.2``. Aucun souci d'impression avec ce fichier !
%    En remplaçant ``eps`` par ``ps`, on n'a plus de souci d'impression.
%
%
% - jmignot 2009-08-25
%
%  * suppression d'appel à contlab et colorbartype sous octave
%
%
% - fplod 2009-08-25T10:40:16Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%  * modif commentaires liée au dernier chgt de Juliette
%
% - jmignot 2009-08-25
%
%  * chgt de noms de fichiers : ``EOF_ devient ``eof_``
%  * remplacement de ``colorbar`` par ``colorbartype``
%    parce pas beau ni en matlab ni en octave (dégradé, noir et blanc dans
%    les .jpg)
%  * contraintes sur les limites x et y de axes car matlab et octave ne
%    fixent pas les même limites par défaut
%  * utisation de ``contlab`` pour ++

%
% - fplod 2009-08-21T15:04:16Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * as octave save default format is ASCII and no default extension
%     while matlab save default format is MAT v5 mat-file (little endian) and
%     default extension is ``.mat``, write more precise save instruction
%
% - fplod 2009-08-21T12:58:38Z aedon.locean-ipsl.upmc.fr (Darwin)
%
%   * add header
%   * répertoire de données relatifs au répertoire courant
%   * suppression de ``clear all`` présent dans ``tpacpandreg_startup.m``
%   * contournement pb octave pour ``contourf`` :
%     Dans le fil de discussion [Pkg-octave-devel] Bug#492223: octave3.0: contourf only accepts vector X and Y if same+++,
%     je comprends que j'ai affaire à un bug de octave 3.0 (donc à surveiller
%     pour les prochaines versions) qui ne veut pas faire de contourf si les
%     variables X et Y n'ont pas la même dimension et c'est le cas ici
%     (lat = 42x1, lon= 70x1 vue avec size).
%     http://lists.alioth.debian.org/pipermail/pkg-octave-devel/2008-July/004672.html
%   * sauvegarde des figures

%
%-
%
N=4; %nb d'EOF a tracer

nomfic='eof_SSTmens_tropAtl'
% mettre 'eof_SSTan_tropAtl' quand on travaille en moyennes annuelles. Ce
% fichier sera rappele dans les programme de regression lineaire
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%              CHARGEMENT ET PREPARATION DES DONNEES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nannees=60;
   repertoire=[ '..' filesep 'DONNEES' filesep ] % a adapter a chaque cas: chemin d acces aux donnees

% Mask
%------
  file='lsmask_tropatl_PacandMedblancs.mat'
   load([repertoire file]);
  ny=size(mask,1); nx=size(mask,2); % dimensions spatiales des tableaux

% Data
%--------
  %SI .MAT
  repertoire2=[repertoire 'MAT' filesep] % a adapter a chaque cas: chemin d acces aux donnees
  file=['skt.mon.tropatl.mat']
  load([repertoire2 file]); % chargement des donnees


  %SI NETCDF
  %repertoire2=[repertoire 'NETCDF' filesep]% a adapter a chaque cas: chemin d acces aux donnees
  % file=['skt.mon.tropatl.nc']
  % f=netcdf([repertoire file],'nowrite');
  % tab=f{'skt'}(:,:,:);


   tab=permute(tab,[2 3 1]);% permutation des dimensions pour organiser le tableau sous la forme [lat lon temps]


%---------------------------------------------------
%  MANIPULATIONS TEMPORELLES: MOYENNE ET ANOMALIES
%---------------------------------------------------

% Construction de moyennes mensuelles saisonnieres
%---------------------------------------------------
tab=reshape(tab,size(tab,1),size(tab,2),12,nannees);
tab=squeeze(mean(tab,4));

%  Construction des anomalies: on soustrait en chque point la moyenne temporelle totale
%--------------------------------------------------------------------------
tab_mean=squeeze(mean(tab,3)); tab_mean=repmat(tab_mean,[1 1 size(tab,3)]);
anom=tab-tab_mean;
ntemp=size(anom,3);

%-------------------------------------------------------------
%  MANIPULATIONS SPATIALES: PONDERATION, ZOOM GEOGRAPHIQUE
%-------------------------------------------------------------

% Pondération
%----------------
%Si on veut éventuellement attribuer un poids différent aux différents
%points spatiaux, on le fait ici.

% Préparation du zoom géographique: sur le masque, j'identifie les points
% supplémentaires que je veux éliminer pour le calcul (critère en longitude
% et/ou latitude par exemple)
%---------------------------------------------------------------------
limite_latS=-20; %limite du domaine d'etuade en latitude
limite_latN=30; %limite du domaine d'etuade en latitude
lat2D=repmat(lat,[1 size(mask,2)]);
mask(lat2D>limite_latN&mask==0)=4;
mask(lat2D<limite_latS&mask==0)=4;
clear lat2D

% Application du masque terre/mer pour pouvoir ensuite enlever les points de terre: on ne veut
% travailler que sur l'océan
%---------------------------------------------------------------------
mask3D=repmat(mask,[1 1 size(tab,3)]);
tab(mask3D~=0)=NaN;

%  Elimination des points NaN
%--------------------------------
anom=reshape(anom,ny*nx,size(anom,3));
maskr=reshape(mask,ny*nx,1);
pt_val=find(maskr==0);
Z=zeros(length(pt_val),ntemp);
Z=anom(pt_val,:);
%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                DECOMPOSITION EN COMPOSANTES PRINCIPALES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%      CALCUL DE LA MATRICE DE COVARIANCE ET DIAGONALISATION

% Z est donc une matrice N(points spatiaux) x M(points temporels) de
% moyenne temporelle nulle sur laquelle on va faire l'analyse en EOF

% matrice de covariance:
  C = Z * Z';
%  eval(['save ' repertoire 'MatCov_SST_tropAtl' num2str(limite_latN) 'N' num2str(limite_latS) 'S.asc C -ascii'])

% diagonalisation
  [E,Lambda] = eig(C);
   % Lambda est la matrice diagonale des valeurs propres et E la matrice dont les colonnes sont les vecteurs propres correspondants.
   % Soit Z*E = E*Lambda
   % Les vecteurs propres sont appelés les EOF. Chaque EOF peut être considéré comme une carte de l'espace physique.


 %%
 %%%% calcul de la matrice donnant les pourcentages de variance
Lambda = diag(Lambda);
pctg_var = round(Lambda / sum(Lambda) * 1000) / 10;


 %On definit les composantes principales de Q comme la decomposition du signal decrit par Q sur les EOF E.

   P = Z' * E;

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%             RECONSTRUCTION DES CARTES PHYSIQUES ET NORMALISATION DES
 %%%             COMPOSANTES PRINCIPALES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%initialisation des tableaux de stockage finaux
TABEOF=zeros(ny,nx,N); %TABLEAU FINAL DES CARTES PHYSIQUES
PC=zeros(ntemp,N);    %TABLEAU FINAL DES COMPOSANTES PRINCIPALES
p_var=zeros(N,1);    %TABLEAU FINAL DES COMPOSANTES PRINCIPALES

% boucle sur les differentes EOF pour reconstruire les cartes, normaliser
% et enregistrer dans le tableau final
for  i=1:1:N

  disp(strcat('normalisation, projection et trace de l eof',num2str(i)))

  pc = normalisation(P(:,end-i+1));
  %
  tabeof_temp=E(:,end-i+1);
  tabeof=NaN+zeros(ny*nx,1);
  tabeof(pt_val)=tabeof_temp;
  tabeof=reshape(tabeof,ny,nx);
  %remarque: autre facon de calculer le tabeof:
  %  tabeof=projection(pc,Q'); tabeof=reshape(tabeof,ny,nx);

      %%% enregistrement
  TABEOF(:,:,i)=tabeof;
  PC(:,i)=pc;
  p_var(i)=pctg_var(end-i+1);

end;

fullfilename=[repertoire nomfic num2str(limite_latN) 'N' num2str(limite_latS) 'S.mat'];
if (run_octave == 0)
 save(fullfilename,'TABEOF','PC','p_var','lon','lat','mask','-v7');
else
 save('-v7',fullfilename,'TABEOF','PC','p_var','lon','lat','mask');
end
clear fullfilename;

 %%               TRACE

close all
%parametres geographiques des graphes
  geo=[-100 20 -30 30]; %[lonmin lonmax latmin latmax];
  yti=[-30:10:30]; xti=[-100:30:20]; %labels de latitude et de longitude
%niveaux de couleurs
  lev=[-2:.02:2];
  clip_lev=[-0.06 0.06];

figure(1);
initfig;
figure(2);
initfig;

for i=1:4

   %subplot(4,2,2*i-1)
   ifigure=1;
   figure(ifigure);
   subplot(2,2,i)
   if (run_octave == 0)
    x=lon;
    y=lat;
   else
    [x,y]=meshgrid(lon,lat);
   end
   z=squeeze(TABEOF(:,:,i));
   contourf(x,y,z,lev);
   caxis(clip_lev);
   hold on
   whorldmap
   % limites geographies et labels de longitude et de latitude
  ylim([geo(3) geo(4)])
  xlim([geo(1) geo(2)])
  set(gca, 'Xtick', xti, 'Ytick', yti)
  if (run_octave == 0)
    [longstr, latstr] = contlab(xti, yti);
    set(gca, 'Xticklabel', longstr, 'Yticklabel', latstr)
  end;
  xlabel('longitude')
  ylabel('latitude')

   %titre de la figure
   title([num2str(p_var(i)) '%'])

  %barre de couleur
  if (run_octave == 0)
    hpal=jet(100); %carte des couleurs: du bleu au rouge par defaut.
    pos1=get(gca,'Position');
    pos=[pos1(1) pos1(2)-.075 pos1(3) .013];
    colorbartype(pos,lev,1,clip_lev,hpal,0);
    title('degC')
  else
    colorbar('East')
  end;
  
   % sauvegarde de la figure
   printer='ps';
   print_printer=['-d', printer];
   fullfilename=['.' filesep nomfic num2str(limite_latN) 'N' num2str(limite_latS) 'S_' , num2str(ifigure) '.' printer];
   print(print_printer,fullfilename);
   clear fullfilename;
   clear printer;

   %subplot(4,2,2*i)
   ifigure=2;
   figure(ifigure);
   subplot(2,2,i)
   plot(PC(:,i),'x-')
   title(['composante principale correspondante (normalisée)'])
   hold on
   plot([1 ntemp],[0 0],'k')
   xlim([1 ntemp])
   xlabel('mois')
   %
   % sauvegarde de la figure
   printer='ps';
   print_printer=['-d', printer];
   fullfilename=['.' filesep nomfic num2str(limite_latN) 'N' num2str(limite_latS) 'S_' , num2str(ifigure) '.' printer];
   print(print_printer,fullfilename);
   clear fullfilename;
   clear printer;
end;

% Graphe pourcentages de variance expliqué
ifigure=3;
figure(ifigure);
bar(flipud(pctg_var))
xlim([0 10])
title('valeurs propres de la matrice de covariance')
ylabel('% de variance expliqué')
xlabel('indice du vecteur propre')

% sauvegarde de la figure
printer='ps';
print_printer=['-d', printer];
fullfilename=['.' filesep nomfic num2str(limite_latN) 'N' num2str(limite_latS) 'S_' , num2str(ifigure) '.' printer];
print(print_printer,fullfilename);
clear fullfilename;
clear printer;