Index: /issm/trunk-jpl/externalpackages/pcatool/install.sh =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install.sh (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install.sh (revision 21589) @@ -0,0 +1,12 @@ +#!/bin/bash + +#Some cleanup +rm -rf install pcatool +mkdir install + +#Download from ISSM server +$ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/pcatool.tar.gz' 'pcatool.tar.gz' + +#Untar into install +cd install +tar -zxvf ../pcatool.tar.gz Index: /issm/trunk-jpl/externalpackages/pcatool/install/Contents.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/Contents.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/Contents.m (revision 21589) @@ -0,0 +1,32 @@ +% PRINCIPAL COMPONENTS ANALYSIS TOOLBOX +% Last Modified: Dec 2007 +% +% This software is provided "as is" without express +% or implied warranty... +% +% +% EOF and SVD +% caleof - Compute EOFs +% calEeof - Compute Extended EOFs +% calCeof - Compute Complex EOFs +% calsvd - Compute SVDs +% confexpvar - Compute statistical confidence of explained variance +% +% VISUALISATION +% dispEOF - Display EOFs +% dispEEOF - Display EEOFs as time sequence map +% dispEEOF_v2 - Display EEOFs as time sequence map with +% the PC time series +% dispCEOF - Display CEOFs +% dispSVD - Display SVDs +% plotm - Plot a 2D geo field +% +% Matrix manipulation +% map2mat - Transform C(TIME,LON,LAT) in D(TIME,PT) +% under a geographic points selection +% mat2map - Invert transform done by map2mat +% +% +%================================================================ +% Author: Guillaume MAZE - LPO/LMD - March 2004 +% gmaze@univ-brest.fr Index: /issm/trunk-jpl/externalpackages/pcatool/install/calCeof.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/calCeof.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/calCeof.m (revision 21589) @@ -0,0 +1,11 @@ +% [e,pc,expvar] = calCeof(M,N,METHOD) +% Compute Complex EOF + + function [e,pc,expvar] = calCeof(M,N,METHOD) + + +% Compute Hilbert transform Mh and create M as: M = M + i*Mh +M = hilbert(M); + +% Compute EOFs + [e,pc,expvar] = caleof(M,N,METHOD); Index: /issm/trunk-jpl/externalpackages/pcatool/install/calEeof.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/calEeof.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/calEeof.m (revision 21589) @@ -0,0 +1,128 @@ +% [EOFs,PC,EXPVAR] = calEeof(M,N,METHOD,NLAG,DT) Compute EEOF +% +% => Compute N Extended EOF of a M(TIME*MAP) matrix. +% METHOD is the EOF method computing (see caleof.m) +% NLAG is the number of map you want +% DT is the increment between each map (not the time step +% of your data) +% +% The function can filter datas with a Lanczos band-pass +% filter. Edit the calEeof.m file and change option in it. +% +% See also: CALEOF +% +% Ref: Weare and Nasstrom 1982 +%================================================================ + +% Guillaume MAZE - LPO/LMD - March 2004 +% gmaze@univ-brest.fr + + function varargout = calEeof(M,N,METHOD,NLAG,DT) + +% --------------------------------- +% Divers +tic +LAG = NLAG*DT; +foll = 1; % Turn this value to 0 if you don t want commentary along computing + +% --------------------------------- +% LAG must be >0, LAG=1 correspond to initial position +if LAG<=0 + LAG = 1; +end + +% --------------------------------- +% Usely for us, M is a TIME*MAP matrix so: +% We get dimensions + [n p]=size(M); +% and put M in the correct form for computing + M = M'; + +% --------------------------------- +% Eventualy time filtering of data +if 0==1 + if (foll~=0),disp('==> Time filtering...');end + dt = 1; % Real time step of your data (diff from DT) + Fc = 1/dt/8; % This the low cut off frequency + Fc2 = 1/dt/2; % and this the high cut off frequency + SIGNAL = M(1,:); + nj = fix(length(SIGNAL)/10); % Nb of points of the window + + for ipt = 1 : p + SIGNAL = M(ipt,:)'; + SIGNALF = lanczos(SIGNAL,Fc2,nj); + SIGNALF = SIGNALF - lanczos(SIGNALF,Fc,nj); + Y(:,ipt) = SIGNALF; + if mod(ipt,10)==0 % We display a plot of filtered data + clf; + plot((0:n-1),SIGNAL,'k',(0:n-1),SIGNALF,'r'); + drawnow; + end + end + M = Y'; +end, clear Fc Fc2 SIGNAL nj SIGNALF Y ipt, close + + +% --------------------------------- +% This is the matrix where we concat all submatrices from M. +% Size of each submatrix is : NLAG*p rows vs n-LAG+1 colums +F = zeros(NLAG*p,n-LAG+1); + +% --------------------------------- +% Let's go for F: +if (foll~=0),disp('==> Forming concatenated matrix...');end +if DT==1 +% CLASSIC CASE +for ilag = 1 : LAG + % Extract submatrix ilag from M + Msub = M(:,ilag:n-LAG+ilag); + % Concat matrix + F( (ilag-1)*p+1 : ilag*p , : ) = Msub; +end +else +% DT>1 +it=0; +for ilag = 1 : DT : LAG + % Extract submatrix ilag from M + Msub = M( : , ilag : n-LAG+ilag); + % Concat matrix + it = it + 1; + F( (it-1)*p+1 : it*p , : ) = Msub; + % imagesc(Msub);pause % for debugging +end +end + +% --------------------------------- +% Compute EOFs by normal way +% (Don't forget that caleof is taking a TIME*MAP argument matrix, +% that's why we have F' !) +if (foll~=0),disp('==> Computing EOFs ...');end +[e,pc,expvar] = caleof(F',N,METHOD); + +% --------------------------------- +% e is the matrix with EOFs inside. We should extract each map +% for different lags from it. +% e is NEOF*(LAG*MAP) and we construct 3D matrix CHP as NEOF*LAG*MAP +if (foll~=0),disp('==> Reshaping vectors ...');end +for ilag = 1 : NLAG + E = e( : , (ilag-1)*p+1 : ilag*p ); + CHP(:,ilag,:) = E; +end + + +% --------------------------------- +% OUTPUT VARIABLES +switch nargout + case 1 + varargout(1) = {CHP} ; + case 2 + varargout(1) = {CHP} ; + varargout(2) = {pc}; + case 3 + varargout(1) = {CHP} ; + varargout(2) = {pc}; + varargout(3) = {expvar}; +end + +% --------------------------------- +if (foll~=0),toc,disp('==> That''s all folks !');end Index: /issm/trunk-jpl/externalpackages/pcatool/install/caleof.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/caleof.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/caleof.m (revision 21589) @@ -0,0 +1,234 @@ +% [EOFs,PC,EXPVAR] = CALEOF(M,N,METHOD) Compute EOF +% +% => Compute the Nth first EOFs of matrix M(TIME,MAP). +% EOFs is a matrix of the form EOFs(N,MAP), PC is the principal +% components matrix ie it has the form PC(N,TIME) and EXPVAR is +% the fraction of total variance "explained" by each EOF ie it has +% the form EXPVAR(N). +% Differents method can be used: +% 1 - The "classic" one, ie we compute eigenvectors of the +% temporal covariance matrix with the eig Matlab function. +% 2 - A faster "classic" one, same as method 1 but we use the +% eigs Matlab function. +% 3 - We compute eigenvectors by the singular value decomposition, +% by used of the svd Matlab function. +% 4 - Same as method 3 but faster by used of the svds Matlab function +% +% See also EIG, EIGS, SVD, SVDS +% +% Ref: L. Hartmann: "Objective Analysis" 2002 +% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - +% Analyses of climatic Data" 1997 +%================================================================ + +% Guillaume MAZE - LPO/LMD - March 2004 +% Revised July 2006 +% gmaze@univ-brest.fr + + +function [e,pc,expvar,L] = caleof(M,N,method); + +% Get dimensions +[n p]=size(M); + +% Temporal covariance is p*p matrix, that why max EOF computable is p, +% so we perform a test on parameter N: +if(N>p) + disp('Warning: N is larger than possible so it''s modified to perform') + disp(' EOFs computing...'); + N = p; +end + + +% Eventualy time filtering of data +if 0==1 + disp('====> Time filtering...') + Fc = 1/20; Fc2 = 1/1; + Fc = 1/7 ; Fc2 = 1/3; + SIGNAL = M(:,1); + nj = fix(length(SIGNAL)/10); % Nombre de points du filtre + for ipt = 1 : p + SIGNAL = M(:,ipt); + SIGNALF = lanczos(SIGNAL,Fc2,nj); + SIGNALF = SIGNALF - lanczos(SIGNALF,Fc,nj); + Y(:,ipt) = SIGNALF; + end + M = Y; +end + + +disp('====> Let''go for EOFs and pc computing...') +switch method + case 1 % CLASSIC METHOD +%================================================================ +% Transform the data matrix in the correct form (map*time) for eig +M = M'; + +% Remove the time mean (ie the mean of each rows of M) +% Rq: detrend remove the mean of columns ie we take M'. +F = detrend(M','constant')'; + +% Covariance Matrix (inner product over space = covariance in time) +R = F * F'; + +% Eigenanalysis of the covariance matrix R +[E,L] = eig(R); + +% Get PC by projecting eigenvectors on original data +Z = E'*F; + +% Make them clear for output +for iN=1:N + e(iN,:) = squeeze( E(:,p-(iN-1)) )'; + pc(iN,:) = squeeze( Z(p-(iN-1),:) ); +end + +% Amount of explained variance (at 0.1%) +dsum = diag(L)./trace(L); +for iN=1:N + expvar(iN)=fix((dsum(p-(iN-1))*100/sum(dsum))*10)/10; +end + +% Plots Original field and reconstructed one +if 0==1 +figure; +subplot(1,2,1);imagesc(abs(M));title('ORIGINAL');cx=caxis; +%subplot(1,2,2);imagesc((E*Z));title('RECONSTRUCTED') +subplot(1,2,2);imagesc(abs(e'*pc));title('RECONSTRUCTED');caxis(cx); +end + + case 2 % RAPID CLASSIC METHOD +%================================================================ +% Remove the time mean of each column +F = detrend(M,'constant'); + +% Covariance Matrix +if n >= p + R = F' * F; +else + R = F * F'; +end + +% Eigen analysis of the square covariance matrix +[E,L] = eigs(R,N); +if n < p + E = F' * E; + sq = [sqrt(diag(L))+eps]'; + sq = sq(ones(1,p),:); + E = E ./ sq; +end + +% Get PC by projecting eigenvectors on original data +if n >= p + Z = (F*E)'; +else + Z = E'*F'; +end + + +% Make them clear for output +for iN=1:N + e(iN,:) = squeeze( E(:,iN) )'; + pc(iN,:) = squeeze( Z(iN,:) ); +end + +% Amount of variance explained a 0.1 pres et en % +dsum=diag(L)./trace(L); +for iN=1:N + expvar(iN)=fix((dsum(iN)*100/sum(dsum))*10)/10; +end + + + case 3 % SVD METHOD +%================================================================ +% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - +% Analyses of climatic Data" 1997 => p18 + +% Assume that M is (time*map) matrix +[n p]=size(M); + +% Remove the mean of each column (ie the time mean in each station records) +F=detrend(M,'constant'); + +% Form the covariance matrix: +R = F'*F; + +% Find eigenvectors and singular values +[C,L,CC] = svd(R); +% Eigenvectors are in CC and the squared diagonal values of L +% are the eigenvalues of the temporal covariance matrix R=F'*F + +% find the PC corresponding to eigenvalue +PC = F*CC; + +% Make them clear for output +for iN=1:N + e(iN,:) = squeeze( CC(:,iN) )'; + pc(iN,:) = squeeze( PC(:,iN) )'; +end + +if 0 +figure; +subplot(1,2,1);imagesc(F);title('ORIGINAL');cx=caxis; +subplot(1,2,2);imagesc(C*L*CC');title('RECONSTRUCTED');caxis(cx); +end + +% Amount of variance explained at 0.1% +dsum=diag(L)./trace(L); +if length(dsum) p18 + +% Assume that M is (time*map) matrix +[n p]=size(M); + +% Remove the mean of each column (ie the time mean in each station records) +F=detrend(M,'constant'); + +% Form the covariance matrix: +R = F' * F; + +% Find eigenvectors and singular values +[C,L,CC,flag] = svds(R,N); +% Eigenvectors are in CC and the squared diagonal values of L +% are the eigenvalues of the temporal covariance matrix R=F'*F +% (Sometimes, CC stops for nul eigenvector, then we need to fill to reach N) +if size(CC,2) Finished !') Index: /issm/trunk-jpl/externalpackages/pcatool/install/calsvd.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/calsvd.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/calsvd.m (revision 21589) @@ -0,0 +1,54 @@ +% [SVD1,SVD2,PC1,PC2,EXPVAR,Lambda] = CALSVD(A,B,N) Compute SVDs +% +% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - +% Analyses of climatic Data" 1997 +%================================================================ +% +% Guillaume MAZE - LPO/LMD - March 2004 +% gmaze@univ-brest.fr + + + function [e1,e2,pc1,pc2,expvar,Lambda] = calsvd(A,B,N); + + +disp('====> Let''go for SVDs and pc computing...') + +% Assume that A is (time*map) matrix +[n p]=size(A); + +% Remove the mean of each column (ie the time mean in each station records) +S=detrend(A,'constant'); +P=detrend(B,'constant'); + +% Form the covariance matrix: +C = S'*P; + +% Find eigenvectors and singular values +[U,Lambda,V] = svds(C,N); + +% PC +a = S*U; +b = P*V; + +% Make them clear for output +for iN=1:N + e1(iN,:) = squeeze( U(:,iN) )'; + pc1(iN,:) = squeeze( a(:,iN) )'; + e2(iN,:) = squeeze( V(:,iN) )'; + pc2(iN,:) = squeeze( b(:,iN) )'; +end + +% Amount of variance explained a 0.1 pres et en % +L2=Lambda.^2; +dsum=diag(L2)/trace(L2); + +for iN=1:N + expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10; +end + + +disp('====> Finished !') + + + + Index: /issm/trunk-jpl/externalpackages/pcatool/install/calsvd2.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/calsvd2.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/calsvd2.m (revision 21589) @@ -0,0 +1,49 @@ +% [SVD1,SVD2,PC1,PC2,EXPVAR,Lambda] = CALSVD2(A,B,N) Compute SVDs +% +% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - +% Analyses of climatic Data" 1997 +%================================================================ +% +% Guillaume MAZE - LPO/LMD - March 2004 +% gmaze@univ-brest.fr + + + function [e1,e2,pc1,pc2,expvar,Lambda,dsumCF] = calsvd2(A,B,N); + + +%================================================================ +% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD - +% Analyses of climatic Data" 1997 => p18 + +% Assume that A is (time*map) matrix +[n p]=size(A); + +% Remove the mean of each column (ie the time mean in each station records) +S=detrend(A,'constant'); +P=detrend(B,'constant'); + +% Form the covariance matrix: +C=S'*P; + +% Find eigenvectors and singular values +[U,Lambda,V] = svds(C,N); + +% PC +a=S*U; +b=P*V; + +% Make them clear for output +for iN=1:N + e1(iN,:) = squeeze( U(:,iN) )'; + pc1(iN,:) = squeeze( a(:,iN) )'; + e2(iN,:) = squeeze( V(:,iN) )'; + pc2(iN,:) = squeeze( b(:,iN) )'; +end + +% Amount of variance explained a 0.1 pres et en % +L2=Lambda.^2; +dsum=diag(L2)/trace(L2); +for iN=1:N + expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10; +end + Index: /issm/trunk-jpl/externalpackages/pcatool/install/confexpvar.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/confexpvar.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/confexpvar.m (revision 21589) @@ -0,0 +1,67 @@ +% [CONF,CONFLEVEL] = CONFEXPVAR(A,B,nsamples,N) +% +% A(Ntime,MmapA), B(Ntime,MmapB) are time-synchronised +% and time-constant sampled data sets. +% MmapA and MmapB are not required to be identic. +% nsamples is what it is (given at least a 100 factor +% will be nice to provide a trustable result). +% N is the max number of eigen vectors to be computed. +% +% CONFEXPVAR Compute the confidence limits for significant +% eigenvalues in a SVD decomposition of data sets A and B. +% This generates nsamples random-in-time data set of B +% and finds percent variance of N SVD modes computed +% with data set A. +% This is exactly the monte-carlo method but based on +% a random-in-time only data set. This considerably +% decrease the cputime. +% +% 2005/11. Guillaume Maze +% gmaze@univ-brest.fr + +function [conf,lev,dsumCF] = confexpvar(A,B,nsamples,N) + +% We try to find the p confidence levels of expvar. +% Keeping as it is the first field, we determine SVDs with the +% second field arbitrary sorted in time. +% 1st field : A +% 2nd field : B + +n = size(A,1); +if size(B,1)~=n + disp('Error, A and B required to be of the same time length !'); + disp('See help confexpvar'); + stop +end + +timer = cputime; +disp('Compute confidence levels of explained variance...'); +for iR=1:nsamples + if iR==fix(nsamples/2) + disp(strcat(num2str(iR),'/',num2str(nsamples))); + disp(['Please, still wait a little time ... around '... + num2str(cputime-timer) ' cpu seconds']); + end + + % Generate a random new time index: + t=randperm(n); + % and sort the 2nd field with it: + Br=B(t,:); + + % Now we compute SVD and expvar: + C=A'*Br; + [U,Lambda,V] = svds(C,N); + L2 = Lambda.^2; dsum = diag(L2)/trace(L2); + L = Lambda; dsum = diag(L)/trace(L); + dsumCF(iR,:) = dsum'; + +end %for iR +disp(['OK, this took ' num2str(cputime-timer) ' cpu seconds']) + +% CF(nsamples,N) contains expvar of the nsamples realisations: +CF=sort(dsumCF,1); + +% Now we just take the decadal significant levels: +lev=fix([nsamples/10:nsamples/10:nsamples]); +conf = CF(lev,:) ; +lev=(lev./max(lev))*100; Index: /issm/trunk-jpl/externalpackages/pcatool/install/dispCEOF.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/dispCEOF.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/dispCEOF.m (revision 21589) @@ -0,0 +1,65 @@ +% dispCEOF(CEOF,EXPVAR,PHI,MOD) Display few CEOFs. +% +% => DISPLAY FEW CEOFs. +% CEOF contains all the CEOFs as nCEOF*X*Y. +% EXPVAR is a matrix with the explained variance of each +% CEOFs in %. This is just for title. +% PHI is phase when ploting maps. +% MOD contains explicitly the CEOFs to display. +% +% Rq: it's using the plotm function +% +%================================================================ + +% Guillaume MAZE - LPO/LMD - June 2004 +% gmaze@univ-brest.fr + +function [] = dispCEOF(CEOF,EXPVAR,PHI,MOD) + +% Number of EEOF to display +NMOD = length(MOD); +NANG = length(PHI); + +% Open figure and first guest +figur;clf;hold on +typg=8; +width = .9/NMOD; +height = .9/NANG; +dleft = (.95-width*NMOD)/4; +left = dleft; +load mapanom2 + +% Let'go : + +% We select a CEOF: +for imod = 1 : NMOD + mod = MOD(imod); + chp=squeeze(CEOF(mod,:,:)); + Rchp=real(chp); + Ichp=-imag(chp); + + % + for ang=1:NANG + C = cos(PHI(ang))*Rchp + sin(PHI(ang))*Ichp; + m(ang)=abs(xtrm(C)); + end + m=max(m); + + % and plot maps + for ang=1:NANG + C = cos(PHI(ang))*Rchp + sin(PHI(ang))*Ichp; + bottom = .95-height*ang; + subplot('Position',[left bottom width height]); + plotm(C,typg); + m_text(10,-70,strcat('\phi=',num2str(PHI(ang)*180/pi),'^o')); +% ylabel(strcat('\phi=',num2str(PHI(ang)*180/pi),'^o'),'Rotation',0); + caxis([-m m]); + colormap(mapanom); + if(ang==1) + titre=strcat('CEOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)'); + title(titre); + end; + end + left = left + width + dleft; + +end Index: /issm/trunk-jpl/externalpackages/pcatool/install/dispEEOF.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/dispEEOF.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/dispEEOF.m (revision 21589) @@ -0,0 +1,57 @@ +% dispEEOF(CHP,EXPVAR,DT,NLAG,MOD) Display few EEOFs. +% +% => DISPLAY FEW EEOFs. +% CHP contains all the EEOFs as EOF*LAG*X*Y. +% EXPVAR is a matrix with the explained variance of each +% EEOFs in %. This is just for title. +% DT is time step between maps. +% NLAG is the number of LAG to display. +% MOD contains explicitly the EEOFs to display. +% +% Rq: it's using the plotm function +% +%================================================================ + +% Guillaume MAZE - LPO/LMD - March 2004 +% gmaze@univ-brest.fr + +function [] = dispEEOF(CHP,EXPVAR,DT,NLAG,MOD) + +% Number of EEOF to display +NMOD = length(MOD); + +% Open figure and first guest +figur; +clf;hold on +typg=8; +width = .9/NMOD; +height = .9/NLAG; +dleft = (.95-width*NMOD)/4; +left = dleft; +load mapanom + +% Let'go : + +% We choose an EEOF +for imod = 1 : NMOD + mod = MOD(imod); + + % and plot maps + for lag=1:NLAG + bottom = .95-height*lag; + subplot('Position',[left bottom width height]); + C = squeeze(real(CHP(mod,lag,:,:))); +% C = C./xtrm(C); % Eventually normalise field + plotm(C,typg); + if(lag==1),caxis([-abs(xtrm(C)) abs(xtrm(C))]);cx=caxis;end; + caxis(cx); + colormap(mapanom); + if(lag==1) + titre=strcat('EOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)'); + titre=strcat(titre,'; DT=',num2str(DT),' '); + title(titre); + end; + end + left = left + width + dleft; + +end Index: /issm/trunk-jpl/externalpackages/pcatool/install/dispEEOF_v2.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/dispEEOF_v2.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/dispEEOF_v2.m (revision 21589) @@ -0,0 +1,76 @@ +% dispEEOF_v2(CHP,PC,EXPVAR,DT,NLAG,MOD) Display few EEOFs with their +% PC time evolution +% +% => DISPLAY FEW EEOFs with their PC time evolution. +% CHP contains all the EEOFs as EOF*LAG*X*Y. +% PC contains all the PCs. +% EXPVAR is a matrix with the explained variance of each +% EEOFs in %. This is just for title. +% DT is time step between maps. +% NLAG is the number of LAG to display. +% MOD contains explicitly the EEOFs to display. +% +% Rq: it's using the plotm function +% See also: dispEEOF +% +%================================================================ + +% Guillaume MAZE - LPO/LMD - June 2004 +% gmaze@univ-brest.fr + +function [] = dispEEOF(CHP,PC,EXPVAR,DT,NLAG,MOD) + +% Number of EEOF to display +NMOD = length(MOD); + +% Open figure and first guest +figur; +clf;hold on +typg=8; +width = .9/NMOD; +height = .9/(NLAG+1); +dleft = (.95-width*NMOD)/4; +left = dleft; +load mapanom + +% Let'go : + +% We choose an EEOF +for imod = 1 : NMOD + mod = MOD(imod); + + % and plot maps + for iplot=1:NLAG+1 + bottom = .95-height*iplot; + subplot('Position',[left bottom width height]); + if iplot==1 % Plot the pc's time serie + plot(PC(imod,:)); + grid on;box on; + set(gca,'XAxisLocation','top');axis tight; + titre=strcat('EEOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)'); + titre=strcat(titre,'; DT=',num2str(DT),' '); + title(titre); + if (imod==1),ax1=get(gca,'YLim');else,set(gca,'YLim',ax1);end + else % Plot the lag maps + lag=iplot-1; + C = squeeze(real(CHP(mod,lag,:,:))); +% C = C./xtrm(C); % Eventually normalise field + plotm(C,typg); + if(lag==1),caxis([-abs(xtrm(C)) abs(xtrm(C))]);cx=caxis;end; + caxis(cx); + colormap(mapanom); + %ylabel(strcat('{\bf',num2str((lag-1)*DT),'}')) + m_text(10,-70,strcat('{\bf',num2str((lag-1)*DT),'}')); + end %if + + end + left = left + width + dleft; + + % Adjust plot width via the plot Aspect ratio + ch=get(gcf,'Children'); + posiPC=get(ch(iplot),'PlotBoxAspectRatio'); + posiMAP=get(ch(iplot-1),'PlotBoxAspectRatio'); + set(ch(iplot),'PlotBoxAspectRatio',posiMAP); + +end + Index: /issm/trunk-jpl/externalpackages/pcatool/install/dispEOF.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/dispEOF.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/dispEOF.m (revision 21589) @@ -0,0 +1,130 @@ +% dispEOF(typg,CHP,N,expvar) Display EOFs +% +% => Display EOFs +% typg is plot type you want (see plotm.m). +% CHP contains EOFs and is of the form: CHP(EOF,X,Y) +% N is number of EOF to display +% expvar contained the explained variance of EOFs +% +% Rq: it''s using the plotm function +% +%================================================================ + +% Guillaume MAZE - LPO/LMD - March 2004 +% gmaze@univ-brest.fr + +function [] = dispEOF(typg,CHP,N,expvar) + +load mapanom2 +global ffile +% -------------------------------------------------- +% AFFICHAGE +% -------------------------------------------------- +if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6) + + CBAROR='horiz'; +if N>12,NN=12;else,NN=N;end +switch NN + case 1,tblmagique=[1]; iw=1;jw=1; + case 2,tblmagique=[1 2]; iw=2;jw=1; + case 3,tblmagique=[1 2 3]; iw=3;jw=1; + case 4,tblmagique=[1 2 3 4]; iw=2;jw=2; + case 5,tblmagique=[1 2 3 4 5]; iw=3;jw=2; + case 6,tblmagique=[1 2 3 4 5 6]; iw=3;jw=2; + case 7,tblmagique=[1 2 3 4 5 6 7]; iw=4;jw=2; + case 8,tblmagique=[1 2 3 4 5 6 7 8]; iw=4;jw=2; + case 9,tblmagique=[1 2 3 4 5 6 7 8 9]; iw=3;jw=3; + case 10,tblmagique=[1 2 3 4 5 6 7 8 9 10]; iw=4;jw=3; + case 11,tblmagique=[1 2 3 4 5 6 7 8 9 10 11]; iw=4;jw=3; + case 12,tblmagique=[1 2 3 4 5 6 7 8 9 10 11 12];iw=4;jw=3; +end +iw=N;jw=1; + +else + + CBAROR='vert'; +if N>12,NN=12;else,NN=N;end +switch NN + case 1,tblmagique=[1]; iw=1;jw=1; + case 2,tblmagique=[1 2]; iw=2;jw=1; + case 3,tblmagique=[1 2 3]; iw=3;jw=1; + case 4,tblmagique=[1 2 3 4]; iw=2;jw=2; + case 5,tblmagique=[1 2 3 4 5]; iw=3;jw=2; + case 6,tblmagique=[1 2 3 4 5 6]; iw=3;jw=2; + case 7,tblmagique=[1 2 3 4 5 6 7]; iw=4;jw=2; + case 8,tblmagique=[1 2 3 4 5 6 7 8]; iw=4;jw=2; + case 9,tblmagique=[1 2 3 4 5 6 7 8 9]; iw=3;jw=3; + case 10,tblmagique=[1 2 3 4 5 6 7 8 9 10]; iw=4;jw=3; + case 11,tblmagique=[1 2 3 4 5 6 7 8 9 10 11]; iw=4;jw=3; + case 12,tblmagique=[1 2 3 4 5 6 7 8 9 10 11 12];iw=4;jw=3; +end + + + +end %if + + + +f1=figur;clf; hold on + +for iN=1:NN + subplot(iw,jw,tblmagique(iN)); hold on +C = squeeze(real(CHP(iN,:,:))); +%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]); +plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]); +if(iN==1);cx=caxis;end;caxis(cx); +colormap(mycolormap(mapanom,21)); + +%c=colorbar(CBAROR); +titre1=strcat('EOF',num2str(iN)); +titre2=strcat(' (',num2str(expvar(iN)),'%)'); +title(strcat(titre1,titre2)); + +end %for iN + +%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2))); +titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%'); +%titre3=strcat('(Analysed file: ',ffile,')'); +%titre=strvcat(titre2,titre3); +suptitle(titre2); + +set(f1,'Position',[378 39 313 647]); % Laptop screen +%set(f1,'Position',[369 55 316 899]); % Desktop screen +%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART')); +set(f1,'Name','<> REAL PART <>'); + + +% ------------------------------------------- +if isreal(CHP)==0 + +f2=figur;clf; hold on + +for iN=1:NN + subplot(iw,jw,tblmagique(iN)); hold on + +C = squeeze(imag(CHP(iN,:,:))); +plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]); +%if(iN==1);cx=caxis;end; +caxis(cx); +colormap(mycolormap(mapanom,21)); + +c=colorbar(CBAROR); +titre1=strcat('EOF',num2str(iN)); +titre2=strcat(' (',num2str(expvar(iN)),'%)'); +title(strcat(titre1,titre2)); + +end %for iNN + +%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2))); +titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%'); +%titre3=strcat('(Analysed file: ',ffile,')'); +%titre=strvcat(titre2,titre3); +suptitle(titre2); + +set(f2,'Position',[378 39 313 647]); % Laptop screen +%set(f2,'Position',[369 55 316 899]); % Desktop screen +set(f2,'Name','<> IMAGINARY PART <>'); + +end + + Index: /issm/trunk-jpl/externalpackages/pcatool/install/dispSVD.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/dispSVD.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/dispSVD.m (revision 21589) @@ -0,0 +1,86 @@ +% dispSVD(typg,CHP1,CHP2,N,expvar) Display SVDs +% +% => Display SVDs +% typg is plot type you want (see plotm.m). +% CHP contains SVDs and is of the form: CHP(iSVD,X,Y) +% N is number of SVD to display +% expvar contained the explained variance of SVDs +% +% Rq: it''s using the plotm function +% +%================================================================ + +% Guillaume MAZE - LPO/LMD - December 2004 +% gmaze@univ-brest.fr + +function [] = dispSVD(typg,CHP1,CHP2,N,expvar) + +load mapanom2 +global ffile1 ffile2 + +if length(N)>1 + mod=N; +else + mod=[1:N]; +end +Nx=length(mod); +% -------------------------------------------------- +% AFFICHAGE +% -------------------------------------------------- +if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6) + CBAROR='horiz'; + iw=Nx;jw=2; +else + CBAROR='vert'; + iw=Nx;jw=2; +end + + +f1=figur;clf; hold on +for iN=1:Nx + subplot(iw,jw,2*iN-1) +C = squeeze(real(CHP1(mod(iN),:,:))); +%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]); +plotm(C,typg); +caxis([-abs(xtrm(C)) abs(xtrm(C))]); +%if(iN==1);cx=caxis;end;caxis(cx); +colormap(mycolormap(mapanom,21)); +c=colorbar(CBAROR); +titre1=strcat('SVD',num2str(mod(iN)),':',blanks(1),ffile1); +titre2=strcat(' (',num2str(expvar(mod(iN))),'%)'); +tt=title(strvcat(strcat(titre1,titre2),char(1))); +set(tt,'fontweight','bold') + + subplot(iw,jw,2*iN) +C = squeeze(real(CHP2(mod(iN),:,:))); +%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]); +plotm(C,typg); +caxis([-abs(xtrm(C)) abs(xtrm(C))]); +%if(iN==1);cx=caxis;end;caxis(cx); +colormap(mycolormap(mapanom,21)); +c=colorbar(CBAROR); +titre1=strcat('SVD',num2str(mod(iN)),':',blanks(1),ffile2); +titre2=strcat(' (',num2str(expvar(mod(iN))),'%)'); +tt=title(strvcat(strcat(titre1,titre2),char(1))); +set(tt,'fontweight','bold') + +end %for iN + + +% $$$ for iN=1:2*N +% $$$ figure(f1); +% $$$ subplot(iw,jw,iN); +% $$$ % colorbar(CBAROR); +% $$$ end +% $$$ refresh + + +%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2))); +titre2=strcat(' Total explained variance:',num2str(sum(expvar(mod))),'%'); +%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')'); +%titre=strvcat(titre2,titre3); +suptitle(titre2); + +set(f1,'Position',[378 39 313 647]); % Laptop screen +%set(f1,'Position',[369 55 316 899]); % Desktop screen +%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART')); Index: /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDlag.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDlag.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDlag.m (revision 21589) @@ -0,0 +1,62 @@ +% dispSVDlag(typg,CHP1,CHP2,N,expvar,LAGS,dt) Display SVDs lagged +% +% => Display SVDs lagged +% typg is plot type you want (see plotm.m). +% CHP contains SVDs and is of the form: CHP(iSVD,X,Y) +% N is number of SVD to display +% expvar contained the explained variance of SVDs +% +% Rq: it''s using the plotm function +% +%================================================================ + +% Guillaume MAZE - LPO/LMD - December 2004 +% gmaze@univ-brest.fr + + function [] = dispSVD(typg,CHP1,CHP2,N,expvar,LAGS,dt) + +load mapanom2 +global ffile1 ffile2 +% -------------------------------------------------- +% AFFICHAGE +% -------------------------------------------------- +if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6) + CBAROR='horiz'; +else + CBAROR='vert'; +end + + CHP1=squeeze(CHP1(:,N,:,:)); + CHP2=squeeze(CHP2(:,N,:,:)); + +Nlag=size(CHP1,1); + +iw=round(Nlag/dt);jw=2; + + +f1=figur;clf; hold on +ilag=0; +for lag=1:dt:Nlag + ilag=ilag+1; + + subplot(iw,jw,2*ilag-1) +C = squeeze(real(CHP1(lag,:,:))); +plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]); +if(ilag==1);cx=caxis;end;caxis(cx); +colormap(mycolormap(mapanom,21)); +c=colorbar(CBAROR); +titre1=strcat('lag:',num2str(LAGS(lag)),'y: ',ffile1); +title(titre1); + + subplot(iw,jw,2*ilag) +C = squeeze(real(CHP2(lag,:,:))); +plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]); +if(ilag==1);cx=caxis;end;caxis(cx); +colormap(mycolormap(mapanom,21)); +c=colorbar(CBAROR); +titre1=strcat('lag:',num2str(LAGS(lag)),'y: ',ffile2); +title(titre1); + +end %for iN + +suptitle(strcat('Lagged SVD, mode:',num2str(N))); Index: /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDnb.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDnb.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDnb.m (revision 21589) @@ -0,0 +1,63 @@ +% dispSVDnb(typg,CHP1,CHP2,N,expvar) Display SVDs +% +% => Display SVDs in white and black +% typg is plot type you want (see plotm.m). +% CHP contains SVDs and is of the form: CHP(iSVD,X,Y) +% N is number of SVD to display +% expvar contained the explained variance of SVDs +% +% Rq: it''s using the plotmNB function +% +%================================================================ + +% Guillaume MAZE - LPO/LMD - December 2004 +% gmaze@univ-brest.fr + +function [] = dispSVDnb(typg,CHP1,CHP2,N,expvar) + +global ffile1 ffile2 +% -------------------------------------------------- +% AFFICHAGE +% -------------------------------------------------- +if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6) + CBAROR='horiz'; + iw=N;jw=2; +else + CBAROR='vert'; + iw=N;jw=2; +end + + +f1=figur;clf; hold on +for iN=1:N + subplot(iw,jw,2*iN-1) +C = squeeze(real(CHP1(iN,:,:))); +%[cs,h]=plotmNB(C,typg); whos cs h +[cs,h]=plotmNB(C,typg); +%clabel(cs,h,'labelspacing',100); +titre1=strcat('SVD',num2str(iN),':',ffile1); +titre2=strcat(' (',num2str(expvar(iN)),'%)'); +title(strcat(titre1,titre2)); + + subplot(iw,jw,2*iN) +C = squeeze(real(CHP2(iN,:,:))); +[cs,h]=plotmNB(C,typg); +%clabel(cs,h,'labelspacing',100); +titre1=strcat('SVD',num2str(iN),':',ffile2); +titre2=strcat(' (',num2str(expvar(iN)),'%)'); +title(strcat(titre1,titre2)); + +end %for iN + + + + +%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2))); +titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%'); +%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')'); +%titre=strvcat(titre2,titre3); +suptitle(titre2); + +set(f1,'Position',[378 39 313 647]); % Laptop screen +%set(f1,'Position',[369 55 316 899]); % Desktop screen +%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART')); Index: /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDnbc.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDnbc.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/dispSVDnbc.m (revision 21589) @@ -0,0 +1,66 @@ +% dispSVDnb(typg,CHP1,CHP2,N,expvar) Display SVDs +% +% => Display SVDs in white and black with correlation form +% typg is plot type you want (see plotm.m). +% CHP contains SVDs and is of the form: CHP(iSVD,X,Y) +% N is number of SVD to display +% expvar contained the explained variance of SVDs +% +% Rq: it''s using the plotmNB function +% +%================================================================ + +% Guillaume MAZE - LPO/LMD - December 2004 +% gmaze@univ-brest.fr + +function [] = dispSVDnbc(typg,CHP1,CHP2,N,expvar) + +global ffile1 ffile2 nbcont +nbcont=[-1:0.2:1]; +nbcont2=[-1:0.4:1]; + +% -------------------------------------------------- +% AFFICHAGE +% -------------------------------------------------- +if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6) + CBAROR='horiz'; + iw=N;jw=2; +else + CBAROR='vert'; + iw=N;jw=2; +end + + +f1=figur;clf; hold on +for iN=1:N + subplot(iw,jw,2*iN-1) +C = squeeze(real(CHP1(iN,:,:))); +[cs,h]=plotmNB(C,typg); +clabel(cs,h,nbcont2,'labelspacing',100,'fontsize',9); +titre1=strcat('SVD',num2str(iN),':',ffile1); +titre2=strcat(' (',num2str(expvar(iN)),'%)'); +title(strcat(titre1,titre2)); + + subplot(iw,jw,2*iN) +C = squeeze(real(CHP2(iN,:,:))); +[cs,h]=plotmNB(C,typg); +clabel(cs,h,nbcont2,'labelspacing',100,'fontsize',9); +%clabel(cs,h,'labelspacing',100); +titre1=strcat('SVD',num2str(iN),':',ffile2); +titre2=strcat(' (',num2str(expvar(iN)),'%)'); +title(strcat(titre1,titre2)); + +end %for iN + + + + +%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2))); +titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%'); +%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')'); +%titre=strvcat(titre2,titre3); +suptitle(titre2); + +set(f1,'Position',[378 39 313 647]); % Laptop screen +%set(f1,'Position',[369 55 316 899]); % Desktop screen +%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART')); Index: /issm/trunk-jpl/externalpackages/pcatool/install/license.txt =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/license.txt (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/license.txt (revision 21589) @@ -0,0 +1,24 @@ +Copyright (c) 2008, Guillaume MAZE +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. Index: /issm/trunk-jpl/externalpackages/pcatool/install/map2mat.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/map2mat.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/map2mat.m (revision 21589) @@ -0,0 +1,40 @@ +% D = MAP2MAT(F,C) Reshaping matrix +% +% => Conversion of a 'map' matrix C(TIME,LON,LAT) into a D(TIME,PT) matrix +% under the mask F(LON,LAT). +% F is a matrix contenaing 1 where you would like to keep the point and +% 0 elsewhere. +% +% Rq: No check is done about the input. +% +% See also: mat2map +%================================================================ + +% March 2004 +% gmaze@univ-brest.fr + +function [D] = map2mat(F,C); + +% Get dimensions +[tps nolon nolat] = size(C); + +% So output matrix will be: +D = zeros(tps,nolon*nolat); + +% point indice +ipt = 0; + +% 'Un-mapping' : +for iy=1:nolat + for ix=1:nolon + if F(ix,iy)>0 + ipt = ipt + 1; + D(:,ipt)=squeeze(C(:,ix,iy)); + end % if + end % for +end %for + +% OUTPUT: +D = squeeze(D(:,1:ipt)); + + Index: /issm/trunk-jpl/externalpackages/pcatool/install/mat2map.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/mat2map.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/mat2map.m (revision 21589) @@ -0,0 +1,41 @@ +% C = MAP2MAT(F,D) Reshaping matrix +% +% => Creation of a 'map' matrix C(TIME,LON,LAT) from D(TIME,PT) +% following mask F(LON,LAT). +% F is a matrix contenaing 1 where you would like to keep the point and +% 0 elsewhere (see mat2map). +% +% Rq: No check is done about the input. +% +% See also: map2mat +%================================================================ + +% March 2004 +% gmaze@univ-brest.fr + +function [C] = mat2map(F,D); + +% Get dimensions +[nolon nolat] = size(F); +[time npt] = size(D); + +% So output 'map' matrix has the form: +C = zeros(time,nolon,nolat); + +% Variables +nul = NaN.*ones(time,1); +ipt = 0 ; + + +% 'mapping' : + for iy=1:nolat + for ix=1:nolon + if F(ix,iy)>0 + ipt = ipt + 1; + C(:,ix,iy) = D(:,ipt); + else + C(:,ix,iy) = nul; + end %if + end %for ix + end %for iy + Index: /issm/trunk-jpl/externalpackages/pcatool/install/plotm.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/plotm.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/plotm.m (revision 21589) @@ -0,0 +1,214 @@ +%PLOTM Trace un champ de donnee geographique +% [CS,H]=PLOTM(C,[ZONE]) +% Trace le champ C(LONGITUDE,LATITUDE) sur une carte +% du monde en projection Mercator par defaut. +% Le parametre optionnel ZONE peut prendre les valeurs +% ZONE = 1 : Projection Stereographique centree sur le pole sud +% ZONE = 2 : Projection Mercator de l'hemisphere sud +% ZONE = 3 : Projection Mercator du globe +% ZONE = 4 : Pas de projection (imagesc) +% ZONE = 5 : Projection Mercator de l'hemisphere sud complet +% ZONE = 6 : Projection Mercator du globe complet +% ZONE = 7 : Projection Stereographique centree sur le pole nord +% ZONE = 8 : Idem 2 mais sans label sur les axes +% ZONE = 9 : Projection Mercator de l'hemisphere sud aux limites du Fmask +% Returns contour matrix CS as described in CONTOURC and a vector H +% of handles to PATCH objects (for use by CLABEL). +% +% La definition des grilles (longitude,latitude) est celle du +% modele par defaut, mais elle peut etre modifiee en passant en +% global les nouvelles variables nommees lon et lat. +% +% NECESSITE LE PACKAGE M_MAP +% +% +% (nov. 2003 -> gmaze@univ-brest.fr) + +function varargout = plotm(C,varargin) + +if (nargin>2)|(nargin<1) + help plotm.m + error('plotm.m : Wrong number or bad parameter(s)') + return +elseif (nargin==2) + arg = varargin{:}; + if (arg>14)|(arg<1) + help plotm.m + error('plotm.m : Wrong number or bad parameter(s)') + return + end +end %if + +%-------------------------------------------------- +% Variabes +%-------------------------------------------------- +[lx,ly]=size(C); +% Si lon et lat ne sont pas definie comme dans le modele +% on prend en global n'importe quel definition +if ((lx~=64)|(ly~=32)) + global lon lat +else + lon= linspace(0.,360.,64); + lat= [85.76 80.27 74.74 69.21 63.68 58.14 52.61 47.07 41.53 ... + 36.00 30.46 24.92 19.38 13.84 8.31 2.77 -2.77 -8.31 ... + -13.84 -19.38 -24.92 -30.46 -36.00 -41.53 -47.07 -52.61 ... + -58.14 -63.68 -69.21 -74.74 -80.27 -85.76]; +end + +% Le nombre de contour est par defaut a 10, mais on peut le +% modifier en declarant nbcont en global +global nbcont +if isempty(nbcont)==1 + nbcont=10; +end + +%----------------------------------------------- +% Projection +%----------------------------------------------- +if (nargin==2) + arg = varargin{:}; + switch arg(1) + case 1 % Projection stereo + m_proj('stereo','lon',0,'lat',-90,'rad',71); + case 2 % Projection mercator (HS) + m_proj('miller','lon',[0 360],'lat',[-80 -20]); + m_proj('miller','lon',[0 360],'lat',[-60 -20]); + case 3 % Projection mercator (Tout) -> Defaut + m_proj('miller','lon',[0 360],'lat',[-80 80]); + case 5 % Projection mercator (HS-complet) + m_proj('miller','lon',[0 360],'lat',[-85 -5]); + case 6 % Projection mercator (Tout-complet) + m_proj('miller','lon',[0 360],'lat',[-85 85]); + case 7 % Projection stereo + m_proj('stereo','lon',0,'lat',90,'rad',71); + case 8 % Projection mercator (HS) sans Label sur les axes + m_proj('miller','lon',[0 360],'lat',[-80 -20]); + m_proj('miller','lon',[0 360],'lat',[-75 -30]); + case 9 % Projection mercator (HS) sans Label sur les axes aux limtes du FMASK + m_proj('miller','lon',[0 360],'lat',[-69 -25]); + case 10 % Projection stereo recentree sur le pole sud + m_proj('stereo','lon',0,'lat',-90,'rad',51); + case 11 % Indien + m_proj('lambert','long',[11 125],'lat',[-60 -30]); + case 12 % Pacifique + m_proj('lambert','long',[120 290],'lat',[-60 -30]); + case 13 % Atlantique + m_proj('lambert','long',[-70 20],'lat',[-60 -30]); + case 14 % Centré Atlantique + m_proj('lambert','long',[-90 90],'lat',[-70 -20],'rectbox','on'); + end +else +% Projection mercator (Tout) -> Defaut + m_proj('mercator','lon',[0 360],'lat',[-80 80]); +end + +%------------------------------------------------ +% Affichage +%------------------------------------------------ +load MODELCOAST +if length(lon)==65 + modelcoast=[modelcoast;modelcoast(1,:)]; +end + +hold on +if (nargin==2) + arg = varargin{:}; + switch arg(1) + case 1 % Projection stereo (Par defaut) + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-30 -40 -50 -60 -70],'yticklabels',[],... + 'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top') + case 2 % Projection mercator (HS) + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-80 -70 -60 -50 -40 -30 -20]); + case 3 % Projection mercator (Tout) + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]); + case 4 % Pas de projection + pcolor(lon,lat,C'); + axis tight;box on +% axis([lon(1) lon(64) lat(32) lat(1)]); +% gc=gca;set(gc,'YDir','reverse'); + case 5 % Projection mercator (Hs-complet) + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + COASTYPE='color';COASTCOLOR=[0 0 0]; +% COASTYPE='patch';COASTCOLOR=[0 0.5 0.6]; +% m_coast(COASTYPE,COASTCOLOR); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-80 -60 -40 -20 0]); + case 6 % Projection mercator (Tout-complet) + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]); + case 7 % Projection stereo + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[30 40 50 60 70],'yticklabels',[],... + 'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top') + case 8 % Projection mercator (HS) sans label sur les axes + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_pcolor(lon,lat,C');shading interp +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick','','xtick',''); + case 9 % Projection mercator (HS) sans label sur les axes aux limtes du FMASK + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick','','xtick',''); + case 10 % Projection stereo recentree sur le pole Sud + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],... + 'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top') + case 11 % Indien + [cs,h]=m_contourf(lon,lat,C',nbcont); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],... + 'xtick',[30 60 90 120 150 180],'XaxisLocation','top') + case 12 % Pacifique + [cs,h]=m_contourf(lon,lat,C',nbcont); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],... + 'xtick',[120 150 180 210 240 270 300],'XaxisLocation','top') + case 13 % Atlantique + [cs,h]=m_contourf(lon,lat,C',nbcont); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],... + 'xtick',[-60 -30 0 30 60],'XaxisLocation','top') + case 14 % Centré Atlantique + [cs,h]=m_contourf(lon,lat,C',nbcont); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('yticklabels',[],'XaxisLocation','top',... + 'fontsize',6,... + 'xtick',[0 60 120 180 -120 -60],'xticklabels',[180 -120 -60 0 60 120]); + + end % switch +else + % Projection Mercator (Par defaut) + [cs,h]=m_contourf(lon,lat,C',nbcont); +% m_coast('color',[0 0 0]); + m_contour(lon,lat,modelcoast',1,'k'); + m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]); +end + + + +% OUTPUT VARIABLES +switch nargout + case 1 + varargout(1) = {cs} ; + case 2 + varargout(1) = {cs} ; + varargout(2) = {h}; +end Index: /issm/trunk-jpl/externalpackages/pcatool/install/test_pcatool.m =================================================================== --- /issm/trunk-jpl/externalpackages/pcatool/install/test_pcatool.m (revision 21589) +++ /issm/trunk-jpl/externalpackages/pcatool/install/test_pcatool.m (revision 21589) @@ -0,0 +1,160 @@ +% Here is a simple script to test caleof.m from PCATOOL package. +% + +% gmaze@mit.edu +% 2006/07/31 +% + +clear + +field = 1; % The best one + +% Generate the data +x = -10 : 1 : 10 ; +y = x; +[X,Y] = meshgrid(x,y); + +switch field + case 1 % H. Bjornson and S.A. Venegas + HH = 12*( 1.2 - 0.35*sqrt( (X-0).^2 + (Y-0).^2 ) ); + HH = HH - mean(mean(HH)); + H(1,:,:) = 1012 + HH; + H(2,:,:) = 1012 + (1012-squeeze(H(1,:,:)));; + H(3,:,:) = 1022.8 - 3.6*Y; + H(4,:,:) = 1012 + (1012-squeeze(H(3,:,:))); + H(5,:,:) = 1001.2 + 3.6*X; + H(6,:,:) = 1012 + (1012-squeeze(H(5,:,:))); + cont = [980:5:1050]; + N = 2; + + case 11 % H. Bjornson and S.A. Venegas modified + % Here it must be only 2 non-zero EOF (45degre oriented ) of similar variance 50% + H = ones(4,length(x),length(y)); + H(1,:,:) = Y; + H(2,:,:) = -Y; + H(3,:,:) = X; + H(4,:,:) = -X; + N = 2; + + case 2 % Hartmann eg: analysis in place + a = [2 4 -6 8 ; 1 2 -3 4]; + H(1,:,:) = a(1,:)'; + H(2,:,:) = a(2,:)'; + + % Get EOFs: + N = 4; + for method = 2 : -1 : 1 + [E,pc,expvar] = caleof(H,N,method); + %E' + end + + return + + case 3 % My signal + np = 4; % Numer of random signal (...and eof) + nt = 10; % Number of timestep + H = zeros(nt,length(x),length(y)); + for ip = 1 : np + xc(ip) = abs(fix(rand(1)*10)); + yc(ip) = abs(fix(rand(1)*10)); + if ip>=fix(np/2) + xc(ip) = -xc(ip); + yc(ip) = -yc(ip); + end + dd(ip) = fix(rand(1)*10); + end %for ip + f = 1/nt; + for it = 1 : nt + H2 = zeros(length(x),length(y)); + for ip = 1 : ip + if it==1,[xc(ip) yc(ip) dd(ip)],end + HH = 12*( 1.2 - 0.35*sqrt( ((X-xc(ip)).^2 + (Y-yc(ip)).^2 )/dd(ip) )); + H2 = HH - mean(mean(HH)); + H(it,:,:) = squeeze(H(it,:,:)) + H2.*cos(pi*it/dd(ip)) ; + end %for ip + end %for it + H = 1012 + H; + cont = [980:5:1050]; + N = 3; + + case 4 % My signal 2 + x = -pi : pi/6 : pi ; + y = x; + [X,Y] = meshgrid(x,y); + HH = cos(X) + cos(Y); + HH2 = cos(Y); + nt = 12; + for it = 1 : nt +% H(it,:,:) = cos(pi*it/nt)*HH; cont=[-2 2]; + H(it,:,:) = 2*cos(pi*it/nt)*HH + 3*cos(pi*it/nt/2)*HH2; cont=[-10 10]; + xtrm(squeeze(H(it,:,:))); + end + cont=[-2 2]; + N = 3; + +end % switch field + +%return + +% Plot field time serie: +if 1 +figure;iw=2;jw=size(H,1)/iw; +set(gcf,'position',[11 533 560 420]); + +for i=1:iw*jw + C = squeeze(H(i,:,:)); + + subplot(iw,jw,i); + pcolor(X,Y,C); + title(num2str(i)); + if i==1,cx=caxis;end + axis square + caxis(cx); + +end %for i +end %fi +%return + +% Get EOFs: +G = map2mat(ones(size(H,2),size(H,3)),H); + +for method = 1 : 4 + [E,pc,expvar] = caleof(G,N,method); + eof = mat2map(ones(size(H,2),size(H,3)),E); + + +figure;iw=1;jw=N+1; +set(gcf,'MenuBar','none'); +posi = [576 0 710 205]; +set(gcf,'position',[posi(1) (method-1)*250+50 posi(3) posi(4)]); + +for i=1:iw*jw + if i<= iw*jw-1 + + C = squeeze(eof(i,:,:)); + subplot(iw,jw,i); + cont = 12; + [cs,h] = contourf(X,Y,C,cont); + clabel(cs,h); + title(strcat('EOF:',num2str(i),'/',num2str(expvar(i)),'%')); + axis square; + %caxis([cont(1) cont(end)]); + + else + subplot(iw,jw,iw*jw); + plot(pc'); + grid on + xlabel('time') + title('PC') + legend(num2str([1:N]'),2); + box on + + + end %if + +end %for i +%suptitle(strcat('METHODE:',num2str(method))); + +end %for method + +