;
; For a test data set, computes EOFs and PCs, then tries:
; (i) to reproduce the PCs by projecting the data onto the EOFs, and
; (ii) to reproduce the original dataset by summing the product of the
; EOFs and the PCs over the leading N EOFs
;
nretain=15
;
; Restore monthly temperature gridded dataset
;
print,'Reading in temperature data'
restore,filename='obs_temp_mon.idlsave'
;  timey,fdmon,xlon,ylat,nyr,nx,ny,fdltm,fdsd,missfrac
;
; Select the spatiotemporal subset to analyse
;
; Pick period first
ksub=where((timey ge 1901) and (timey le 1930),nsub)
fdsub=reform(fdmon(*,*,*,ksub),nx*ny,nmon*nsub)
; Now extract only those boxes with full data
fdnum=total(finite(fdsub),2)
kl=where(fdnum eq nmon*nsub,nbox)
fdsub=fdsub(kl,*)
nval=nmon*nsub
print,'Time values = ',nval
print,'Space values = ',nbox
; Now convert to anomalies
mkanomaly,fdsub
;
; Compute the EOFs
;
print,'Computing EOFs'
myeof1d,fdsub,ev,ea,lam,lampct,lamcum,$
  nretain=nretain,correlation=1
;
; Project data onto EOFs to estimate the PCs
;
mknormal,fdsub
allea=fltarr(nval,nretain)
for iretain = 0 , nretain-1 do begin
  print,iretain,format='($,I4)'
  for i = 0 , nval-1 do begin
    allea(i,iretain)=total(fdsub(*,i)*ev(*,iretain),/nan)
  endfor
endfor
print
;
; Now plot them
;
loadct,39
multi_plot,nrow=3,ncol=1,layout='large'
if !d.name eq 'X' then window,ysize=850
;
lam=[!values.f_nan,lam]
lampct=[!values.f_nan,lampct]
lamcum=[0,lamcum]
xt=timey
plot,lam,/xstyle,xtitle='Mode',ytitle='Eigenvalue',$
  psym=10,xrange=[0,20],/ylog
plot,lampct,/xstyle,xtitle='Mode',ytitle='% variance explained',$
  psym=10,xrange=[0,20]
plot,lamcum,/xstyle,xtitle='Mode',ytitle='Cumulative % variance explained',$
  psym=10,xrange=[0,20]
;
map=def_map(/npolar)  &  map.limit(0)=0.
coast=def_coast(/get_device)
labels=def_labels(/off)
;
levels=findgen(21)*0.02-0.2
levels(0)=-0.3
levels(20)=0.3
c_colors=indgen(20)+50
def_1color,cr,cg,cb,50,color='dblue'
def_1color,cr,cg,cb,59,color='lgreen'
def_1color,cr,cg,cb,60,color='lyellow'
def_1color,cr,cg,cb,69,color='dred'
def_smearcolor,fromto=[50,59]
def_smearcolor,fromto=[60,69]
;
th=10.
;
for i = 0 , nretain-1 do begin
  pause
  ;
  yt=ea(*,i)
  yt2=allea(*,i)
  filter_cru,th,tsin=yt,tslow=tslow,/nan
  plot,yt,title='Mode'+string(i+1,format='(I2)'),$
    /xstyle,$
    xtitle='Year',ytitle='Amplitude'
  ;
  oplot,yt2,color=65
  filter_cru,th,tsin=yt,tslow=tslow,/nan
;  oplot,xt,tslow,thick=5,linestyle=1
;  tslow(ml)=!values.f_nan
;  oplot,xt,tslow,thick=5
  ;
  oplot,!x.crange,[0,0],linestyle=1
  ;
;  fd=reform(ev(*,*,i))
;  inter_boxfd,fd,xlon,ylat,$
;    map=map,coast=coast,labels=labels,$
;    levels=levels,c_colors=c_colors
endfor
;
end