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