How to use Ferret External Functions for computing EOF's

Below is an annotated version of the script ef_eof_demo.jnl

! ef_eof_demo.jnl (8/2000)

! Description: Demonstration of computing EOFs using the
! External Functions EOF_SPACE, EOF_TFUNC, EOF_STAT

The functions implements Chelton's '82 method for finding EOFs of gappy time series. All the functions perform the same computation but return different portions of the results. They take time to run as the computation inverts a matrix of size [NT by (NX*NY)]

For all functions the arguments are as follows:

!
! EOF functions Example 1: Define a function of x,y,time using trig
! functions.  Decompose into spatial and time EOF's and display statistics.
! 
yes? LET time = t[t=1-jan-1990:10-jan-1990:24]   ! 24 hour resolution time axis
 
yes? DEFINE AXIS/x=0:10:0.5 x10
yes? DEFINE AXIS/y=0:10:0.5 y10
yes? DEFINE GRID/x=x10/y=y10/t=time g10x10
yes? SET GRID g10x10
 
yes? LET fcn1 = 15.* sin(omega1*t)*cos(r)/(r+1)
yes? LET fcn2 = 20.* sin(omega2*t)*(sin(s)-.2*sin(q))/(s+1)
 
yes? LET r = ((6.*(xpts-x0)^2 + 7.*(ypts-y0)^2)^0.5)
yes? LET s = (((xpts-x1)^2 + 2*(ypts-y1)^2)^0.5)
yes? LET q = ((3*(xpts-x0)^2 + (ypts-y1)^2)^0.5)
 
yes? LET x0 = 2
yes? LET y0 = 4
yes? LET x1 = 5
yes? LET y1 = 7
 
yes? LET omega1 = 1/10*2*3.14159
yes? LET omega2 = 2 * omega1
 
yes? LET sample_function = fcn1 + fcn2
yes? LET xpts = x
yes? LET ypts = y
 

!  Compute the statistics on the EOF solution.  EOF_STAT returns:
!  for J=1  the number of EOFs returned
!  for J=2  the percent variance explained by each EOF
!  for J=3  the eigenvalue for each EOF


yes? LET estat = eof_stat(sample_function, 0.1)
 
yes? LIST/I=1/J=1 estat
             EOF_STAT(SAMPLE_FUNCTION, 0.1)
             X: 1
             Y: 1
          2.000
 
yes? LIST/I=1:4/J=2 estat
             EOF_STAT(SAMPLE_FUNCTION, 0.1)
             2    
 1   / 1:  93.59
 2   / 2:   6.41
 3   / 3:   0.00
 4   / 4:   0.00
 
yes? LIST/I=1:4/J=3 estat
             EOF_STAT(SAMPLE_FUNCTION, 0.1)
             3   
 1   / 1:  1104.
 2   / 2:    76.
 3   / 3:     0.
 4   / 4:     0.
 
!  Plot the original function (averaged over time) and its spatial 
!  decomposition by EOF_SPACE
 
yes? DEFINE VIEW/xlim=0.,.33/ylim=.6,1./text=0.2 vul
yes? DEFINE VIEW/xlim=.33,.66/ylim=.6,1./text=0.2 vuc
yes? DEFINE VIEW/xlim=.66,1./ylim=.6,1./text=0.2 vur
 
yes? DEFINE VIEW/xlim=0.,.33/ylim=.1,.5/text=0.2 vll
yes? DEFINE VIEW/xlim=.33,.66/ylim=.1,.5/text=0.2 vlc
yes? DEFINE VIEW/xlim=.66,1./ylim=.1,.5/text=0.2 vlr
 
yes? SET VIEW vul; CONTOUR/TITLE="FCN1" fcn1[l=1:10@ave]
yes? SET VIEW vuc; CONTOUR/TITLE="FCN2" fcn2[l=1:10@ave]
yes? SET VIEW vur; CONTOUR/TITLE="FCN1 + FCN2" sample_function[l=1:10@ave]
 

yes? LET exy = eof_space(sample_function, 0.1)
 
yes? SET VIEW vll; CONTOUR/L=1/TITLE="EOF 1" exy
yes? SET VIEW vlc; CONTOUR/L=2/TITLE="EOF 2" exy




 
yes? CANCEL VIEW


! Now plot the time amplitude functions.
 
yes? LET etim = eof_tfunc(sample_function, 0.1)
 
yes? SET VIEW ul
yes? PLOT/I=1/TITLE=taf1 etim
 
yes? SET VIEW ur
yes? PLOT/I=2/TITLE=taf2 etim
 
