; ; Performs PCA on correlations between MXD and monthly temperature and ; precipitation. For precipitation, should use the partial correlations ; found when temperature is held constant. ; Uses covariance matrix, since higher correlations should have more ; weight. ; ivar=0 ; 0=precip partial correlations, 1=temperature correlations ; 2=precip correlations, 3=cloud partial, 4=cloud isqu=0 ; 0=r, 1=r**2 (but converted back to r afterwards, and sign kept) ; nplot=-2 ; number of EOFs to plot, -ve values suppress plotting of eigenvalues ; trv=1 ; selects tree-ring-variable: 0=MXD 1=TRW case trv of 0: fnadd='mxd' 1: fnadd='trw' 2: fnadd='mxd-trw' 3: fnadd='trwmxd' 4: fnadd='trwmxdmult' endcase titadd=strupcase(fnadd) ; if ivar eq 0 then msgn=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1] $ else msgn=[1,-1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1] ; ; Restore all correlations ; if ivar eq 3 then fncld='_cldtem' else fncld='' restore,filename=fnadd+'_moncorr'+fncld+'.idlsave' ; allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon ; ; Perform PCA analysis ; nret=15 case ivar of 0: data1=reform(allp(*,0:15,0)) ; precip (PARTIAL CORREL)! 1: data1=reform(allr(*,0:15,1)) ; temperature 2: data1=reform(allr(*,0:15,0)) ; precip 3: data1=reform(allp(*,0:15,0)) ; cloud (PARTIAL CORREL)! 4: data1=reform(allr(*,0:15,0)) ; cloud endcase ; ; Compute overall mean ; totval=total(data1,/nan) totnum=total(finite(data1)) totval=totval/float(totnum) ; ; Optionally convert to r^2 but keeping sign ; if isqu eq 1 then begin kpdat=data1 data1=data1^2 neglist=where(kpdat lt 0.,nneg) if nneg gt 0 then data1(neglist)=-data1(neglist) endif ; myeof1d,data1,ev,ea,lam,lampct,lamcum,nretain=nret ;,/correlation ; ; 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 ; ; May need to convert from r^2 back to r, allowing negatives still ; if isqu eq 1 then begin kpdat=ea ea=sqrt(abs(ea)) neglist=where(kpdat lt 0.,nneg) if nneg gt 0 then ea(neglist)=-ea(neglist) endif ; ; Now prepare for plotting ; 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,/helvetica,/bold,font_size=12 endelse def_1color,20,color='red' def_1color,21,color='deepblue' def_1color,22,color='lgrey' 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 ; case ivar of 0: mtit='monthly precipitation , holding temperature constant )' 1: mtit='monthly temperature )' 2: mtit='monthly precipitation )' 3: mtit='monthly cloud cover , holding temperature constant )' 4: mtit='monthly cloud cover )' endcase ; if nplot ge 0 then begin ; plot,findgen(nret)+1.,lam,xtitle='Mode',$ ytitle='Eigenvalue',$ title='Covariance matrix PCA of!Cr ( '+titadd+' , '+mtit ; 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 ; !p.multi[0]=0 nret=abs(nplot) ; for i = 0 , nret-1 do begin ; pause inter_boxfd,/nodata,map=map,coast=coast,labels=labels,$ title=titadd+' EOF mode #'+string(i+1,format='(I2)')+' ('+$ string(round(lampct(i)),format='(I3)')+' % )' ; for j = 0 , nchron-1 do begin onev=ev(j,i)*msgn(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 ; if (ivar eq 1) or (ivar eq 4) then ytit='C' else ytit='Partial c' if ivar eq 1 then mtit='temperature' else mtit='precipitation' if ivar ge 3 then mtit='cloud cover' ; xxx=findgen(18) if i eq 0 then yadd=totval else yadd=0. cpl_barts,xxx,[!values.f_nan,reform(ea(*,i))*msgn(i),!values.f_nan]+yadd,$ /xstyle,ytitle=ytit+'orrelation anomaly',$ /ystyle,yrange=[-0.7,0.7],$ ymargin=[9,9],outline=-1,bar_color=22,$ xticks=15,xtickv=findgen(16)+1.,xtickname=corrname(0:15),$ title='Monthly '+mtit if i eq 0 then begin if !d.name eq 'PS' then device,font_size=6 plots,!x.crange,[yadd,yadd],thick=4 xyouts,!x.crange(1)+0.2,yadd-0.02,'Mean' if !d.name eq 'PS' then device,font_size=12 endif ; endfor ; end