\SET MODE VERIFY ! ef_eof_demo.jnl ( acm 9/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. ! * EOF_SPACE returns EOF eigenvectors: spatial EOF in x and y with the ! same units as the incoming data ! * EOF_TFUNC returns EOF time amplitude functions; dimensionless ! * EOF_STAT returns statistics on the EOF computation: number of EOFs ! scaled and returned; Percent variance explained by each EOF; and the ! eigenvalues. ! ! For all functions the arguments are as follows: ! ! * 1st argument: Input data field, a function of x, y, and time; may be a ! function of z ! * 2nd argument: The minimum percent variance explained by the EOF's that ! are computed and scaled, e.g. to return EOFs that explain at least 2% ! variance, use 2.0 ! * 3nd argument: The fraction of each time series that must be present to ! include it in the calculations, e.g. to use all time series that have ! at least half the data present, use 0.5 ! ! ! EOF functions Example 1: define a function of x,y,time using trig ! functions. Decompose into spatial and time EOF's and display statistics. ! PAUSE LET time = t[t=1-jan-1990:10-jan-1990:24] ! 24 hour resolution time axis DEFINE AXIS/x=0:10:0.5 x10 DEFINE AXIS/y=0:10:0.5 y10 DEFINE GRID/x=x10/y=y10/t=time g10x10 SET GRID g10x10 LET fcn1 = 15.* sin(omega1*t)*cos(r)/(r+1) LET fcn2 = 20.* sin(omega2*t)*(sin(s)-.2*sin(q))/(s+1) LET r = ((6.*(xpts-x0)^2 + 7.*(ypts-y0)^2)^0.5) LET s = (((xpts-x1)^2 + 2*(ypts-y1)^2)^0.5) LET q = ((3*(xpts-x0)^2 + (ypts-y1)^2)^0.5) LET x0 = 2 LET y0 = 4 LET x1 = 5 LET y1 = 7 LET omega1 = 1/10*2*3.14159 LET omega2 = 2 * omega1 LET sample_function = fcn1 + fcn2 LET xpts = x 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 LET estat = eof_stat(sample_function, 0.1) LIST/I=1/J=1 estat LIST/I=1:4/J=2 estat LIST/I=1:4/J=3 estat ! Plot the original function (averaged over time) and its spatial ! decomposition by EOF_SPACE PAUSE DEFINE VIEW/xlim=0.,.33/ylim=.6,1./text=0.2 vul DEFINE VIEW/xlim=.33,.66/ylim=.6,1./text=0.2 vuc DEFINE VIEW/xlim=.66,1./ylim=.6,1./text=0.2 vur DEFINE VIEW/xlim=0.,.33/ylim=.1,.5/text=0.2 vll DEFINE VIEW/xlim=.33,.66/ylim=.1,.5/text=0.2 vlc DEFINE VIEW/xlim=.66,1./ylim=.1,.5/text=0.2 vlr SET VIEW vul; CONTOUR/TITLE="FCN1" fcn1[l=1:10@ave] SET VIEW vuc; CONTOUR/TITLE="FCN2" fcn2[l=1:10@ave] SET VIEW vur; CONTOUR/TITLE="FCN1 + FCN2" sample_function[l=1:10@ave] LET exy = eof_space(sample_function, 0.1) SET VIEW vll; CONTOUR/L=1/TITLE="EOF 1" exy SET VIEW vlc; CONTOUR/L=2/TITLE="EOF 2" exy PAUSE CANCEL VIEW ! Now plot the time amplitude functions. LET etim = eof_tfunc(sample_function, 0.1) SET VIEW ul PLOT/I=1/TITLE=taf1 etim SET VIEW ur PLOT/I=2/TITLE=taf2 etim ! should be all bad flags... SET VIEW ll PLOT/I=3/TITLE=taf3 etim PAUSE CANCEL DATA/ALL CANCEL VARIABLE/ALL CANCEL REGION CANCEL VIEW ! A second example, using the COADS climatology data. ! USE coads_climatology SET REGION/X=67w:1w/Y=11S:11N ! Compute and save the spatial EOF functions. These have the same units ! as the data. LET eof_xyfcn = eof_space(sst, 0.5) SAVE/CLOBBER/FILE=sst_clim_eof_space.cdf eof_xyfcn CANCEL DATA/ALL CANCEL VARIABLE/ALL CANCEL REGION USE sst_clim_eof_space.cdf SET VIEW ul; fill/l=1/TITLE="eof 1" eof_xyfcn; go land SET VIEW ur; fill/l=2/TITLE="eof 2" eof_xyfcn; go land SET VIEW ll; fill/l=3/TITLE="eof 3" eof_xyfcn; go land PAUSE CANCEL DATA/ALL CANCEL VARIABLE/ALL CANCEL REGION ! Compute the statistics on the EOFs: number of EOFs scaled ! and returned; Percent variance explained by each EOF; and the ! eigenvalues. USE coads_climatology SET REGION/X=67w:1w/Y=11S:11N LET eofstat = eof_stat(sst[X=67w:1w,Y=11S:11N], 0.5) LET nout = eofstat[i=1,j=1] LET pcts = eofstat[i=1:`nout`,j=2] LET eigenv = eofstat[i=1:`nout`,j=3] LIST/I=1:10 nout LIST/I=1:10 pcts LIST/I=1:10 eigenv PAUSE ! Compute and save time amplitude functions: Note they are dimensionless. USE coads_climatology 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`] CANCEL VIEW CANCEL VARIABLE/ALL USE sst_clim_eof_tfunc.cdf SET VIEW ul PLOT/I=1/TITLE="time function 1"/VLIMITS=-2:2:0.5 eoftime SET VIEW ur PLOT/I=2/TITLE="time function 2"/VLIMITS=-2:2:0.5 eoftime SET VIEW ll PLOT/I=3/TITLE="time function 3"/VLIMITS=-2:2:0.5 eoftime \SET MODE/LAST VERIFY