; ; Reads in site-by-site MXD and temperature series in 5 yr blocks, all ; correctly normalised etc. Rotated PCA is performed to obtain the 'decline' ; signal! ; ;---------------------------------------------------- ; restore,filename='briffa_sep98_decline1.idlsave' ; allmxd,alltemp,mxd5yr,temp5yr,nchron,nyr,nblock,statlat,statlon ; blockyr,timey ; nret=4 diff5yr=mxd5yr-temp5yr ; myeof1d,diff5yr,ev,ea,lam,lampct,lamcum,nretain=nret,correlation=1 myeof1d_rot,diff5yr,ev,ea,lam,lampct,lamcum,$ evrot,earot,lamrot,lampctrot,lamcumrot,$ nretain=nret,correlation=1 ev=evrot ea=earot lam=lamrot lampct=lampctrot lamcum=lamcumrot ; ; Scale EOFs and PCs so that highest value of each EOF is +1 ; for i = 0 , nret-1 do begin dummy=max(abs(ev(*,i)),imax) maxval=ev(imax,i) ev(*,i)=ev(*,i)/maxval ea(*,i)=ea(*,i)*maxval endfor ; ; Sort modes into order of descending variance explained ; i=reverse(sort(lampct)) lam=lam(i) lampct=lampct(i) lamcum(0)=lampct(0) for j = 1 , nret-1 do lamcum(j)=lamcum(j-1)+lampct(j) ev=ev(*,i) ea=ea(*,i) ; ; Now prepare for plotting ; loadct,39 multi_plot,nrow=2,ncol=2,layout='large' if !d.name eq 'X' then begin window,ysize=850 !p.font=-1 endif else begin !p.font=0 device,/helvetica,/bold,font_size=12 endelse def_1color,20,color='red' def_1color,21,color='deepblue' def_1color,22,color='lgrey' def_1color,95,color='dgreen' cpl_usersym,/circle,/fill ; map=def_map(/npolar) & map.limit(0)=25. coast=def_coast(/get_device) & coast.double=1 labels=def_labels(/off) ; ; Plot variance stats first ; if 2 eq 3 then begin plot,findgen(nret)+1.,lam,xtitle='Mode',$ ytitle='Eigenvalue',$ title='Correlation matrix PCA of (MXD minus TEMP)' ; plot,findgen(nret)+1.,alog(lam),xtitle='Mode',$ ytitle='Log of eigenvalue' ; plot,findgen(nret+1),[!values.f_nan,lampct],xtitle='Mode',$ ytitle='Explained variance (%)',psym=10 ; plot,findgen(nret+1),[0.,lamcum],xtitle='Mode',$ ytitle='Cumulative explained variance (%)',psym=10,yrange=[0,100] endif ; ; Now plot some of the retained EOFs and PCs ; nret=4 ; for i = 1 , 2 do begin ; pause inter_boxfd,/nodata,map=map,coast=coast,labels=labels,$ title='EOF mode #'+string(i+1,format='(I2)')+' ('+$ string(round(lampct(i)),format='(I3)')+' % )' ; for j = 0 , nchron-1 do begin onev=ev(j,i) if onev lt 0 then scol=21 else scol=20 ssiz=abs(onev)*0.8 ; scaling factor to get good sizes ssiz=ssiz+0.1 ; minimum size! map_plots,statlon(j),statlat(j),psym=8,color=scol,symsize=ssiz endfor ; plot,blockyr,ea(*,i),/xstyle,thick=2,$ xtitle='Year',ymargin=[10,10],xmargin=[4,4] oplot,blockyr,smooth(ea(*,i),5,/nan,/edge_truncate),thick=8,color=95 ; endfor ; end