source: trunk/src/affiche_EOF_PC_3D.m @ 325

function [EOF2D,PC]=affiche_EOF_PC_3D(mat_3D,date, lon,lat,num_eof,anomalies,mask,titre_fig,ydir,fac,aff_moy,aff_eof)

%
%+
%
% ===================
% affiche_EOF_PC_3D.m
% ===================
%
% .. function:: affiche_EOF_PC_3D(mat_3D,date, lon,lat,num_eof,anomalies,mask,titre_fig,ydir,fac)
%
% DESCRIPTION
% ===========
%
% matrice 3D de la forme (lon,lat,date)
%
% dates de la taille size(mat_3D,3)
%
% lon, lat, matrices de la taille [size(mat_3D,1),size(mat_3D,2)]
%
% num_eof: le numero des EOF a afficher. PAr defaut : les 3 premieres
%
% anomalies : booleen a mettre a true si on travaille sur des champs
% d'anomalies
%
% aff_moy : [min max step] pour les couleurs de matrice moyenne OU nombre
%           pour couleurs matlab OU [] (defaut) pour les rentrer a la main   
%
% aff_eof : [min max step] pour les couleurs des figures d'eof OU nombre 
%           pour couleurs matlab OU [] (defaut) pour les rentrer a la main      
%
% EOF faits sur le temps, PC 2D
%
% Affichage sur 1 seule page : mean mat_3D, EOF 1-3 et PC 1-3
% 
% EVOLUTIONS
% ==========
%
% $Id$
%
% $URL$
%
% - fplod 20110705T125628Z cratos.locean-ipsl.upmc.fr (Linux)
%
%   * minimal header
%
%-
  nb_eof=3;
sx=size(mat_3D,1);
sy=size(mat_3D,2);
load cmgco_traitdecote.mat Xgco Ygco

 if(nargin<3)
     lon=1:sy;
end
 if(nargin<4)
     lat=1:sx;
end
 if(nargin<5)
     num_eof=[1 2 3];
end
 if(nargin<6)
  anomalies=true;
end
 if(nargin<7)
     mask=ones(sx,sy);
%mask=ones(size(mat_3D,1),size(mat_3D,2));
end
 if(nargin<8)
     titre_fig='';
end
 if(nargin<9)
     ydir='normal';
     %sinon: 'reverse';
end
if(nargin<10)
  fac=[1 1 1];
end
if(nargin<11)
  aff_moy=[];
end
if(nargin<11)
  aff_eof=[];
end 

%preparation de la matrice
for(i=1:size(mat_3D,3))
     mat_3D(:,:,i)=mat_3D(:,:,i).*mask;
end
mat_2D=reshape(mat_3D,sx*sy,size(mat_3D,3));

mat_2D=mat_2D';%'

%elimination des points tq pour tout t mat_2D =NaN
mean_m2D=meanNaN(mat_2D,1);
mat_2D_ssnan=mat_2D(:,~isnan(mean_m2D));
clear mat_2D
%elimination des dates tq 1er ou dernier element de matrice = NaN
iidata=find(~isnan(mat_2D_ssnan(:,1)) & ~isnan(mat_2D_ssnan(:,end)));
mat_2D_ssnan=mat_2D_ssnan(iidata,:);
date=date(iidata);

%soustraction de la moyenne temporelle
size(mat_2D_ssnan);

for(i=1:size(mat_2D_ssnan,2))
     mat_2D_ssnan(:,i)= mat_2D_ssnan(:,i)-meanNaN(mat_2D_ssnan(:,i));
end


