; ; Calibrates the age-banded density data against observed temperature ; for a selected region (or quasi-hemispheric mean), using various ; options for focussing on low-frequency term only, and various ; options for quantifying uncertainty. ; doplot=2 ; 0=none, 1=time series, 2=scatter of parameters, ; 3=mean parameters over timscales ; ; Select which version of the age-banded records to use. fnband='_50-300' ; age bands used. Choices: '_50-300' ''=20-550 fnmint='_min2' ; minimum trees in each band. Choices: '_min2' ''=all used fnepsw='' ; weighting used for all. Choices: '_eps' ''=no weighting ; dosave=0 ; save results or not? 0=no, 1=IDLSAVE, 2=IDLSAVE & ASCII ; ; Define periods and time scales for calibration ; calper=[1871,1958] ; calibration period ; caltscale=[-10,-9,-8,-7,-6,-5,-4,-3,-2,1,2,3,4,5,6,7,8,9,10] ; smoothing to apply prior to calibration ;caltscale=[9] caltscale=[1] ncalts=n_elements(caltscale) ; errtscale=[1] nerrts=n_elements(errtscale) ; ; Select region(s) to calibrate (0-8=regions, 9=all, 10=nh) ; listreg=[9] ndoreg=n_elements(listreg) ; ; Select sequence of fixed grids to calibrate, ending in ; the "best coverage" grid. Also list the times between ; which to interpolate the calibration error statistics, ; for generating a smoothly varying error through time. ; ; All available fixed grids allfixedyr=[1400,1500,1600,1700,1800,1950,1988] nfixed=n_elements(allfixedyr) ; ; Which to use and in what order iord=[0,1,2,3,4,6,5] ; order to process (keep 1950 last!) iord=[5] nord=n_elements(iord) ; ; When to apply the values, and which value to use useyr= [1400,1500,1600,1700,1800,1900,1983,1994] usefix=[1400,1500,1600,1700,1800,1950,1950,1988] ; ; Prepare for plotting ; loadct,39 multi_plot,nrow=4,ncol=2,layout='large',/landscape if !d.name eq 'X' then begin wintim,xsize=850,ysize=750 !p.font=-1 endif else begin !p.font=0 device,/helvetica,/bold,font_size=9 endelse def_1color,10,color='lblue' def_1color,11,color='red' def_1color,12,color='vlblue' ; range0=[-2.5,-2.0,-3.5,-2.5,-1.5,-0.5,-2.5,-2.0,-3.0,-1.5] range1=[ 1.5, 2.0, 2.5, 1.5, 1.0, 1.0, 1.5, 1.5, 1.5, 1.0] ; ; Get the hemispheric temperature series ; fn='datastore/compbest_mxd_fixed1950.idlsave' print,'Observed NH and ALL temperature' print,fn print restore,filename=fn nhtemp=reform(comptemp(*,3)) alltemp=reform(comptemp(*,2)) ; Get rid of pre-1871 temperatures knh=where(timey lt 1871) nhtemp(knh)=!values.f_nan alltemp[knh]=!values.f_nan ; ; Get all the age-banded MXD series (regions, hemi and back to 1400) ; fn='/cru/u2/f055/treeharry/banding/bandregnh'+$ fnband+fnmint+fnepsw+'_fixed.idlsave' print,'Full tree time series from full coverage data' print,fn print restore,filename=fn ; Gets: timey,nyr,nhts,allts,alleps,nreg,regname yrmxd=timey fullmxd=allts ; Combine hemispheric MXD with regional MXD nfull=n_elements(yrmxd) newfm=fltarr(nfull,nreg+2)*!values.f_nan newfm(*,0:nreg-1)=transpose(fullmxd(*,*)) newfm(*,nreg)=nhts(*) ; predictor is same for ALL and NH, so repeat newfm(*,nreg+1)=nhts(*) fullmxd=newfm ; ; Get regional temperature series ; fn='datastore/regbest_mxd_fixed1950.idlsave' print,'Observed regional temperature' print,fn print restore,filename=fn ; Gets: nreg,regname,regtemp,rawtemp,timey,bestmxd etc. ; ; Identify overlap between instrumental and MXD series ; nyr=n_elements(timey) kcomp=where( (yrmxd ge timey(0)) and (yrmxd le timey(nyr-1)) ,ncomp ) ; ; Combine hemispheric and regional temperatures ; newrt=fltarr(nyr,nreg+2)*!values.f_nan newrt(*,0:nreg-1)=regtemp(*,*) newrt(*,nreg)=alltemp(*) newrt(*,nreg+1)=nhtemp(*) regtemp=newrt tempyear=timey tempnyr=nyr ; ; Need to read in all the series generated by various fixed grids, ; then cut down to just the period that overlaps with observations ; for ifixed = 0 , nfixed-1 do begin ; fyr=allfixedyr(ifixed) fnadd=string(fyr,format='(I4)') if fnadd eq '1950' then fnadd='' ; ; Get all the age-banded MXD series (regions, hemi and back to 1400) ; fn='/cru/u2/f055/treeharry/banding/bandregnh'+$ fnband+fnmint+fnepsw+'_fixed'+fnadd+'.idlsave' if ifixed eq 0 then print,'Tree time series from fixed grid coverage data' print,fn if ifixed eq nfixed-1 then print restore,filename=fn ; Gets: timey,nyr,nhts,allts,alleps,nreg,regname ; if ifixed eq 0 then begin newrm=fltarr(tempnyr,nreg+2,nfixed)*!values.f_nan endif newrm(0:ncomp-1,0:nreg-1,ifixed)=transpose(allts(*,kcomp)) newrm(0:ncomp-1,nreg,ifixed)=nhts(kcomp) newrm(0:ncomp-1,nreg+1,ifixed)=nhts(kcomp) ; endfor ; bestmxd=newrm regname=[regname,'ALL','NH'] nreg=nreg+2 nyr=tempnyr timey=tempyear ; ; Create arrays for storing calibrations and uncertainties ; calregts=fltarr(nfull,ndoreg)*!values.f_nan calregse=fltarr(nfull,ndoreg)*!values.f_nan ; ; Repeat for each region separately ; for jreg = 0 , ndoreg-1 do begin ireg=listreg[jreg] print,regname[ireg] ; meanfitcoeff=fltarr(nord,7,ncalts)*!values.f_nan meanautocoeff=fltarr(nord,4,ncalts)*!values.f_nan ; ; Repeat for each fixed grid calculation ; for jfixed = 0 , nord-1 do begin ifixed=iord(jfixed) fyr=allfixedyr(ifixed) print,'FIXED for',fyr ; ; Extract series ; tstemp=reform(regtemp(*,ireg)) tsmxd=reform(bestmxd(*,ireg,ifixed)) tsfull=reform(fullmxd(*,ireg)) ; for icalts = 0 , ncalts-1 do begin ; ; If requested, filter the series prior to calibration ; Two possible options are (1) running smoothing (e.g. gaussian weighting), ; and (2) discrete means (not running means). The latter has the disadvantage ; that some data is being thrown away (by taking e.g. non-overlapping decadal ; mean values), but the advantage that the filtering does not impose any ; extra autocorrelation on the data. The solution used here is to do (2), but ; repeated for all options of the starting point for computing discrete means, ; and then the regression parameters are averaged over all outcomes. ; if abs(caltscale[icalts]) le 1 then begin ; ; No smoothing required ; nser=1 tempser=reform(tstemp,tempnyr,nser) mxdcalser=reform(tsmxd,tempnyr,nser) mxdallser=reform(tsfull,nfull,nser) timeyser=timey ; endif ; if caltscale[icalts] lt -1 then begin ; ; Simple Gaussian-weighted filtering is done ; filter_cru,/nan,-caltscale[icalts],tsin=tstemp,tslow=tempser filter_cru,/nan,-caltscale[icalts],tsin=tsmxd,tslow=mxdcalser filter_cru,/nan,-caltscale[icalts],tsin=tsfull,tslow=mxdallser nser=1 tempser=reform(tempser,tempnyr,nser) mxdcalser=reform(mxdcalser,tempnyr,nser) mxdallser=reform(mxdallser,nfull,nser) timeyser=timey ; endif ; if caltscale[icalts] gt 1 then begin ; ; Non-overlapping running means, for all registration positions ; nser=caltscale[icalts] tempser=mkallmeans(tstemp,caltscale[icalts]) mxdcalser=mkallmeans(tsmxd,caltscale[icalts]) mxdallser=mkallmeans(tsfull,caltscale[icalts]) ; dummy=reform(tempser[*,0]) ndummy=n_elements(dummy) timeyser=timey[0]+findgen(ndummy)*caltscale[icalts]+0.5*(caltscale[icalts]-1) ; if (fyr eq '1950') and (doplot eq 1) then begin for i = 0 , caltscale[icalts]-1 do begin pause plot,timey,tstemp,/xstyle oplot,timeyser+i,tempser[*,i],psym=10 oplot,timey,tsmxd*0.4,color=240 oplot,timeyser+i,mxdcalser[*,i]*0.4,psym=10,color=240 endfor endif ; endif ; ; Now calibrate each of the generated series ; allfitcoeff=reform(fltarr(7,nser),7,nser)*!values.f_nan allautocoeff=reform(fltarr(4,nser),4,nser)*!values.f_nan help,allfitcoeff,allautocoeff for iser = 0 , nser-1 do begin ; ; Identify calibration and verification subsets ; x=reform(tempser[*,iser]) y=reform(mxdcalser[*,iser]) ty=timeyser+iser calkl=where(finite(x+y) and (ty ge calper(0)) and (ty le calper(1)),ncal) if jfixed eq 0 then begin print,'Calibration period:'+string(calper,format='(2I5)')+$ ' length:'+string(ncal,format='(I5)') endif ; ; Calibrate the MXD ; mkcalibrate,y(calkl),x(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff ; allfitcoeff[*,iser]=fitcoeff allautocoeff[*,iser]=autocoeff ; if iser+jfixed eq 0 then print,' Fyr SEalpha SEbeta rmse a_resid r alpha beta' fixcoeff=[fitcoeff([3,4,6]),autocoeff(2),fitcoeff([0,1,2])] print,fyr,fixcoeff,format='(I4,3F8.4,2F8.2,2F8.4)' ; ; Generate the calibrated record and the residuals ; xpred=fitcoeff[1]+fitcoeff[2]*y xresid=xpred-x if doplot eq 1 then begin pause plot,ty,x,/xstyle,psym=10 oplot,ty,xpred,color=240,psym=10 pause plot,ty,xresid,/xstyle,psym=10 endif ; ; Analyse the autoregressive and autocorrelation structure of residuals ; maxlag=10+20*(caltscale[icalts] lt 3) maxlag=maxlag < (ncal-1) xlag=findgen(maxlag+1) maxord=7 tjo_tsauto,xresid[calkl],maxlag=maxlag,maxorder=maxord,acf=acf,$ akaike=akaike,printakaike=(nser eq 1),wantorder=[1],arwant=ar1,$ bestorder=bestorder,arbest=arbest print,'BEST ORDER =',bestorder if iser eq 0 then allakaike=fltarr(nser,2,maxord+1)*!values.f_nan allakaike[iser,*,*]=akaike[*,*] ; ; Plot autocorrelation function ; if doplot eq 2 then begin pause noer=((iser gt 0) and (doplot eq 0)) if noer gt 0 then !p.multi[0]=!p.multi[0]+1 plot,xlag,acf,psym=def_sym(10),yrange=[-0.4,1],/ystyle,$ /xstyle,noclip=1,color=20+iser*20 oplot,!x.crange,[0,0],linestyle=1 ; ; Compute and plot the AR1 model ; acf_ar1=tjo_ar2acf(ar1[0,2],maxlag=maxlag) oplot,xlag,acf_ar1,color=20+iser*20 ; ; Compute and plot the best order AR model ; if bestorder gt 0 then begin acf_arbest=tjo_ar2acf(arbest[2:*],maxlag=maxlag) oplot,xlag,acf_arbest,color=20+iser*20,thick=3 endif ; endif ; endfor ; meanakaike=total(allakaike,1)/float(nser) print print,'Mean Akaike results' print print,'Autoregressive model order selection' print print,' Lag: m r(m) se2(m) AIC(m)' for ilag = 0 , maxord do begin print,ilag,acf[ilag],meanakaike[*,ilag],format='(I8,F8.3,F8.4,F8.3)' endfor print ; meanfitcoeff[jfixed,*,icalts]=total(allfitcoeff,2)/float(nser) meanautocoeff[jfixed,*,icalts]=total(allautocoeff,2)/float(nser) if doplot eq 2 then begin pause plot,allfitcoeff[0,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='r',/ynozero oplot,!x.crange,replicate(meanfitcoeff[jfixed,0,icalts],2) pause plot,allfitcoeff[1,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='intercept',/ynozero oplot,!x.crange,replicate(meanfitcoeff[jfixed,1,icalts],2) pause plot,allfitcoeff[2,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='slope',/ynozero oplot,!x.crange,replicate(meanfitcoeff[jfixed,2,icalts],2) pause plot,allfitcoeff[3,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='SE intercept',/ynozero oplot,!x.crange,replicate(meanfitcoeff[jfixed,3,icalts],2) pause plot,allfitcoeff[4,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='SE slope',/ynozero oplot,!x.crange,replicate(meanfitcoeff[jfixed,4,icalts],2) pause plot,allfitcoeff[6,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],title='rmse',/ynozero oplot,!x.crange,replicate(meanfitcoeff[jfixed,6,icalts],2) pause plot,allautocoeff[2,*],psym=def_sym(10),/xstyle,xrange=[-1,nser],$ title='residuals autocorrelation',/ynozero oplot,!x.crange,replicate(meanautocoeff[jfixed,2,icalts],2) endif ; endfor ; endfor ; if doplot eq 3 then begin !x.range=[min(caltscale)-1,max(caltscale)+1] !x.style=1 pause plot,caltscale,meanfitcoeff[0,0,*],psym=def_sym(10),title='r',/ystyle,yrange=[0.3,0.9] for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,0,*],psym=def_sym(10),color=50+30*jfixed pause plot,caltscale,meanfitcoeff[0,1,*],psym=def_sym(10),title='intercept',/ystyle,yrange=[-0.4,-0.14] for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,1,*],psym=def_sym(10),color=50+30*jfixed pause plot,caltscale,meanfitcoeff[0,2,*],psym=def_sym(10),title='slope',/ystyle,yrange=[0.1,0.5] for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,2,*],psym=def_sym(10),color=50+30*jfixed oplot,!x.crange,[0.42,0.42],linestyle=1 pause plot,caltscale,meanfitcoeff[0,3,*],psym=def_sym(10),title='SE intercept',/ystyle,yrange=[0.02,0.09] for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,3,*],psym=def_sym(10),color=50+30*jfixed pause plot,caltscale,meanfitcoeff[0,4,*],psym=def_sym(10),title='SE slope',/ystyle,yrange=[0.02,0.16] for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,4,*],psym=def_sym(10),color=50+30*jfixed pause ar1=reform(meanautocoeff[0,2,*]) plot,caltscale,meanfitcoeff[0,4,*]*sqrt((1.+abs(ar1))/(1.-abs(ar1))),$ psym=def_sym(10),title='SE slope (badly adjusted for AR1)',/ystyle,yrange=[0.04,0.3] for jfixed = 0 , nord-1 do begin ar1=reform(meanautocoeff[jfixed,2,*]) plots,caltscale,meanfitcoeff[jfixed,4,*]*sqrt((1.+abs(ar1))/(1.-abs(ar1))),$ psym=def_sym(10),color=50+30*jfixed endfor pause plot,caltscale,meanfitcoeff[0,6,*],psym=def_sym(10),title='rmse',/ystyle,yrange=[0.11,0.32] for jfixed = 0 , nord-1 do plots,caltscale,meanfitcoeff[jfixed,6,*],psym=def_sym(10),color=50+30*jfixed pause plot,caltscale,meanautocoeff[0,2,*],psym=def_sym(10),$ title='residuals autocorrelation',/ystyle,yrange=[-0.2,1] for jfixed = 0 , nord-1 do plots,caltscale,meanautocoeff[jfixed,2,*],psym=def_sym(10),color=50+30*jfixed for jfixed = 0 , nord-1 do xyouts,-10,-0.1+0.15*jfixed,string(allfixedyr[iord[jfixed]],format='(I4)'),$ color=50+30*jfixed oplot,!x.crange,[0.1,0.1],linestyle=1 !x.range=[0,0] !x.style=0 endif ; endfor ; if 1 eq 0 then begin ; if ncal gt 30 then begin ; ; Generate calibrated record and uncertainties etc. ; Use the extended MXD record (after checking it matches the short one!) ; And only for the best fixed grid ; if fyr eq 1950 then begin xxx=tsfull(kcomp) yyy=tsmxd tse=total( (xxx-yyy)^2 , /nan ) ; print,tse if (tse gt 0.03) or (finite(tse) eq 0) then message,'Series do not match' calmxd=tsfull*fitcoeff(2)+fitcoeff(1) calregts(*,ireg)=calmxd(*) ; ; To compute errors, need the predictor in terms of anomalies wrt to ; its mean over the calibration period predanom=tsfull mkanomaly,predanom,yrmxd,refperiod=calper,refmean=predmean print,'Predictor mean over cal. per. is',predmean ; ; Generate time dependent standard error ; filter_cru,/nan,thalf,tsin=tstemp,tslow=templow filter_cru,/nan,thalf,tsin=calmxd,tslow=callow if dosmooth ne 0 then begin filter_cru,/nan,thalf,tsin=predanom,tslow=lowpred,weight=filwt nwt=n_elements(filwt) allwt=[reverse(filwt(1:nwt-1)),filwt] prets=callow endif else begin prets=calmxd endelse fixcoeff(*,nfixed)=fixcoeff(*,nfixed-1) ; move 1988 values to end fixcoeff(*,nfixed-1)=fixcoeff(*,nfixed-2) ; repeat 1950 values for 1983 for jjj = nfixed-3 , 0 , -1 do begin ; fill in missing values if finite(fixcoeff(0,jjj)) eq 0 then begin fixcoeff(*,jjj)=fixcoeff(*,jjj+1) endif endfor if finite(fixcoeff(0,nfixed)) eq 0 then begin ; and end value fixcoeff(0,nfixed)=fixcoeff(0,nfixed-1) endif for iyr = 0 , nfull-1 do begin onefc3=interpol(reform(fixcoeff(0,*)),useyr,yrmxd(iyr)) onefc4=interpol(reform(fixcoeff(1,*)),useyr,yrmxd(iyr)) onefc6=interpol(reform(fixcoeff(2,*)),useyr,yrmxd(iyr)) oneac2=interpol(reform(fixcoeff(3,*)),useyr,yrmxd(iyr)) onefc3=onefc3(0) onefc4=onefc4(0) onefc6=onefc6(0) oneac2=oneac2(0) ; ; Now compute the one standard error for each reconstructed value ; (optionally follow procedure for smoothed data), using error ; parameters interpolated from various fixed grid series. ; ; SE due to uncertainty in (a) slope and (b) intercept coefficients, ; and (c) due to unexplained variance ; ; Scale up the standard error on the regression coefficient if there is ; autocorrelation in the residuals ; serialfac=(1.+abs(oneac2))/(1.-abs(oneac2)) sebeta=onefc4*serialfac ; if dosmooth ne 0 then begin se_a=lowpred(iyr)*sebeta varfac=serialfac*sqrt(total(allwt^2)) < 1. ; don't let errors get bigger se_c=onefc6*varfac endif else begin se_a=predanom(iyr)*sebeta se_c=onefc6 endelse se_b=onefc3 serr=sqrt( se_a^2 + se_b^2 + se_c^2 ) calregse(iyr,ireg)=serr ; endfor serr=calregse(*,ireg) ; if regname(ireg) eq nhtit then !p.multi=[0,1,2,0,0] pause ; plot,timey,tstemp,/nodata,$ /xstyle,xrange=[1400,1994],xtitle='Year',$ ytitle='Temperature anomaly (!Uo!NC wrt 1961-90)',$ /ystyle,yrange=[range0(ireg),range1(ireg)],$ title=regname(ireg) ; xfill=[yrmxd,reverse(yrmxd)] yfill=[prets-serr*2.,reverse(prets+serr*2.)] kl=where(finite(yfill),nkeep) yfill=yfill(kl) xfill=xfill(kl) polyfill,xfill,yfill,color=12 ; xfill=[yrmxd,reverse(yrmxd)] yfill=[prets-serr,reverse(prets+serr)] kl=where(finite(yfill),nkeep) yfill=yfill(kl) xfill=xfill(kl) polyfill,xfill,yfill,color=10 ; axis,yaxis=0,/ystyle,ytickformat='nolabels' axis,yaxis=1,/ystyle,ytickformat='nolabels' ; if dosmooth eq 0 then begin oplot,timey,tstemp oplot,yrmxd,calmxd,color=11 endif ; oplot,!x.crange,[0.,0.],linestyle=1 oplot,timey,templow,thick=3 oplot,yrmxd,callow,thick=3,color=11 ; if (regname(ireg) eq nhtit) and (dosmooth eq 0) then begin plot,timey,tstemp,thick=3,$ /xstyle,xrange=[1850,2000],xtitle='Year',$ ytitle='Temperature anomaly (!Uo!NC wrt 1961-90)',$ title=regname(ireg) oplot,yrmxd,calmxd,color=11,thick=3 oplot,yrmxd,calmxd+serr,color=10 oplot,yrmxd,calmxd-serr,color=10 oplot,yrmxd,calmxd+serr*2.,color=10 oplot,yrmxd,calmxd-serr*2.,color=10 oplot,!x.crange,[0.,0.],linestyle=1 endif ; endif endif endif ; if dosave gt 0 then begin ; ; Save standard errors (note they apply only to the time scale chosen!) ; timey=yrmxd save,filename='bandtemp_mxd_errors.idlsave',$ calregse,thalf,dosmooth,nreg,regname,timey ; endif ; end