;
; Computes EOFs of infilled calibrated MXD gridded dataset.
; Can use corrected or uncorrected MXD data (i.e., corrected for the decline).
; Do not usually rotate, since this loses the common volcanic and global
; warming signal, and results in regional-mean series instead.
; Generally use the correlation matrix EOFs.
;
usefilt=0       ; 0=unfiltered, 1=pre-high-pass-filtered, -1=pre-low-pass-filt
thalf=10.       ; cutoff for high/low-pass filter
userot=0        ; 0=unrotated,   1=rotated
usecorr=1       ; 0=covariance 1=correlation
usedecline=1    ; 0=uncorrected 1=decline-corrected
useper=[1400,1976]
nretain=12
;
; Restore MXD gridded dataset
;
print,'Reading in MXD data'
restore,filename='calibmxd5_abdlow.idlsave'
;  g,mxdyear,mxdnyr,fdcalibu,fdcalibc,mxdfd2,timey,fdseas
if usedecline eq 1 then fdcalibu=fdcalibc
;
; Extract required period and locate boxes with complete data
;
print,'Finding spatial and temporal coverage'
kmxd=where((mxdyear ge useper(0)) and (mxdyear le useper(1)),mxdnyr)
mxdyear=mxdyear(kmxd)
fdcalibu=fdcalibu(*,*,kmxd)
fdmask=total(fdcalibu,3)*0.+1.
kmask=where(finite(fdmask),nmask)
print,'NYR=',mxdnyr,'  NBOX=',nmask
fdcalibu=reform(fdcalibu,g.nx*g.ny,mxdnyr)
;
; Optionally pre-filter (high-pass) each grid-box time series
;
if usefilt ne 0 then begin
  print,'Filtering data'
  for i = 0 , nmask-1 do begin 
    print,i,format='($,I4)'
    onets=reform(fdcalibu(kmask(i),*))
    filter_cru,thalf,/nan,tsin=onets,tshigh=onehi,tslow=onelo
    if usefilt eq 1 then fdcalibu(kmask(i),*)=onehi(*) $
                    else fdcalibu(kmask(i),*)=onelo(*)
  endfor
  print
endif
;
; Compute them
;
print,'Computing EOFs'
myeof1d_rot,fdcalibu,ev,ea,lam,lampct,lamcum,$
  evrot,earot,lamrot,lampctrot,lamcumrot,$
  nretain=nretain,correlation=usecorr
;
; Choose rotated or unrotated EOFs
;
if userot eq 1 then begin
  ev=evrot & ea=earot & lam=lamrot & lampct=lampctrot & lamcum=lamcumrot
  ; sort into descending order of variance explained
  isort=reverse(sort(lampct))
  ev=ev(*,isort)
  ea=ea(*,isort)
  lam=lam(isort)
  lampct=lampct(isort)
  for i = 0 , nretain-1 do lamcum(i)=total(lampct(0:i))
endif
;
; Now plot them
;
loadct,39
multi_plot,nrow=3,ncol=2,layout='large'
if !d.name eq 'X' then begin
  window,ysize=850
  !p.font=-1
endif else begin
  !p.font=0
  device,/bold,/helvetica,font_size=14
endelse
;
lampct=[!values.f_nan,lampct]
lamcum=[0,lamcum]
xt=mxdyear
plot,lampct,/xstyle,xtitle='Mode',ytitle='% variance explained',$
  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]
plot,[0,1],/nodata,xstyle=5,ystyle=5
;
if usedecline eq 0 then ttt='' else ttt=' (corrected)'
xyouts,0,0.9,'PCA of MXD gridded data'+ttt
ttt=string(useper,format='(2I5)')
xyouts,0,0.6,'Period:'+ttt
case usefilt of
  -1: ttt='Pre-filtered ( > '+string(thalf,format='(I3)')+' )'
  0: ttt='Unfiltered'
  1: ttt='Pre-filtered ( < '+string(thalf,format='(I3)')+' )'
endcase
xyouts,0,0.3,ttt
if usecorr eq 0 then tt1='Covariance' else tt1='Correlation'
if userot eq 0 then tt2='unrotated' else tt2='rotated'
xyouts,0,0.,tt1+' matrix,'+tt2
;
map=def_map(/npolar)  &  map.limit(0)=25.
coast=def_coast(/get_device)
labels=def_labels(/off)
;
levels=findgen(21)*0.03-0.3
levels(0)=-1.
levels(10)=0.
levels(20)=1.
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=reform(ea(*,i))
  filter_cru,th,tsin=yt,tslow=tslow,/nan
  plot,xt,yt,title='Mode'+string(i+1,format='(I2)'),$
    /xstyle,$
    xtitle='Year',ytitle='Amplitude'
  oplot,xt,tslow,thick=5
  oplot,!x.crange,[0,0],linestyle=1
  ;
  fd=reform(ev(*,i),g.nx,g.ny)
  inter_boxfd,fd,g.x,g.y,$
    map=map,coast=coast,labels=labels,$
    levels=levels,c_colors=c_colors
endfor
;
; Now save the EOFs for later use
;
fn='mxdmodes_abdlow_'+string(useper,format='(2I4)')
if userot eq 0 then fn=fn+'unr' else fn=fn+'rot'
if usecorr eq 0 then fn=fn+'cov' else fn=fn+'cor'
if usedecline eq 0 then fn=fn+'una' else fn=fn+'adj'
case usefilt of
  -1: fn=fn+'l'+string(thalf,format='(I2)')
  0: fn=fn+'raw'
  1: fn=fn+'h'+string(thalf,format='(I2)')
endcase
;
save,filename=fn+'.idlsave',$
  g,mxdyear,mxdnyr,usefilt,thalf,g,ea,ev,lampct,$
  usedecline,useper,usecorr,nretain,userot
;
end