source: issm/oecreview/Archive/21337-21723/ISSM-21588-21589.diff@ 21726

Last change on this file since 21726 was 21726, checked in by Mathieu Morlighem, 8 years ago

CHG added Archive/21337-21723
File size: 53.9 KB

  • ../trunk-jpl/externalpackages/pcatool/install.sh

     
     1#!/bin/bash
     2
     3#Some cleanup
     4rm -rf install pcatool
     5mkdir install
     6
     7#Download from ISSM server
     8$ISSM_DIR/scripts/DownloadExternalPackage.py 'http://issm.jpl.nasa.gov/files/externalpackages/pcatool.tar.gz' 'pcatool.tar.gz'
     9
     10#Untar  into install
     11cd install
     12tar -zxvf  ../pcatool.tar.gz

  • ../trunk-jpl/externalpackages/pcatool/install/license.txt 

    Property changes on: ../trunk-jpl/externalpackages/pcatool/install.sh
___________________________________________________________________
Added: svn:executable
   + *

Cannot display: file marked as a binary type.
svn:mime-type = application/octet-stream

Property changes on: ../trunk-jpl/externalpackages/pcatool/pcatool.tar.gz
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream
     
     1Copyright (c) 2008, Guillaume MAZE
     2All rights reserved.
     3
     4Redistribution and use in source and binary forms, with or without
     5modification, are permitted provided that the following conditions are
     6met:
     7
     8    * Redistributions of source code must retain the above copyright
     9      notice, this list of conditions and the following disclaimer.
     10    * Redistributions in binary form must reproduce the above copyright
     11      notice, this list of conditions and the following disclaimer in
     12      the documentation and/or other materials provided with the distribution
     13
     14THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
     15AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     16IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     17ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
     18LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
     19CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
     20SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
     21INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
     22CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
     23ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
     24POSSIBILITY OF SUCH DAMAGE.

  • ../trunk-jpl/externalpackages/pcatool/install/dispEEOF_v2.m

     
     1% dispEEOF_v2(CHP,PC,EXPVAR,DT,NLAG,MOD) Display few EEOFs with their
     2%                                     PC time evolution
     3%
     4% => DISPLAY FEW EEOFs with their PC time evolution.
     5% CHP contains all the EEOFs as EOF*LAG*X*Y.
     6% PC contains all the PCs.
     7% EXPVAR is a matrix with the explained variance of each
     8%  EEOFs in %. This is just for title.
     9% DT is time step between maps.
     10% NLAG is the number of LAG to display.
     11% MOD contains explicitly the EEOFs to display.
     12%
     13% Rq: it's using the plotm function
     14% See also: dispEEOF
     15%
     16%================================================================
     17
     18%  Guillaume MAZE - LPO/LMD - June 2004
     19%  gmaze@univ-brest.fr
     20
     21function [] = dispEEOF(CHP,PC,EXPVAR,DT,NLAG,MOD)
     22
     23% Number of EEOF to display     
     24NMOD = length(MOD);
     25
     26% Open figure and first guest
     27figur;
     28clf;hold on
     29typg=8;
     30width  = .9/NMOD;
     31height = .9/(NLAG+1);
     32dleft = (.95-width*NMOD)/4;
     33left  = dleft;
     34load mapanom
     35
     36% Let'go :
     37
     38% We choose an EEOF
     39for imod = 1 : NMOD
     40    mod = MOD(imod);
     41   
     42    % and plot maps
     43    for iplot=1:NLAG+1
     44      bottom = .95-height*iplot;
     45      subplot('Position',[left bottom width height]);
     46      if iplot==1 % Plot the pc's time serie
     47         plot(PC(imod,:));
     48         grid on;box on;
     49         set(gca,'XAxisLocation','top');axis tight;
     50         titre=strcat('EEOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)');
     51         titre=strcat(titre,'; DT=',num2str(DT),' ');
     52         title(titre);
     53         if (imod==1),ax1=get(gca,'YLim');else,set(gca,'YLim',ax1);end
     54      else % Plot the lag maps
     55         lag=iplot-1;
     56         C = squeeze(real(CHP(mod,lag,:,:)));
     57%         C = C./xtrm(C); % Eventually normalise field
     58         plotm(C,typg);
     59         if(lag==1),caxis([-abs(xtrm(C)) abs(xtrm(C))]);cx=caxis;end;
     60         caxis(cx);
     61         colormap(mapanom);
     62         %ylabel(strcat('{\bf',num2str((lag-1)*DT),'}'))
     63         m_text(10,-70,strcat('{\bf',num2str((lag-1)*DT),'}'));
     64      end %if
     65
     66    end
     67    left = left + width + dleft;
     68
     69    % Adjust plot width via the plot Aspect ratio
     70    ch=get(gcf,'Children');
     71    posiPC=get(ch(iplot),'PlotBoxAspectRatio');
     72    posiMAP=get(ch(iplot-1),'PlotBoxAspectRatio');
     73    set(ch(iplot),'PlotBoxAspectRatio',posiMAP);
     74
     75end
     76

  • ../trunk-jpl/externalpackages/pcatool/install/confexpvar.m

     
     1% [CONF,CONFLEVEL] = CONFEXPVAR(A,B,nsamples,N)
     2%
     3% A(Ntime,MmapA), B(Ntime,MmapB) are time-synchronised
     4%  and time-constant sampled data sets.
     5% MmapA and MmapB are not required to be identic.
     6% nsamples is what it is (given at least a 100 factor
     7%  will be nice to provide a trustable result).
     8% N is the max number of eigen vectors to be computed.
     9%
     10% CONFEXPVAR Compute the confidence limits for significant
     11% eigenvalues in a SVD decomposition of data sets A and B.
     12% This generates nsamples random-in-time data set of B
     13% and finds percent variance of N SVD modes computed
     14% with data set A.
     15% This is exactly the monte-carlo method but based on
     16% a random-in-time only data set. This considerably
     17% decrease the cputime.
     18%
     19% 2005/11. Guillaume Maze
     20% gmaze@univ-brest.fr
     21
     22function [conf,lev,dsumCF] = confexpvar(A,B,nsamples,N)
     23 
     24% We try to find the p confidence levels of expvar.
     25% Keeping as it is the first field, we determine SVDs with the
     26% second field arbitrary sorted in time.
     27% 1st field : A
     28% 2nd field : B
     29
     30n = size(A,1);
     31if size(B,1)~=n
     32  disp('Error, A and B required to be of the same time length !');
     33  disp('See help confexpvar');
     34  stop
     35end
     36 
     37timer = cputime;
     38disp('Compute confidence levels of explained variance...');
     39for iR=1:nsamples
     40    if iR==fix(nsamples/2)
     41       disp(strcat(num2str(iR),'/',num2str(nsamples)));
     42       disp(['Please, still wait a little time ... around '...
     43             num2str(cputime-timer) ' cpu seconds']);
     44    end
     45   
     46    % Generate a random new time index:
     47    t=randperm(n);
     48    % and sort the 2nd field with it:
     49    Br=B(t,:);
     50
     51    % Now we compute SVD and expvar:
     52    C=A'*Br;
     53    [U,Lambda,V] = svds(C,N);
     54    L2 = Lambda.^2;   dsum = diag(L2)/trace(L2);
     55    L  = Lambda;      dsum =  diag(L)/trace(L);
     56    dsumCF(iR,:) = dsum';
     57
     58end %for iR
     59disp(['OK, this took ' num2str(cputime-timer) ' cpu seconds'])
     60
     61% CF(nsamples,N) contains expvar of the nsamples realisations:
     62CF=sort(dsumCF,1);
     63
     64% Now we just take the decadal significant levels:
     65lev=fix([nsamples/10:nsamples/10:nsamples]);
     66conf = CF(lev,:) ;
     67lev=(lev./max(lev))*100;

  • ../trunk-jpl/externalpackages/pcatool/install/calsvd.m

     
     1% [SVD1,SVD2,PC1,PC2,EXPVAR,Lambda] = CALSVD(A,B,N) Compute SVDs
     2%
     3% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD -
     4%      Analyses of climatic Data" 1997
     5%================================================================
     6%
     7%  Guillaume MAZE - LPO/LMD - March 2004
     8%  gmaze@univ-brest.fr
     9
     10
     11   function [e1,e2,pc1,pc2,expvar,Lambda] = calsvd(A,B,N);
     12
     13
     14disp('====> Let''go for SVDs and pc computing...')
     15
     16% Assume that A is (time*map) matrix
     17[n p]=size(A);
     18
     19% Remove the mean of each column (ie the time mean in each station records)
     20S=detrend(A,'constant');
     21P=detrend(B,'constant');
     22
     23% Form the covariance matrix:
     24C = S'*P;
     25
     26% Find eigenvectors and singular values
     27[U,Lambda,V] = svds(C,N);
     28
     29% PC
     30a = S*U;
     31b = P*V;
     32
     33% Make them clear for output
     34for iN=1:N
     35    e1(iN,:) = squeeze( U(:,iN) )';
     36   pc1(iN,:) = squeeze( a(:,iN) )';
     37    e2(iN,:) = squeeze( V(:,iN) )';
     38   pc2(iN,:) = squeeze( b(:,iN) )';
     39end
     40
     41% Amount of variance explained a 0.1 pres et en %
     42L2=Lambda.^2;
     43dsum=diag(L2)/trace(L2);
     44
     45for iN=1:N
     46   expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10;
     47end
     48
     49
     50disp('====> Finished !')
     51
     52
     53
     54

  • ../trunk-jpl/externalpackages/pcatool/install/dispSVD.m

     
     1% dispSVD(typg,CHP1,CHP2,N,expvar) Display SVDs
     2%
     3% => Display SVDs
     4% typg is plot type you want (see plotm.m).
     5% CHP contains SVDs and is of the form: CHP(iSVD,X,Y)
     6% N is number of SVD to display
     7% expvar contained the explained variance of SVDs
     8%
     9% Rq: it''s using the plotm function
     10%
     11%================================================================
     12
     13%  Guillaume MAZE - LPO/LMD - December 2004
     14%  gmaze@univ-brest.fr
     15
     16function [] = dispSVD(typg,CHP1,CHP2,N,expvar)
     17
     18load mapanom2
     19global ffile1 ffile2
     20
     21if length(N)>1
     22  mod=N;
     23else
     24  mod=[1:N];
     25end 
     26Nx=length(mod);
     27% --------------------------------------------------
     28% AFFICHAGE
     29% --------------------------------------------------
     30if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
     31   CBAROR='horiz';
     32   iw=Nx;jw=2;
     33else
     34   CBAROR='vert';
     35   iw=Nx;jw=2;
     36end
     37
     38
     39f1=figur;clf; hold on
     40for iN=1:Nx
     41     subplot(iw,jw,2*iN-1)
     42C = squeeze(real(CHP1(mod(iN),:,:)));
     43%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     44plotm(C,typg);
     45caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     46%if(iN==1);cx=caxis;end;caxis(cx);
     47colormap(mycolormap(mapanom,21));
     48c=colorbar(CBAROR);
     49titre1=strcat('SVD',num2str(mod(iN)),':',blanks(1),ffile1);
     50titre2=strcat(' (',num2str(expvar(mod(iN))),'%)');
     51tt=title(strvcat(strcat(titre1,titre2),char(1)));
     52set(tt,'fontweight','bold')
     53
     54     subplot(iw,jw,2*iN)
     55C = squeeze(real(CHP2(mod(iN),:,:)));
     56%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     57plotm(C,typg);
     58caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     59%if(iN==1);cx=caxis;end;caxis(cx);
     60colormap(mycolormap(mapanom,21));
     61c=colorbar(CBAROR);
     62titre1=strcat('SVD',num2str(mod(iN)),':',blanks(1),ffile2);
     63titre2=strcat(' (',num2str(expvar(mod(iN))),'%)');
     64tt=title(strvcat(strcat(titre1,titre2),char(1)));
     65set(tt,'fontweight','bold')
     66
     67end %for iN
     68
     69
     70% $$$ for iN=1:2*N
     71% $$$   figure(f1);
     72% $$$   subplot(iw,jw,iN);
     73% $$$ %  colorbar(CBAROR);
     74% $$$ end
     75% $$$ refresh
     76
     77
     78%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
     79titre2=strcat(' Total explained variance:',num2str(sum(expvar(mod))),'%');
     80%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')');
     81%titre=strvcat(titre2,titre3);
     82suptitle(titre2);
     83
     84set(f1,'Position',[378 39 313 647]); % Laptop screen
     85%set(f1,'Position',[369 55 316 899]); % Desktop screen
     86%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART'));

  • ../trunk-jpl/externalpackages/pcatool/install/caleof.m

     
     1% [EOFs,PC,EXPVAR] = CALEOF(M,N,METHOD) Compute EOF
     2%
     3% => Compute the Nth first EOFs of matrix M(TIME,MAP).
     4% EOFs is a matrix of the form EOFs(N,MAP), PC is the principal
     5% components matrix ie it has the form PC(N,TIME) and EXPVAR is
     6% the fraction of total variance "explained" by each EOF ie it has
     7% the form EXPVAR(N).
     8% Differents method can be used:
     9% 1 - The "classic" one, ie we compute eigenvectors of the
     10%     temporal covariance matrix with the eig Matlab function.
     11% 2 - A faster "classic" one, same as method 1 but we use the
     12%     eigs Matlab function.
     13% 3 - We compute eigenvectors by the singular value decomposition,
     14%     by used of the svd Matlab function.
     15% 4 - Same as method 3 but faster by used of the svds Matlab function
     16%
     17% See also EIG, EIGS, SVD, SVDS
     18%
     19% Ref: L. Hartmann: "Objective Analysis" 2002
     20% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD -
     21%      Analyses of climatic Data" 1997
     22%================================================================
     23
     24%  Guillaume MAZE - LPO/LMD - March 2004
     25%  Revised July 2006
     26%  gmaze@univ-brest.fr
     27
     28
     29function [e,pc,expvar,L] = caleof(M,N,method);
     30
     31% Get dimensions
     32[n p]=size(M);
     33
     34% Temporal covariance is p*p matrix, that why max EOF computable is p,
     35% so we perform a test on parameter N:
     36if(N>p)
     37 disp('Warning: N is larger than possible so it''s modified to perform')
     38 disp('         EOFs computing...');
     39 N = p;
     40end
     41
     42
     43% Eventualy time filtering of data
     44if 0==1
     45   disp('====> Time filtering...')
     46   Fc  = 1/20; Fc2 = 1/1;
     47   Fc  = 1/7 ; Fc2 = 1/3;
     48   SIGNAL = M(:,1);
     49   nj = fix(length(SIGNAL)/10); % Nombre de points du filtre
     50   for ipt = 1 : p
     51       SIGNAL = M(:,ipt);
     52       SIGNALF = lanczos(SIGNAL,Fc2,nj);
     53       SIGNALF = SIGNALF - lanczos(SIGNALF,Fc,nj);
     54       Y(:,ipt) = SIGNALF;
     55   end
     56   M = Y;
     57end
     58
     59
     60disp('====> Let''go for EOFs and pc computing...')
     61switch method
     62    case 1 % CLASSIC METHOD
     63%================================================================
     64% Transform the data matrix in the correct form (map*time) for eig
     65M = M';
     66
     67% Remove the time mean (ie the mean of each rows of M)
     68% Rq: detrend remove the mean of columns ie we take M'.
     69F = detrend(M','constant')';
     70
     71% Covariance Matrix (inner product over space = covariance in time)
     72R = F * F';
     73
     74% Eigenanalysis of the covariance matrix R
     75[E,L] = eig(R);
     76
     77% Get PC by projecting eigenvectors on original data
     78Z = E'*F;
     79
     80% Make them clear for output
     81for iN=1:N
     82    e(iN,:) = squeeze( E(:,p-(iN-1)) )';
     83   pc(iN,:) = squeeze( Z(p-(iN-1),:) );
     84end
     85
     86% Amount of explained variance (at 0.1%)
     87dsum = diag(L)./trace(L);
     88for iN=1:N
     89   expvar(iN)=fix((dsum(p-(iN-1))*100/sum(dsum))*10)/10;
     90end
     91
     92% Plots Original field and reconstructed one
     93if 0==1
     94figure;
     95subplot(1,2,1);imagesc(abs(M));title('ORIGINAL');cx=caxis;
     96%subplot(1,2,2);imagesc((E*Z));title('RECONSTRUCTED')
     97subplot(1,2,2);imagesc(abs(e'*pc));title('RECONSTRUCTED');caxis(cx);
     98end
     99
     100    case 2 % RAPID CLASSIC METHOD
     101%================================================================
     102% Remove the time mean of each column
     103F = detrend(M,'constant');
     104
     105% Covariance Matrix
     106if n >= p
     107   R = F' * F;
     108else
     109   R = F * F';
     110end
     111
     112% Eigen analysis of the square covariance matrix
     113[E,L] = eigs(R,N);
     114if n < p
     115  E = F' * E;
     116  sq = [sqrt(diag(L))+eps]';
     117  sq = sq(ones(1,p),:);
     118  E = E ./ sq;
     119end
     120
     121% Get PC by projecting eigenvectors on original data
     122if n >= p
     123   Z = (F*E)';
     124else
     125   Z =  E'*F';
     126end
     127
     128
     129% Make them clear for output
     130for iN=1:N
     131    e(iN,:) = squeeze( E(:,iN) )';
     132   pc(iN,:) = squeeze( Z(iN,:) );
     133end
     134
     135% Amount of variance explained a 0.1 pres et en %
     136dsum=diag(L)./trace(L);
     137for iN=1:N
     138   expvar(iN)=fix((dsum(iN)*100/sum(dsum))*10)/10;
     139end
     140
     141
     142    case 3 % SVD METHOD
     143%================================================================
     144% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD -
     145% Analyses of climatic Data" 1997 => p18
     146
     147% Assume that M is (time*map) matrix
     148[n p]=size(M);
     149
     150% Remove the mean of each column (ie the time mean in each station records)
     151F=detrend(M,'constant');
     152
     153% Form the covariance matrix:
     154R = F'*F;
     155
     156% Find eigenvectors and singular values
     157[C,L,CC] = svd(R);
     158% Eigenvectors are in CC and the squared diagonal values of L
     159% are the eigenvalues of the temporal covariance matrix R=F'*F
     160
     161% find the PC corresponding to eigenvalue
     162PC = F*CC;
     163
     164% Make them clear for output
     165for iN=1:N
     166    e(iN,:) = squeeze( CC(:,iN) )';
     167   pc(iN,:) = squeeze( PC(:,iN) )';
     168end
     169
     170if 0
     171figure;
     172subplot(1,2,1);imagesc(F);title('ORIGINAL');cx=caxis;
     173subplot(1,2,2);imagesc(C*L*CC');title('RECONSTRUCTED');caxis(cx);
     174end
     175
     176% Amount of variance explained at 0.1%
     177dsum=diag(L)./trace(L);
     178if length(dsum)<N % L was not squared
     179  dsum = [dsum ;zeros(N-length(dsum),1)];
     180end
     181for iN = 1 : N
     182   expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10;
     183end
     184
     185
     186    case 4 % FAST SVD METHOD
     187%================================================================
     188% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD -
     189% Analyses of climatic Data" 1997 => p18
     190
     191% Assume that M is (time*map) matrix
     192[n p]=size(M);
     193
     194% Remove the mean of each column (ie the time mean in each station records)
     195F=detrend(M,'constant');
     196
     197% Form the covariance matrix:
     198R = F' * F;
     199
     200% Find eigenvectors and singular values
     201[C,L,CC,flag] = svds(R,N);
     202% Eigenvectors are in CC and the squared diagonal values of L
     203% are the eigenvalues of the temporal covariance matrix R=F'*F
     204% (Sometimes, CC stops for nul eigenvector, then we need to fill to reach N)
     205if size(CC,2)<N
     206  CC = [CC  zeros(size(CC,1),N-size(CC,2)+1)];
     207end
     208
     209% find the PC corresponding to eigenvalue
     210PC = F*CC;
     211% Which is similar to: C*L
     212
     213% Make them clear for output
     214for iN=1:N
     215    e(iN,:) = squeeze( CC(:,iN) )';
     216   pc(iN,:) = squeeze( PC(:,iN) )';
     217end
     218
     219% Amount of variance explained a 0.1 pres et en %
     220dsum=diag(L)./trace(L);
     221if length(dsum)<N % L was not squared
     222  dsum = [dsum ;zeros(N-length(dsum),1)];
     223end
     224for iN=1:N
     225   expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10;
     226end
     227
     228%figure;
     229%subplot(1,2,1);imagesc(M);title('ORIGINAL');cx=caxis;
     230%subplot(1,2,2);imagesc((e'*pc)');title('RECONSTRUCTED');caxis(cx);
     231
     232
     233end % switch method
     234disp('====> Finished !')

  • ../trunk-jpl/externalpackages/pcatool/install/dispEOF.m

     
     1% dispEOF(typg,CHP,N,expvar) Display EOFs
     2%
     3% => Display EOFs
     4% typg is plot type you want (see plotm.m).
     5% CHP contains EOFs and is of the form: CHP(EOF,X,Y)
     6% N is number of EOF to display
     7% expvar contained the explained variance of EOFs
     8%
     9% Rq: it''s using the plotm function
     10%
     11%================================================================
     12
     13%  Guillaume MAZE - LPO/LMD - March 2004
     14%  gmaze@univ-brest.fr
     15
     16function [] = dispEOF(typg,CHP,N,expvar)
     17
     18load mapanom2
     19global ffile
     20% --------------------------------------------------
     21% AFFICHAGE
     22% --------------------------------------------------
     23if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
     24
     25   CBAROR='horiz';
     26if N>12,NN=12;else,NN=N;end
     27switch NN
     28  case 1,tblmagique=[1];                         iw=1;jw=1;
     29  case 2,tblmagique=[1 2];                       iw=2;jw=1;
     30  case 3,tblmagique=[1 2 3];                     iw=3;jw=1;
     31  case 4,tblmagique=[1 2 3 4];                   iw=2;jw=2;
     32  case 5,tblmagique=[1 2 3 4 5];                 iw=3;jw=2;
     33  case 6,tblmagique=[1 2 3 4 5 6];               iw=3;jw=2;
     34  case 7,tblmagique=[1 2 3 4 5 6 7];             iw=4;jw=2;
     35  case 8,tblmagique=[1 2 3 4 5 6 7 8];           iw=4;jw=2;
     36  case 9,tblmagique=[1 2 3 4 5 6 7 8 9];         iw=3;jw=3;
     37 case 10,tblmagique=[1 2 3 4 5 6 7 8 9 10];      iw=4;jw=3;
     38 case 11,tblmagique=[1 2 3 4 5 6 7 8 9 10 11];   iw=4;jw=3;
     39 case 12,tblmagique=[1 2 3 4 5 6 7 8 9 10 11 12];iw=4;jw=3;
     40end
     41iw=N;jw=1;
     42
     43else
     44
     45   CBAROR='vert';
     46if N>12,NN=12;else,NN=N;end
     47switch NN
     48  case 1,tblmagique=[1];                         iw=1;jw=1;
     49  case 2,tblmagique=[1 2];                       iw=2;jw=1;
     50  case 3,tblmagique=[1 2 3];                     iw=3;jw=1;
     51  case 4,tblmagique=[1 2 3 4];                   iw=2;jw=2;
     52  case 5,tblmagique=[1 2 3 4 5];                 iw=3;jw=2;
     53  case 6,tblmagique=[1 2 3 4 5 6];               iw=3;jw=2;
     54  case 7,tblmagique=[1 2 3 4 5 6 7];             iw=4;jw=2;
     55  case 8,tblmagique=[1 2 3 4 5 6 7 8];           iw=4;jw=2;
     56  case 9,tblmagique=[1 2 3 4 5 6 7 8 9];         iw=3;jw=3;
     57 case 10,tblmagique=[1 2 3 4 5 6 7 8 9 10];      iw=4;jw=3;
     58 case 11,tblmagique=[1 2 3 4 5 6 7 8 9 10 11];   iw=4;jw=3;
     59 case 12,tblmagique=[1 2 3 4 5 6 7 8 9 10 11 12];iw=4;jw=3;
     60end
     61
     62
     63
     64end %if
     65
     66
     67
     68f1=figur;clf; hold on
     69
     70for iN=1:NN
     71     subplot(iw,jw,tblmagique(iN)); hold on
     72C = squeeze(real(CHP(iN,:,:)));
     73%plotm(C./xtrm(C),typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     74plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     75if(iN==1);cx=caxis;end;caxis(cx);
     76colormap(mycolormap(mapanom,21));
     77
     78%c=colorbar(CBAROR);
     79titre1=strcat('EOF',num2str(iN));
     80titre2=strcat(' (',num2str(expvar(iN)),'%)');
     81title(strcat(titre1,titre2));
     82
     83end %for iN
     84
     85%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
     86titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%');
     87%titre3=strcat('(Analysed file: ',ffile,')');
     88%titre=strvcat(titre2,titre3);
     89suptitle(titre2);
     90
     91set(f1,'Position',[378 39 313 647]); % Laptop screen
     92%set(f1,'Position',[369 55 316 899]); % Desktop screen
     93%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART'));
     94set(f1,'Name','<> REAL PART <>');
     95
     96
     97% -------------------------------------------
     98if isreal(CHP)==0
     99
     100f2=figur;clf; hold on
     101
     102for iN=1:NN
     103     subplot(iw,jw,tblmagique(iN)); hold on
     104
     105C = squeeze(imag(CHP(iN,:,:)));
     106plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     107%if(iN==1);cx=caxis;end;
     108caxis(cx);
     109colormap(mycolormap(mapanom,21));
     110
     111c=colorbar(CBAROR);
     112titre1=strcat('EOF',num2str(iN));
     113titre2=strcat(' (',num2str(expvar(iN)),'%)');
     114title(strcat(titre1,titre2));
     115
     116end %for iNN
     117
     118%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
     119titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%');
     120%titre3=strcat('(Analysed file: ',ffile,')');
     121%titre=strvcat(titre2,titre3);
     122suptitle(titre2);
     123
     124set(f2,'Position',[378 39 313 647]); % Laptop screen
     125%set(f2,'Position',[369 55 316 899]); % Desktop screen
     126set(f2,'Name','<> IMAGINARY PART <>');
     127
     128end
     129
     130

  • ../trunk-jpl/externalpackages/pcatool/install/plotm.m

     
     1%PLOTM  Trace un champ de donnee geographique
     2% [CS,H]=PLOTM(C,[ZONE])
     3%   Trace le champ C(LONGITUDE,LATITUDE) sur une carte
     4%   du monde en projection Mercator par defaut.
     5%   Le parametre optionnel ZONE peut prendre les valeurs
     6%     ZONE = 1 : Projection Stereographique centree sur le pole sud
     7%     ZONE = 2 : Projection Mercator de l'hemisphere sud
     8%     ZONE = 3 : Projection Mercator du globe
     9%     ZONE = 4 : Pas de projection (imagesc)
     10%     ZONE = 5 : Projection Mercator de l'hemisphere sud complet
     11%     ZONE = 6 : Projection Mercator du globe complet
     12%     ZONE = 7 : Projection Stereographique centree sur le pole nord
     13%     ZONE = 8 : Idem 2 mais sans label sur les axes
     14%     ZONE = 9 : Projection Mercator de l'hemisphere sud aux limites du Fmask
     15%   Returns contour matrix CS as described in CONTOURC and a vector H
     16%   of handles to PATCH objects (for use by CLABEL).
     17%
     18%   La definition des grilles (longitude,latitude) est celle du
     19%   modele par defaut, mais elle peut etre modifiee en passant en
     20%   global les nouvelles variables nommees lon et lat.
     21%
     22% NECESSITE LE PACKAGE M_MAP
     23%
     24%
     25% (nov. 2003 -> gmaze@univ-brest.fr)
     26
     27function varargout = plotm(C,varargin)
     28
     29if (nargin>2)|(nargin<1)
     30     help plotm.m
     31     error('plotm.m : Wrong number or bad parameter(s)')
     32     return
     33elseif (nargin==2)
     34   arg = varargin{:};
     35   if (arg>14)|(arg<1)
     36     help plotm.m
     37     error('plotm.m : Wrong number or bad parameter(s)')
     38     return
     39   end
     40end %if
     41
     42%--------------------------------------------------
     43% Variabes
     44%--------------------------------------------------
     45[lx,ly]=size(C);
     46% Si lon et lat ne sont pas definie comme dans le modele
     47% on prend en global n'importe quel definition
     48if ((lx~=64)|(ly~=32))
     49   global lon lat
     50else
     51   lon= linspace(0.,360.,64);
     52   lat= [85.76 80.27 74.74 69.21 63.68 58.14 52.61 47.07 41.53 ...
     53      36.00 30.46 24.92 19.38 13.84 8.31 2.77 -2.77 -8.31 ...
     54      -13.84 -19.38 -24.92 -30.46 -36.00 -41.53 -47.07 -52.61 ...
     55      -58.14 -63.68 -69.21 -74.74 -80.27 -85.76];
     56end
     57
     58% Le nombre de contour est par defaut a 10, mais on peut le
     59% modifier en declarant nbcont en global
     60global nbcont
     61if isempty(nbcont)==1
     62   nbcont=10;
     63end
     64
     65%-----------------------------------------------
     66% Projection
     67%-----------------------------------------------
     68if (nargin==2)
     69   arg = varargin{:};
     70   switch arg(1)
     71     case 1 % Projection stereo
     72       m_proj('stereo','lon',0,'lat',-90,'rad',71);
     73     case 2 % Projection mercator (HS)
     74       m_proj('miller','lon',[0 360],'lat',[-80 -20]);
     75       m_proj('miller','lon',[0 360],'lat',[-60 -20]);
     76     case 3 % Projection mercator (Tout)  -> Defaut
     77       m_proj('miller','lon',[0 360],'lat',[-80 80]);
     78     case 5 % Projection mercator (HS-complet)
     79       m_proj('miller','lon',[0 360],'lat',[-85 -5]);
     80     case 6 % Projection mercator (Tout-complet)
     81       m_proj('miller','lon',[0 360],'lat',[-85 85]);
     82     case 7 % Projection stereo
     83       m_proj('stereo','lon',0,'lat',90,'rad',71);
     84     case 8 % Projection mercator (HS) sans Label sur les axes
     85       m_proj('miller','lon',[0 360],'lat',[-80 -20]);
     86       m_proj('miller','lon',[0 360],'lat',[-75 -30]);
     87     case 9 % Projection mercator (HS) sans Label sur les axes aux limtes du FMASK
     88       m_proj('miller','lon',[0 360],'lat',[-69 -25]);
     89     case 10 % Projection stereo recentree sur le pole sud
     90       m_proj('stereo','lon',0,'lat',-90,'rad',51);
     91     case 11 % Indien
     92       m_proj('lambert','long',[11 125],'lat',[-60 -30]);
     93     case 12 % Pacifique
     94       m_proj('lambert','long',[120 290],'lat',[-60 -30]);
     95     case 13 % Atlantique
     96       m_proj('lambert','long',[-70 20],'lat',[-60 -30]);
     97     case 14 % Centré Atlantique       
     98       m_proj('lambert','long',[-90 90],'lat',[-70 -20],'rectbox','on');
     99   end
     100else
     101% Projection mercator (Tout)  -> Defaut
     102  m_proj('mercator','lon',[0 360],'lat',[-80 80]);
     103end
     104
     105%------------------------------------------------
     106% Affichage
     107%------------------------------------------------
     108load MODELCOAST
     109if length(lon)==65
     110  modelcoast=[modelcoast;modelcoast(1,:)];
     111end
     112
     113hold on
     114if (nargin==2)
     115   arg = varargin{:};
     116   switch arg(1)
     117     case 1 % Projection stereo (Par defaut)
     118      [cs,h]=m_contourf(lon,lat,C',nbcont);
     119%      m_coast('color',[0 0 0]);
     120      m_contour(lon,lat,modelcoast',1,'k');
     121      m_grid('ytick',[-30 -40 -50 -60 -70],'yticklabels',[],...
     122             'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top')
     123     case 2 % Projection mercator (HS)
     124      [cs,h]=m_contourf(lon,lat,C',nbcont);
     125%      m_coast('color',[0 0 0]);
     126      m_contour(lon,lat,modelcoast',1,'k');
     127      m_grid('ytick',[-80 -70 -60 -50 -40 -30 -20]);
     128     case 3 % Projection mercator (Tout)
     129      [cs,h]=m_contourf(lon,lat,C',nbcont);
     130%      m_coast('color',[0 0 0]);
     131      m_contour(lon,lat,modelcoast',1,'k');
     132      m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]);
     133     case 4 % Pas de projection
     134      pcolor(lon,lat,C');
     135      axis tight;box on
     136%      axis([lon(1) lon(64) lat(32) lat(1)]);
     137%      gc=gca;set(gc,'YDir','reverse');
     138     case 5 % Projection mercator (Hs-complet)
     139      [cs,h]=m_contourf(lon,lat,C',nbcont);
     140%      m_coast('color',[0 0 0]);
     141      COASTYPE='color';COASTCOLOR=[0 0 0];
     142%      COASTYPE='patch';COASTCOLOR=[0 0.5 0.6];
     143%      m_coast(COASTYPE,COASTCOLOR);
     144      m_contour(lon,lat,modelcoast',1,'k');
     145      m_grid('ytick',[-80 -60 -40 -20 0]);
     146     case 6 % Projection mercator (Tout-complet)
     147      [cs,h]=m_contourf(lon,lat,C',nbcont);
     148%      m_coast('color',[0 0 0]);
     149      m_contour(lon,lat,modelcoast',1,'k');
     150      m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]);
     151     case 7 % Projection stereo
     152      [cs,h]=m_contourf(lon,lat,C',nbcont);
     153%      m_coast('color',[0 0 0]);
     154      m_contour(lon,lat,modelcoast',1,'k');
     155      m_grid('ytick',[30 40 50 60 70],'yticklabels',[],...
     156             'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top')
     157     case 8 % Projection mercator (HS) sans label sur les axes
     158      [cs,h]=m_contourf(lon,lat,C',nbcont);
     159%      m_pcolor(lon,lat,C');shading interp
     160%      m_coast('color',[0 0 0]);
     161      m_contour(lon,lat,modelcoast',1,'k');
     162      m_grid('ytick','','xtick','');
     163     case 9 % Projection mercator (HS) sans label sur les axes aux limtes du FMASK
     164      [cs,h]=m_contourf(lon,lat,C',nbcont);
     165%      m_coast('color',[0 0 0]);
     166      m_contour(lon,lat,modelcoast',1,'k');
     167      m_grid('ytick','','xtick','');
     168     case 10 % Projection stereo recentree sur le pole Sud
     169      [cs,h]=m_contourf(lon,lat,C',nbcont);
     170%      m_coast('color',[0 0 0]);
     171      m_contour(lon,lat,modelcoast',1,'k');
     172      m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],...
     173             'xtick',[0 60 120 180 -120 -60],'XaxisLocation','top')
     174     case 11 % Indien
     175      [cs,h]=m_contourf(lon,lat,C',nbcont);
     176      m_contour(lon,lat,modelcoast',1,'k');
     177      m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],...
     178             'xtick',[30 60 90 120 150 180],'XaxisLocation','top')
     179     case 12 % Pacifique
     180      [cs,h]=m_contourf(lon,lat,C',nbcont);
     181      m_contour(lon,lat,modelcoast',1,'k');
     182      m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],...
     183             'xtick',[120 150 180 210 240 270 300],'XaxisLocation','top')
     184     case 13 % Atlantique
     185      [cs,h]=m_contourf(lon,lat,C',nbcont);
     186      m_contour(lon,lat,modelcoast',1,'k');
     187      m_grid('ytick',[-40 -50 -60 -70],'yticklabels',[],...
     188             'xtick',[-60 -30 0 30 60],'XaxisLocation','top')
     189     case 14 % Centré Atlantique
     190      [cs,h]=m_contourf(lon,lat,C',nbcont);
     191      m_contour(lon,lat,modelcoast',1,'k');
     192      m_grid('yticklabels',[],'XaxisLocation','top',...
     193             'fontsize',6,...
     194             'xtick',[0 60 120 180 -120 -60],'xticklabels',[180 -120 -60 0 60 120]);
     195           
     196   end % switch
     197else
     198  %  Projection Mercator (Par defaut)
     199     [cs,h]=m_contourf(lon,lat,C',nbcont);
     200%     m_coast('color',[0 0 0]);
     201     m_contour(lon,lat,modelcoast',1,'k');
     202     m_grid('ytick',[-80 -60 -40 -20 0 20 40 60 80]);
     203end
     204
     205
     206
     207% OUTPUT VARIABLES
     208switch nargout
     209  case 1
     210   varargout(1) = {cs} ;
     211  case 2
     212   varargout(1) = {cs} ;
     213   varargout(2) = {h};
     214end

  • ../trunk-jpl/externalpackages/pcatool/install/calCeof.m

     
     1% [e,pc,expvar] = calCeof(M,N,METHOD)
     2%   Compute Complex EOF
     3
     4    function [e,pc,expvar] = calCeof(M,N,METHOD)
     5
     6
     7% Compute Hilbert transform Mh and create M as: M = M + i*Mh
     8M = hilbert(M);
     9
     10% Compute EOFs
     11 [e,pc,expvar] = caleof(M,N,METHOD);

  • ../trunk-jpl/externalpackages/pcatool/install/dispCEOF.m

     
     1% dispCEOF(CEOF,EXPVAR,PHI,MOD) Display few CEOFs.
     2%
     3% => DISPLAY FEW CEOFs.
     4% CEOF contains all the CEOFs as nCEOF*X*Y.
     5% EXPVAR is a matrix with the explained variance of each
     6%  CEOFs in %. This is just for title.
     7% PHI is phase when ploting maps.
     8% MOD contains explicitly the CEOFs to display.
     9%
     10% Rq: it's using the plotm function
     11%
     12%================================================================
     13
     14%  Guillaume MAZE - LPO/LMD - June 2004
     15%  gmaze@univ-brest.fr
     16
     17function [] = dispCEOF(CEOF,EXPVAR,PHI,MOD)
     18
     19% Number of EEOF to display     
     20NMOD = length(MOD);
     21NANG = length(PHI);
     22
     23% Open figure and first guest
     24figur;clf;hold on
     25typg=8;
     26width  = .9/NMOD;
     27height = .9/NANG;
     28dleft = (.95-width*NMOD)/4;
     29left  = dleft;
     30load mapanom2
     31
     32% Let'go :
     33
     34% We select a CEOF:
     35for imod = 1 : NMOD
     36    mod = MOD(imod);
     37    chp=squeeze(CEOF(mod,:,:));
     38    Rchp=real(chp);
     39    Ichp=-imag(chp);
     40
     41    %
     42    for ang=1:NANG
     43      C = cos(PHI(ang))*Rchp + sin(PHI(ang))*Ichp;
     44      m(ang)=abs(xtrm(C));
     45    end
     46    m=max(m);
     47   
     48    % and plot maps
     49    for ang=1:NANG
     50      C = cos(PHI(ang))*Rchp + sin(PHI(ang))*Ichp;     
     51      bottom = .95-height*ang;
     52      subplot('Position',[left bottom width height]);
     53      plotm(C,typg);
     54      m_text(10,-70,strcat('\phi=',num2str(PHI(ang)*180/pi),'^o'));
     55%      ylabel(strcat('\phi=',num2str(PHI(ang)*180/pi),'^o'),'Rotation',0);
     56      caxis([-m m]);
     57      colormap(mapanom);
     58      if(ang==1)
     59         titre=strcat('CEOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)');
     60         title(titre);
     61      end;
     62    end
     63    left = left + width + dleft;
     64
     65end

  • ../trunk-jpl/externalpackages/pcatool/install/dispSVDnb.m

     
     1% dispSVDnb(typg,CHP1,CHP2,N,expvar) Display SVDs
     2%
     3% => Display SVDs in white and black
     4% typg is plot type you want (see plotm.m).
     5% CHP contains SVDs and is of the form: CHP(iSVD,X,Y)
     6% N is number of SVD to display
     7% expvar contained the explained variance of SVDs
     8%
     9% Rq: it''s using the plotmNB function
     10%
     11%================================================================
     12
     13%  Guillaume MAZE - LPO/LMD - December 2004
     14%  gmaze@univ-brest.fr
     15
     16function [] = dispSVDnb(typg,CHP1,CHP2,N,expvar)
     17
     18global ffile1 ffile2
     19% --------------------------------------------------
     20% AFFICHAGE
     21% --------------------------------------------------
     22if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
     23   CBAROR='horiz';
     24   iw=N;jw=2;
     25else
     26   CBAROR='vert';
     27   iw=N;jw=2;
     28end
     29
     30
     31f1=figur;clf; hold on
     32for iN=1:N
     33     subplot(iw,jw,2*iN-1)
     34C = squeeze(real(CHP1(iN,:,:)));
     35%[cs,h]=plotmNB(C,typg); whos cs h
     36[cs,h]=plotmNB(C,typg);
     37%clabel(cs,h,'labelspacing',100);
     38titre1=strcat('SVD',num2str(iN),':',ffile1);
     39titre2=strcat(' (',num2str(expvar(iN)),'%)');
     40title(strcat(titre1,titre2));
     41
     42     subplot(iw,jw,2*iN)
     43C = squeeze(real(CHP2(iN,:,:)));
     44[cs,h]=plotmNB(C,typg);
     45%clabel(cs,h,'labelspacing',100);
     46titre1=strcat('SVD',num2str(iN),':',ffile2);
     47titre2=strcat(' (',num2str(expvar(iN)),'%)');
     48title(strcat(titre1,titre2));
     49
     50end %for iN
     51
     52
     53
     54
     55%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
     56titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%');
     57%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')');
     58%titre=strvcat(titre2,titre3);
     59suptitle(titre2);
     60
     61set(f1,'Position',[378 39 313 647]); % Laptop screen
     62%set(f1,'Position',[369 55 316 899]); % Desktop screen
     63%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART'));

  • ../trunk-jpl/externalpackages/pcatool/install/Contents.m

     
     1% PRINCIPAL COMPONENTS ANALYSIS TOOLBOX
     2% Last Modified: Dec 2007
     3%
     4% This software is provided "as is" without express
     5% or implied warranty...
     6%
     7%
     8% EOF and SVD
     9%    caleof       - Compute EOFs
     10%    calEeof      - Compute Extended EOFs
     11%    calCeof      - Compute Complex  EOFs
     12%    calsvd       - Compute SVDs
     13%    confexpvar   - Compute statistical confidence of explained variance
     14%
     15% VISUALISATION
     16%    dispEOF      - Display EOFs
     17%    dispEEOF     - Display EEOFs as time sequence map
     18%    dispEEOF_v2  - Display EEOFs as time sequence map with
     19%                   the PC time series
     20%    dispCEOF     - Display CEOFs
     21%    dispSVD      - Display SVDs
     22%    plotm        - Plot a 2D geo field
     23%
     24% Matrix manipulation
     25%    map2mat      - Transform C(TIME,LON,LAT) in D(TIME,PT)
     26%                   under a geographic points selection
     27%    mat2map      - Invert transform done by map2mat
     28%
     29%
     30%================================================================
     31% Author: Guillaume MAZE - LPO/LMD - March 2004
     32%         gmaze@univ-brest.fr

  • ../trunk-jpl/externalpackages/pcatool/install/calEeof.m

     
     1% [EOFs,PC,EXPVAR] = calEeof(M,N,METHOD,NLAG,DT) Compute EEOF
     2%
     3% => Compute N Extended EOF of a M(TIME*MAP) matrix.
     4% METHOD is the EOF method computing (see caleof.m)
     5% NLAG is the number of map you want
     6% DT is the increment between each map (not the time step
     7% of your data)
     8%
     9% The function can filter datas with a Lanczos band-pass
     10% filter. Edit the calEeof.m file and change option in it.
     11%
     12% See also: CALEOF
     13%
     14% Ref: Weare and Nasstrom 1982
     15%================================================================
     16
     17%  Guillaume MAZE - LPO/LMD - March 2004
     18%  gmaze@univ-brest.fr
     19
     20    function varargout = calEeof(M,N,METHOD,NLAG,DT)
     21
     22% ---------------------------------
     23% Divers
     24tic
     25LAG = NLAG*DT;
     26foll = 1; % Turn this value to 0 if you don t want commentary along computing
     27
     28% ---------------------------------
     29% LAG must be >0, LAG=1 correspond to initial position
     30if LAG<=0
     31   LAG = 1;
     32end
     33
     34% ---------------------------------
     35% Usely for us, M is a TIME*MAP matrix so:
     36% We get dimensions
     37  [n p]=size(M);
     38% and put M in the correct form for computing
     39  M = M';
     40
     41% ---------------------------------
     42% Eventualy time filtering of data
     43if 0==1
     44   if (foll~=0),disp('==> Time filtering...');end
     45   dt = 1;         % Real time step of your data (diff from DT)
     46   Fc  = 1/dt/8;   % This the low cut off frequency
     47   Fc2 = 1/dt/2;   % and this the high cut off frequency
     48   SIGNAL = M(1,:);
     49   nj = fix(length(SIGNAL)/10); % Nb of points of the window
     50
     51   for ipt = 1 : p
     52       SIGNAL = M(ipt,:)';
     53       SIGNALF = lanczos(SIGNAL,Fc2,nj);
     54       SIGNALF = SIGNALF - lanczos(SIGNALF,Fc,nj);
     55       Y(:,ipt) = SIGNALF;
     56       if mod(ipt,10)==0 % We display a plot of filtered data
     57          clf;
     58          plot((0:n-1),SIGNAL,'k',(0:n-1),SIGNALF,'r');         
     59          drawnow;
     60       end
     61   end
     62   M = Y';
     63end, clear Fc Fc2 SIGNAL nj SIGNALF Y ipt, close
     64
     65
     66% ---------------------------------
     67% This is the matrix where we concat all submatrices from M.
     68% Size of each submatrix is : NLAG*p rows vs n-LAG+1 colums
     69F = zeros(NLAG*p,n-LAG+1);
     70
     71% ---------------------------------
     72% Let's go for F:
     73if (foll~=0),disp('==> Forming concatenated matrix...');end
     74if DT==1
     75% CLASSIC CASE
     76for ilag = 1 : LAG
     77       % Extract submatrix ilag from M
     78       Msub = M(:,ilag:n-LAG+ilag);
     79       % Concat matrix
     80       F( (ilag-1)*p+1 : ilag*p , : ) = Msub;
     81end
     82else
     83% DT>1
     84it=0;
     85for ilag = 1 : DT : LAG
     86       % Extract submatrix ilag from M
     87         Msub = M( : , ilag : n-LAG+ilag);
     88       % Concat matrix
     89         it = it + 1;
     90         F( (it-1)*p+1 : it*p , : ) = Msub;
     91         % imagesc(Msub);pause % for debugging
     92end
     93end
     94
     95% ---------------------------------
     96% Compute EOFs by normal way
     97% (Don't forget that caleof is taking a TIME*MAP argument matrix,
     98% that's why we have F' !)
     99if (foll~=0),disp('==> Computing EOFs ...');end
     100[e,pc,expvar] = caleof(F',N,METHOD);
     101
     102% ---------------------------------
     103% e is the matrix with EOFs inside. We should extract each map
     104% for different lags from it.
     105% e is NEOF*(LAG*MAP) and we construct 3D matrix CHP as NEOF*LAG*MAP
     106if (foll~=0),disp('==> Reshaping vectors ...');end
     107for ilag = 1 : NLAG
     108    E = e( : , (ilag-1)*p+1 : ilag*p );
     109    CHP(:,ilag,:) = E;
     110end
     111
     112
     113% ---------------------------------
     114% OUTPUT VARIABLES
     115switch nargout
     116  case 1
     117   varargout(1) = {CHP} ;
     118  case 2
     119   varargout(1) = {CHP} ;
     120   varargout(2) = {pc};
     121  case 3
     122   varargout(1) = {CHP} ;
     123   varargout(2) = {pc};
     124   varargout(3) = {expvar};
     125end
     126
     127% ---------------------------------
     128if (foll~=0),toc,disp('==> That''s all folks !');end

  • ../trunk-jpl/externalpackages/pcatool/install/calsvd2.m

     
     1% [SVD1,SVD2,PC1,PC2,EXPVAR,Lambda] = CALSVD2(A,B,N) Compute SVDs
     2%
     3% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD -
     4%      Analyses of climatic Data" 1997
     5%================================================================
     6%
     7%  Guillaume MAZE - LPO/LMD - March 2004
     8%  gmaze@univ-brest.fr
     9
     10
     11   function [e1,e2,pc1,pc2,expvar,Lambda,dsumCF] = calsvd2(A,B,N);
     12
     13
     14%================================================================
     15% Ref: H. Bjornson and S.A. Venegas: "A manual for EOF and SVD -
     16% Analyses of climatic Data" 1997 => p18
     17
     18% Assume that A is (time*map) matrix
     19[n p]=size(A);
     20
     21% Remove the mean of each column (ie the time mean in each station records)
     22S=detrend(A,'constant');
     23P=detrend(B,'constant');
     24
     25% Form the covariance matrix:
     26C=S'*P;
     27
     28% Find eigenvectors and singular values
     29[U,Lambda,V] = svds(C,N);
     30
     31% PC
     32a=S*U;
     33b=P*V;
     34
     35% Make them clear for output
     36for iN=1:N
     37    e1(iN,:) = squeeze( U(:,iN) )';
     38   pc1(iN,:) = squeeze( a(:,iN) )';
     39    e2(iN,:) = squeeze( V(:,iN) )';
     40   pc2(iN,:) = squeeze( b(:,iN) )';
     41end
     42
     43% Amount of variance explained a 0.1 pres et en %
     44L2=Lambda.^2;
     45dsum=diag(L2)/trace(L2);
     46for iN=1:N
     47   expvar(iN)=fix( ( dsum(iN)*100/sum(dsum) )*10 ) /10;
     48end
     49

  • ../trunk-jpl/externalpackages/pcatool/install/dispEEOF.m

     
     1% dispEEOF(CHP,EXPVAR,DT,NLAG,MOD) Display few EEOFs.
     2%
     3% => DISPLAY FEW EEOFs.
     4% CHP contains all the EEOFs as EOF*LAG*X*Y.
     5% EXPVAR is a matrix with the explained variance of each
     6%  EEOFs in %. This is just for title.
     7% DT is time step between maps.
     8% NLAG is the number of LAG to display.
     9% MOD contains explicitly the EEOFs to display.
     10%
     11% Rq: it's using the plotm function
     12%
     13%================================================================
     14
     15%  Guillaume MAZE - LPO/LMD - March 2004
     16%  gmaze@univ-brest.fr
     17
     18function [] = dispEEOF(CHP,EXPVAR,DT,NLAG,MOD)
     19
     20% Number of EEOF to display     
     21NMOD = length(MOD);
     22
     23% Open figure and first guest
     24figur;
     25clf;hold on
     26typg=8;
     27width  = .9/NMOD;
     28height = .9/NLAG;
     29dleft = (.95-width*NMOD)/4;
     30left  = dleft;
     31load mapanom
     32
     33% Let'go :
     34
     35% We choose an EEOF
     36for imod = 1 : NMOD
     37    mod = MOD(imod);
     38   
     39    % and plot maps
     40    for lag=1:NLAG
     41      bottom = .95-height*lag;
     42      subplot('Position',[left bottom width height]);
     43      C = squeeze(real(CHP(mod,lag,:,:)));
     44%      C = C./xtrm(C); % Eventually normalise field
     45      plotm(C,typg);
     46      if(lag==1),caxis([-abs(xtrm(C)) abs(xtrm(C))]);cx=caxis;end;
     47      caxis(cx);
     48      colormap(mapanom);
     49      if(lag==1)
     50         titre=strcat('EOF',num2str(mod),'(',num2str(EXPVAR(mod)),'%)');
     51         titre=strcat(titre,'; DT=',num2str(DT),' ');
     52         title(titre);
     53      end;
     54    end
     55    left = left + width + dleft;
     56
     57end

  • ../trunk-jpl/externalpackages/pcatool/install/dispSVDnbc.m

     
     1% dispSVDnb(typg,CHP1,CHP2,N,expvar) Display SVDs
     2%
     3% => Display SVDs in white and black with correlation form
     4% typg is plot type you want (see plotm.m).
     5% CHP contains SVDs and is of the form: CHP(iSVD,X,Y)
     6% N is number of SVD to display
     7% expvar contained the explained variance of SVDs
     8%
     9% Rq: it''s using the plotmNB function
     10%
     11%================================================================
     12
     13%  Guillaume MAZE - LPO/LMD - December 2004
     14%  gmaze@univ-brest.fr
     15
     16function [] = dispSVDnbc(typg,CHP1,CHP2,N,expvar)
     17
     18global ffile1 ffile2 nbcont
     19nbcont=[-1:0.2:1];
     20nbcont2=[-1:0.4:1];
     21
     22% --------------------------------------------------
     23% AFFICHAGE
     24% --------------------------------------------------
     25if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
     26   CBAROR='horiz';
     27   iw=N;jw=2;
     28else
     29   CBAROR='vert';
     30   iw=N;jw=2;
     31end
     32
     33
     34f1=figur;clf; hold on
     35for iN=1:N
     36     subplot(iw,jw,2*iN-1)
     37C = squeeze(real(CHP1(iN,:,:)));
     38[cs,h]=plotmNB(C,typg);
     39clabel(cs,h,nbcont2,'labelspacing',100,'fontsize',9);
     40titre1=strcat('SVD',num2str(iN),':',ffile1);
     41titre2=strcat(' (',num2str(expvar(iN)),'%)');
     42title(strcat(titre1,titre2));
     43
     44     subplot(iw,jw,2*iN)
     45C = squeeze(real(CHP2(iN,:,:)));
     46[cs,h]=plotmNB(C,typg);
     47clabel(cs,h,nbcont2,'labelspacing',100,'fontsize',9);
     48%clabel(cs,h,'labelspacing',100);
     49titre1=strcat('SVD',num2str(iN),':',ffile2);
     50titre2=strcat(' (',num2str(expvar(iN)),'%)');
     51title(strcat(titre1,titre2));
     52
     53end %for iN
     54
     55
     56
     57
     58%titre1=strcat(num2str(lat(il1)),'/',num2str(lat(il2)));
     59titre2=strcat(' Total explained variance:',num2str(sum(expvar)),'%');
     60%titre3=strcat('(Analysed files: ',ffile1,' vs ',ffile2,')');
     61%titre=strvcat(titre2,titre3);
     62suptitle(titre2);
     63
     64set(f1,'Position',[378 39 313 647]); % Laptop screen
     65%set(f1,'Position',[369 55 316 899]); % Desktop screen
     66%set(f1,'Name',strcat(ffile,'<>',titre1,'<> REAL PART'));

  • ../trunk-jpl/externalpackages/pcatool/install/test_pcatool.m

     
     1% Here is a simple script to test caleof.m from PCATOOL package.
     2%
     3
     4% gmaze@mit.edu
     5% 2006/07/31
     6%
     7
     8clear
     9
     10field = 1; % The best one
     11
     12% Generate the data
     13x = -10 : 1 : 10 ;
     14y = x;
     15[X,Y] = meshgrid(x,y);
     16
     17switch field
     18 case 1   % H. Bjornson and S.A. Venegas
     19   HH = 12*( 1.2 - 0.35*sqrt( (X-0).^2 + (Y-0).^2 ) );
     20   HH = HH - mean(mean(HH));
     21   H(1,:,:) = 1012 + HH;
     22   H(2,:,:) = 1012 + (1012-squeeze(H(1,:,:)));;
     23   H(3,:,:) = 1022.8 - 3.6*Y;
     24   H(4,:,:) = 1012 + (1012-squeeze(H(3,:,:)));
     25   H(5,:,:) = 1001.2 + 3.6*X;
     26   H(6,:,:) = 1012 + (1012-squeeze(H(5,:,:)));
     27  cont = [980:5:1050];
     28  N = 2;
     29 
     30 case 11   % H. Bjornson and S.A. Venegas modified
     31   % Here it must be only 2 non-zero EOF (45degre oriented ) of similar variance 50%
     32   H = ones(4,length(x),length(y));
     33   H(1,:,:) =  Y;
     34   H(2,:,:) = -Y;
     35   H(3,:,:) = X;
     36   H(4,:,:) = -X;
     37   N = 2;
     38
     39 case 2 % Hartmann eg:  analysis in place
     40  a = [2 4 -6 8 ; 1 2 -3 4];
     41  H(1,:,:) = a(1,:)';
     42  H(2,:,:) = a(2,:)';
     43 
     44  % Get EOFs:
     45  N = 4;
     46  for method = 2 : -1 : 1
     47    [E,pc,expvar] = caleof(H,N,method);
     48    %E'
     49  end
     50
     51  return
     52 
     53 case 3  % My signal
     54  np = 4; % Numer of random signal (...and eof)
     55  nt = 10; % Number of timestep
     56  H = zeros(nt,length(x),length(y));
     57  for ip = 1 : np
     58    xc(ip) = abs(fix(rand(1)*10));
     59    yc(ip) = abs(fix(rand(1)*10));
     60    if ip>=fix(np/2)
     61      xc(ip) = -xc(ip);
     62      yc(ip) = -yc(ip);
     63    end
     64    dd(ip) = fix(rand(1)*10);
     65  end %for ip
     66  f = 1/nt;
     67  for it = 1 : nt
     68    H2 = zeros(length(x),length(y));
     69    for ip = 1 : ip
     70        if it==1,[xc(ip) yc(ip) dd(ip)],end
     71        HH = 12*( 1.2 - 0.35*sqrt(  ((X-xc(ip)).^2 + (Y-yc(ip)).^2 )/dd(ip) ));
     72        H2 = HH - mean(mean(HH));
     73        H(it,:,:) = squeeze(H(it,:,:)) + H2.*cos(pi*it/dd(ip)) ;
     74    end %for ip
     75  end %for it
     76  H = 1012 + H;
     77  cont = [980:5:1050];
     78  N = 3;
     79 
     80 case 4 % My signal 2
     81  x = -pi : pi/6 : pi ;
     82  y = x;
     83  [X,Y] = meshgrid(x,y);
     84  HH  = cos(X) + cos(Y);
     85  HH2 = cos(Y);
     86  nt = 12;
     87  for it = 1 : nt
     88%     H(it,:,:) = cos(pi*it/nt)*HH; cont=[-2 2];
     89     H(it,:,:) = 2*cos(pi*it/nt)*HH + 3*cos(pi*it/nt/2)*HH2; cont=[-10 10];
     90     xtrm(squeeze(H(it,:,:)));
     91  end
     92  cont=[-2 2];
     93  N = 3;
     94 
     95end % switch field
     96
     97%return
     98
     99% Plot field time serie:
     100if 1
     101figure;iw=2;jw=size(H,1)/iw;
     102set(gcf,'position',[11 533 560 420]);
     103
     104for i=1:iw*jw
     105  C = squeeze(H(i,:,:));
     106 
     107  subplot(iw,jw,i);
     108  pcolor(X,Y,C);
     109  title(num2str(i));
     110  if i==1,cx=caxis;end
     111  axis square
     112  caxis(cx);
     113
     114end %for i
     115end %fi
     116%return
     117
     118% Get EOFs:
     119G = map2mat(ones(size(H,2),size(H,3)),H);
     120
     121for method = 1 : 4
     122  [E,pc,expvar] = caleof(G,N,method); 
     123  eof = mat2map(ones(size(H,2),size(H,3)),E);
     124 
     125
     126figure;iw=1;jw=N+1;
     127set(gcf,'MenuBar','none');
     128posi = [576 0 710 205];
     129set(gcf,'position',[posi(1) (method-1)*250+50 posi(3) posi(4)]);
     130
     131for i=1:iw*jw
     132  if i<= iw*jw-1
     133   
     134  C = squeeze(eof(i,:,:));
     135  subplot(iw,jw,i);
     136  cont = 12;
     137  [cs,h] = contourf(X,Y,C,cont);
     138  clabel(cs,h);
     139  title(strcat('EOF:',num2str(i),'/',num2str(expvar(i)),'%'));
     140  axis square;
     141  %caxis([cont(1) cont(end)]);
     142 
     143  else
     144  subplot(iw,jw,iw*jw);
     145  plot(pc');
     146  grid on
     147  xlabel('time')
     148  title('PC')
     149  legend(num2str([1:N]'),2);
     150  box on
     151 
     152
     153  end %if
     154 
     155end %for i
     156%suptitle(strcat('METHODE:',num2str(method)));
     157
     158end %for method
     159
     160

  • ../trunk-jpl/externalpackages/pcatool/install/dispSVDlag.m

     
     1% dispSVDlag(typg,CHP1,CHP2,N,expvar,LAGS,dt) Display SVDs lagged
     2%
     3% => Display SVDs lagged
     4% typg is plot type you want (see plotm.m).
     5% CHP contains SVDs and is of the form: CHP(iSVD,X,Y)
     6% N is number of SVD to display
     7% expvar contained the explained variance of SVDs
     8%
     9% Rq: it''s using the plotm function
     10%
     11%================================================================
     12
     13%  Guillaume MAZE - LPO/LMD - December 2004
     14%  gmaze@univ-brest.fr
     15
     16   function [] = dispSVD(typg,CHP1,CHP2,N,expvar,LAGS,dt)
     17
     18load mapanom2
     19global ffile1 ffile2
     20% --------------------------------------------------
     21% AFFICHAGE
     22% --------------------------------------------------
     23if (typg~=1)&(typg~=7)&(typg~=3)&(typg~=6)
     24   CBAROR='horiz';
     25else
     26   CBAROR='vert';
     27end
     28
     29   CHP1=squeeze(CHP1(:,N,:,:));
     30   CHP2=squeeze(CHP2(:,N,:,:));
     31
     32Nlag=size(CHP1,1);
     33
     34iw=round(Nlag/dt);jw=2;
     35
     36
     37f1=figur;clf; hold on
     38ilag=0;
     39for lag=1:dt:Nlag
     40  ilag=ilag+1;
     41
     42    subplot(iw,jw,2*ilag-1)
     43C = squeeze(real(CHP1(lag,:,:)));
     44plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     45if(ilag==1);cx=caxis;end;caxis(cx);
     46colormap(mycolormap(mapanom,21));
     47c=colorbar(CBAROR);
     48titre1=strcat('lag:',num2str(LAGS(lag)),'y: ',ffile1);
     49title(titre1);
     50
     51     subplot(iw,jw,2*ilag)
     52C = squeeze(real(CHP2(lag,:,:)));
     53plotm(C,typg);caxis([-abs(xtrm(C)) abs(xtrm(C))]);
     54if(ilag==1);cx=caxis;end;caxis(cx);
     55colormap(mycolormap(mapanom,21));
     56c=colorbar(CBAROR);
     57titre1=strcat('lag:',num2str(LAGS(lag)),'y: ',ffile2);
     58title(titre1);
     59
     60end %for iN
     61
     62suptitle(strcat('Lagged SVD, mode:',num2str(N)));

  • ../trunk-jpl/externalpackages/pcatool/install/map2mat.m

     
     1% D = MAP2MAT(F,C) Reshaping matrix
     2%
     3% => Conversion of a 'map' matrix C(TIME,LON,LAT) into a D(TIME,PT) matrix
     4% under the mask F(LON,LAT).
     5% F is a matrix contenaing 1 where you would like to keep the point and
     6%  0 elsewhere.
     7%
     8% Rq: No check is done about the input.
     9%
     10% See also: mat2map
     11%================================================================
     12
     13% March 2004
     14% gmaze@univ-brest.fr
     15
     16function [D] = map2mat(F,C);
     17
     18% Get dimensions
     19[tps nolon nolat] = size(C);
     20
     21% So output matrix will be:
     22D = zeros(tps,nolon*nolat);
     23
     24% point indice
     25ipt = 0;
     26
     27% 'Un-mapping' :
     28for iy=1:nolat
     29  for ix=1:nolon
     30      if F(ix,iy)>0
     31         ipt = ipt + 1;
     32         D(:,ipt)=squeeze(C(:,ix,iy));
     33      end % if
     34  end % for
     35end %for
     36
     37% OUTPUT:
     38D = squeeze(D(:,1:ipt));
     39
     40

  • ../trunk-jpl/externalpackages/pcatool/install/mat2map.m

     
     1% C = MAP2MAT(F,D) Reshaping matrix
     2%
     3% => Creation of a 'map' matrix C(TIME,LON,LAT) from D(TIME,PT)
     4% following mask F(LON,LAT).
     5% F is a matrix contenaing 1 where you would like to keep the point and
     6%  0 elsewhere (see mat2map).
     7%
     8% Rq: No check is done about the input.
     9%
     10% See also: map2mat
     11%================================================================
     12
     13% March 2004
     14% gmaze@univ-brest.fr
     15
     16function [C] = mat2map(F,D);
     17
     18% Get dimensions
     19[nolon nolat] = size(F);
     20[time npt] = size(D);
     21
     22% So output 'map' matrix has the form:
     23C = zeros(time,nolon,nolat);
     24
     25% Variables
     26nul = NaN.*ones(time,1);
     27ipt = 0 ;
     28
     29
     30% 'mapping' :
     31  for iy=1:nolat
     32      for ix=1:nolon
     33         if F(ix,iy)>0
     34            ipt = ipt + 1;
     35            C(:,ix,iy) = D(:,ipt);
     36         else
     37            C(:,ix,iy) = nul;
     38         end %if
     39      end %for ix
     40  end %for iy
     41
Note: See TracBrowser for help on using the repository browser.