! should be all bad flags...
yes? SET VIEW ll
yes? PLOT/I=3/TITLE=taf3 etim





 
yes? CANCEL DATA/ALL
yes? CANCEL VARIABLE/ALL
yes? CANCEL REGION

!  A second example, using the COADS climatology data.
!

yes? USE coads_climatology
yes? SET REGION/X=67w:1w/Y=11S:11N

!  Compute and save the spatial EOF functions.  These have the same units
!  as the data.
 
yes? LET eof_xyfcn = eof_space(sst, 0.5)
yes? SAVE/CLOBBER/FILE=sst_clim_eof_space.cdf eof_xyfcn
yes? LISTing to file sst_clim_eof_space.cdf
 
yes? CANCEL DATA/ALL
yes? CANCEL VARIABLE/ALL
yes? CANCEL REGION

yes? USE sst_clim_eof_space.cdf
yes? SET VIEW ul; fill/l=1/TITLE="eof 1" eof_xyfcn; go land
yes? SET VIEW ur; fill/l=2/TITLE="eof 2" eof_xyfcn; go land
yes? SET VIEW ll; fill/l=3/TITLE="eof 3" eof_xyfcn; go land
 
 


yes? CANCEL DATA/ALL
yes? CANCEL VARIABLE/ALL
yes? CANCEL REGION

!  Compute the statistics on the EOFs: number of EOFs scaled
!  and returned; Percent variance explained by each EOF; and the
!  eigenvalues.
 

yes? USE coads_climatology
yes? SET REGION/X=67w:1w/Y=11S:11N
yes? LET eofstat = eof_stat(sst[X=67w:1w,Y=11S:11N], 0.5)

yes? let nout = eofstat[i=1,j=1]
yes? let pcts = eofstat[i=1:`nout`,j=2]
 !-> DEFINE VARIABLE pcts = eofstat[i=1:3,j=2]
yes? let eigenv = eofstat[i=1:`nout`,j=3]
 !-> DEFINE VARIABLE eigenv = eofstat[i=1:3,j=3]


yes? LIST nout
             VARIABLE : EOFSTAT[I=1,J=1]
             DATA SET : COADS Monthly Climatology (1946-1989)
             FILENAME : coads_climatology.des
             FILEPATH : /home/ja9/tmap/fer_dsets/descr/
             X        : 1
             Y        : 1
          285.0

yes? LIST/I=1:6 pcts
             VARIABLE : EOFSTAT[I=1:285,J=2]
             DATA SET : COADS Monthly Climatology (1946-1989)
             FILENAME : coads_climatology.des
             FILEPATH : /home/ja9/tmap/fer_dsets/descr/
             SUBSET   : 6 points (X)
             Y        : 2
               2    
               2
 1    /  1:  86.62
 2    /  2:   5.97
 3    /  3:   3.89
 4    /  4:   1.53
 5    /  5:   0.67
 6    /  6:   0.38

yes? LIST/I=1:6 eigenv
             VARIABLE : EOFSTAT[I=1:285,J=3]
             DATA SET : COADS Monthly Climatology (1946-1989)
             FILENAME : coads_climatology.des
             FILEPATH : /home/ja9/tmap/fer_dsets/descr/
             SUBSET   : 6 points (X)
             Y        : 3
               3    
               3
 1    /  1:  249.4
 2    /  2:   17.2
 3    /  3:   11.2
 4    /  4:    4.4
 5    /  5:    1.9
 6    /  6:    1.1

 yes? CANCEL REGION
 
!  Compute and save time amplitude functions: Note they are dimensionless.
 
yes? USE coads_climatology
yes? LET eoftime = eof_tfunc(sst[X=67w:1w,Y=11S:11N], 0.5)
 
SAVE/CLOBBER/FILE=sst_clim_eof_tfunc.cdf eoftime[i=1:`nout`]
 !-> LIST/FORMAT=CDF/CLOBBER/FILE=sst_clim_eof_tfunc.cdf eoftime[i=1:3]
 LISTing to file sst_clim_eof_tfunc.cdf
 
yes? CANCEL VARIABLE/ALL
yes? CANCEL VIEW
yes? USE sst_clim_eof_tfunc.cdf

yes? SET VIEW ul
yes? PLOT/I=1/TITLE="time function 1"/VLIMITS=-2:2:0.5 eoftime
yes? SET VIEW ur
yes? PLOT/I=2/TITLE="time function 2"/VLIMITS=-2:2:0.5 eoftime
yes? SET VIEW ll
yes? PLOT/I=3/TITLE="time function 3"/VLIMITS=-2:2:0.5 eoftime



 


oar.pmel.contact_ferret@noaa.gov 30 Apr, 1999