[U,S,V] = svd(mat_2D_ssnan',0); %'
clear mat_2D_ssnan
std_u=nstd(U,0,1);
diag_s=diag(S);
pc_var_exp=(diag_s.^2/sum(diag_s.^2))*100;
std_v=nstd(V,0,1);


for(i=1:nb_eof)
              PC(:,i)=  V(:,num_eof(i))/std_v(num_eof(i))*fac(i);     
              EOF(:,i)=U(:,num_eof(i))*std_v(num_eof(i))*diag_s(num_eof(i))*fac(i); 
end

%Plot des EOF
figure_a4p

hf=subplot(3,2,1);
if(anomalies)
  mat_3D_moy=squeeze(nstd(mat_3D,1,3));
else    
  mat_3D_moy=squeeze(meanNaN(mat_3D,3));
end     
if(min(size(lon))==1)
  imagesc(lon,lat,mat_3D_moy)
  set(gca,'ydir',ydir)
  hc=colorbar('WestOutside');
else
  if(length(aff_moy)==3)
    ech_c=aff_moy(1):aff_moy(3):aff_moy(2);
    num_cont=contourf_CB(lon,lat,mat_3D_moy,ech_c);
  else
    num_cont=contourf_CB(lon,lat,mat_3D_moy); 
  end
  set(gca,'ydir',ydir)
  CB_discrete(num_cont,'WestOutside');
  if(length(aff_moy)==0)  
    q =input('Changer l''affichage de la matrice moyenne (o/n) ? ', 's');
    if (q=='o')
      min_coul=input('valeur min : ');
      max_coul=input('valeur max : ');
      stp_coul=input('step : ');
      hf=subplot(3,2,1);              
      num_cont=contourf_CB(lon,lat,mat_3D_moy,[min_coul:stp_coul: max_coul]);
      set(gca,'ydir',ydir)
      CB_discrete([min_coul-stp_coul:stp_coul: max_coul+stp_coul],'WestOutside');
      
    end
  end
end

hold on
plot(Xgco,Ygco,'k.','MarkerSize',1)

if(anomalies)
  title('std of the field')
else    
  title('mean field')
end
ajoute_date_fig(gcf,3);



for(i=1:3)
  subplot(3,2,i+1)
  EOFnan=NaN(sx*sy,1);
  EOFnan(~isnan(mean_m2D))=EOF(:,i);
  EOF2D(:,:,i)=reshape(EOFnan,sx,sy);
  
  if(min(size(lon))==1)
    
    imagesc(lon,lat,EOF2D(:,:,i))
    set(gca,'ydir','normal')
    if(rem(i,2)==1)
      hc=colorbar('EastOutside');
    else
      hc=colorbar('WestOutside');
    end
  else
    if(length(aff_eof)==3)
      ech_c=aff_eof(1):aff_eof(3):aff_eof(2);
      num_cont=contourf_CB(lon,lat,squeeze(EOF2D(:,:,i)),ech_c);
      else  
      num_cont=contourf_CB(lon,lat,squeeze(EOF2D(:,:,i)));
    end
    set(gca,'ydir',ydir)
    if(rem(i,2)==1)
      CB_discrete(num_cont,'EastOutside');
    else
      CB_discrete(num_cont,'WestOutside');
    end
  
  end
  
  eval ([' title(''EOF ',num2str(num_eof(i)),': var exp = ', num2str(pc_var_exp(num_eof(i)),2),'%'')'])
  hold on
  plot(Xgco,Ygco,'k.','MarkerSize',1)
  
end
if(length(aff_eof)==0)
q =input('Changer les parametres de l''affichage (o/n) ? ', 's');
if (q=='o')
  min_coul=input('valeur min : ');
  max_coul=input('valeur max : ');
  if(min(size(lon))>1)
    stp_coul=input('step : ');
  end
  for(i=1:3)
    subplot(3,2,i+1)
    if(min(size(lon))==1)
      
      imagesc(lon,lat,EOF2D(:,:,i),[min_coul max_coul])
      set(gca,'ydir',ydir)
      caxis([min_coul max_coul])
      if(rem(i,2)==1)
        hc=colorbar('EastOutside');
      else
        hc=colorbar('WestOutside');
      end
    else
      num_cont=contourf_CB(lon,lat,squeeze(EOF2D(:,:,i)),[min_coul:stp_coul: max_coul]);
      set(gca,'ydir',ydir)
      if(rem(i,2)==1)
        CB_discrete([min_coul-stp_coul:stp_coul: max_coul+stp_coul],'EastOutside');
      else
        CB_discrete([min_coul-stp_coul:stp_coul: max_coul+stp_coul],'WestOutside');
      end
    end
    eval ([' title(''EOF ',num2str(num_eof(i)),': var exp = ', num2str(pc_var_exp(num_eof(i)),2),'%'')'])
    hold on
    plot(Xgco,Ygco,'k.','MarkerSize',1)
    
  end     
  
end
end
              subplot(3,2,5:6)  
              %Plot des PC
              
              plot(date, PC(:,1),date, PC(:,2),date, PC(:,3))
              hold on
              plot(date,zeros(size(date)),'k')
              if(date(end)-date(1)>1000)
              datetick('x','yyyy') 
              else
              datetick('x','mmmyy')
              end
              axis([date(1) date(end) -4 4]);
              
              legend(['PC',num2str(num_eof(1))],['PC',num2str(num_eof(2))],['PC',num2str(num_eof(3))])
              title(titre_fig)