; ; 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