source: trunk/src/affiche_EOF_PC_3Dsscycle.m @ 327

function [EOF2D,PC,varargout]=affiche_eofpc_3Drev(mat_3D,date, lon,lat,titre_fig,anomalies,mask,ydir,fac)
%
%+
%
% ==========================
% affiche_EOF_PC_3Dsscycle.m
% ==========================
%
% .. function:: affiche_eofpc_3Drev(mat_3D,date, lon,lat,titre_fig,anomalies,mask,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)]
%
% anomalies : booleen a mettre a true si on travaille sur des champs d''anomalies
% EOF faits sur le temps, PC 2D
%
% Affichage sur 1 seule page : mean mat_3D, EOF 1-3 et PC 1-3
%
% Output: EOF et PC; Output optionnel: pc_var_exp si spécifié  
%
% EVOLUTIONS
% ==========
%
% $Id: affiche_EOF_PC_3Dsscycle.m 325 2011-08-04 15:30:01Z pinsard $
%
% $URL$
%
% - fplod 20110705T125628Z cratos.locean-ipsl.upmc.fr (Linux)
%
%   * minimal header
%
%-
%

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)
     titre_fig='';
end
 if(nargin<6)
  anomalies=true;
end
 if(nargin<7)
     mask=ones(sx,sy);
end
 if(nargin<8)
     ydir='normal';
end
 if(nargin<9)
     fac=[1 1 1];
 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: 1)
              PC(:,i)=  V(:,i)/std_v(i)*fac(i);     
              EOF(:,i)=U(:,i)*std_v(i)*diag_s(i)*fac(i); 
end

%test d elimination du cycle saisonnier (eof 1)
 
 %if(nargin<5)
 mat_2D_ssnan=mat_2D_ssnan-PC(:,1)*(EOF(:,1))';%'


[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: 3)
              PC(:,i)=  V(:,i)/std_v(i)*fac(i);     
              EOF(:,i)=U(:,i)*std_v(i)*diag_s(i)*fac(i);
end
%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
              num_cont=contourf_CB(lon,lat,mat_3D_moy);
              set(gca,'ydir',ydir)
              CB_discrete(num_cont,'WestOutside');

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

              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
              num_cont=contourf_CB(lon,lat,squeeze(EOF2D(:,:,i)));
              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(i),': var exp = ', num2str(pc_var_exp(i),2),'%'')'])
              hold on
              plot(Xgco,Ygco,'k.','MarkerSize',1)
              end

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(i),': var exp = ', num2str(pc_var_exp(i),2),'%'')'])
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')
              datetick('x','yyyy') 
              axis([date(1) date(end) -4 4]);
              
              legend('PC1','PC2','PC3')
              title(titre_fig)
              
if(nargout==3)
  pc_var_exp=pc_var_exp(1:3);
  varargout={pc_var_exp};
end