; ; 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 isqu=0 ; 0=r, 1=r**2 (but converted back to r afterwards, and sign kept) ; 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' 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] pp='('+[ ['a','c','e','g','i','k'] , $ ['b','d','f','h','j','l'] ]+') ' ; ; Restore all correlations ; restore,filename=fnadd+'_moncorr.idlsave' ; allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon ; ; Perform PCA analysis ; nret=16 if ivar eq 1 then data1=reform(allr(*,0:15,1)) $ ; temperature else data1=reform(allp(*,0:15,0)) ; precip (PARTIAL CORREL)! ; ; 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 ; ; 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=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' 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 ivar eq 1 then mtit='monthly temperature )' $ else mtit='monthly precipitation , holding temperature constant )' ; 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] ; ; Now plot some of the retained EOFs and PCs ; nret=6 nret=2 nret=1 ; for i = 0 , nret-1 do begin ; pause inter_boxfd,/nodata,map=map,coast=coast,labels=labels,$ title=pp(i,0)+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 if onev lt 0 then ssym=15 else ssym=10 ssiz=abs(onev)*0.8 ; scaling factor to get good sizes ssiz=ssiz+0.1 ; minimum size! map_plots,statlon(j),statlat(j),psym=def_sym(ssym),symsize=ssiz endfor ; if ivar eq 1 then ytit='C' else ytit='Partial c' if ivar eq 1 then mtit='temperature' else mtit='precipitation' ; 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),$ xtitle='Monthly '+mtit,title=pp(i,1) 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