;
; Computes EOFs of observed Apr-Sep land air temperature, to test and compare
; different methods.
; (1) Compare co-located versus full coverage
; (2) Compare high-pass-filtered versus no pre-filtering
;
usecoloc=1      ; 0=full data, 1=co-located with MXD sites only
usefilt=0       ; 0=unfiltered,  1=pre-filtered
thalf=40.       ; cutoff for high-pass filter
nretain=12
;
; Restore Apr-Sep land air temperature gridded dataset
;
print,'Reading in temperature data'
restore,filename='obs_lat_as.idlsave'
;  timey,fdseas,xlon,ylat,nyr,nx,ny,fdltm,fdsd,missfrac
;
; Restore regional breakdown of grid boxes
;
restore,filename='../reg_mxdboxes.idlsave'
;  nreg,boxlists,boxnumbs,regname,g,boxregs,boxprec,boxtemp
;
; Compute required mask to apply to data
;
print,'Masking out unused data'
if usecoloc eq 1 then fdmask=total(finite(boxlists),3) $
                 else fdmask=total(finite(boxregs),3)
kl=where(fdmask gt 0,nbox)
print,'Retaining maximum of boxes = ',nbox
fdmask=fltarr(nx,ny)*!values.f_nan
fdmask(kl)=1.
for iyr = 0 , nyr-1 do fdseas(*,*,iyr)=fdseas(*,*,iyr)*fdmask(*,*)
;
; Optionally pre-filter (high-pass) each grid-box time series
;
if usefilt eq 1 then begin
  print,'Filtering data: rows=',nx
  for i = 0 , nx-1 do begin 
    print,i,format='($,I4)'
    for j = 0 , ny-1 do begin
      onets=reform(fdseas(i,j,*))
      dummy=where(finite(onets),ngot)
      if ngot gt 5. then begin
        filter_cru,thalf,/nan,tsin=onets,tshigh=onehi
        fdseas(i,j,*)=onehi(*)
      endif
    endfor
  endfor
  print
endif
;
; Select the years to analyse
;
ksub=where(missfrac lt 0.29,nsub)
if nsub ne 66 then message,'Ooops!'
;
fdsub=fdseas(*,*,ksub)
yrsub=timey(ksub)
;
; Compute them
;
print,'Computing EOFs'
myeof2d_rot,fdsub,ev,ea,lam,lampct,lamcum,$
  nretain=nretain,correlation=1
;
; Try to infill the amplitude timeseries for those years not used in the EOF
; analysis (indeed, apply it to all the data to see if it gives the correct
; values!).
;
; First remove long-term mean (not 61-90 mean!)
print,'Infilling time series gaps'
fdanom=reform(fdseas,nx*ny,nyr)
mkanomaly,fdanom
fdanom=reform(fdanom,nx,ny,nyr)
allea=fltarr(nyr,nretain)
for iretain = 0 , nretain-1 do begin
  for iyr = 0 , nyr-1 do begin
    allea(iyr,iretain)=total(fdanom(*,*,iyr)*ev(*,*,iretain),/nan)
  endfor
endfor
fillea=allea
fillea(ksub,*)=ea(*,*)
;
; Now plot them
;
loadct,39
multi_plot,nrow=3,ncol=2,layout='large'
if !d.name eq 'X' then window,ysize=850
;
lampct=[!values.f_nan,lampct]
lamcum=[0,lamcum]
xt=timey
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
xyouts,0,0.9,'PCA of observed Apr-Sep temperature'
if usecoloc eq 1 then ttt='Co-located with MXD' else ttt='Full regions'
xyouts,0,0.6,ttt
if usefilt eq 1 then ttt='Pre-filtered ( < '+string(thalf,format='(I3)')+' )' $
                else ttt='Unfiltered'
xyouts,0,0.3,ttt
xyouts,0,0.,'Correlation matrix, unrotated'
;
map=def_map(/npolar)  &  map.limit(0)=15.
coast=def_coast(/get_device)
labels=def_labels(/off)
;
levels=[-0.5,findgen(19)*0.02-0.2,0.5]
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=fltarr(nyr)*!values.f_nan
  yt(ksub)=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'
  ;
  ml=where(finite(yt) eq 0)
  yt(ml)=allea(ml,i)
  oplot,xt,yt,linestyle=1
  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
;
if usecoloc eq 0 then begin
  save,filename='olat_modes.idlsave',$
    timey,nretain,usefilt,usecoloc,thalf,xlon,ylat,nyr,nx,ny,fillea,ev
endif
;
end