; ; Calibrates and plots regional MXD against regional temperature ; Computes standard errors on the reconstructions too. ; DOES IT FOR AGE-BANDED DENSITY REGIONS AND NH ; ; Define periods for calibration ; calper=[1881,1960] ; calibration period verper=[1961,1994] ; verification period ; dosmooth=0 ; 0=plot (& do errors for) annual data, 1=smoothed data thalf=10 donh=0 ; 0=do NH 1=do ALL ; if donh eq 0 then nhtit='NH' else nhtit='ALL' ; ; Get the hemispheric temperature series ; restore,filename='datastore/compbest_mxd_fixed1950.idlsave' if donh eq 0 then nhtemp=reform(comptemp(*,3)) $ else nhtemp=reform(comptemp(*,2)) ; Get rid of pre-1871 temperatures knh=where(timey lt 1871) nhtemp(knh)=!values.f_nan ; ; Get all the age-banded MXD series (regions, hemi and back to 1400) ; restore,filename='/cru/u2/f055/treeharry/banding/bandregnh_fixed.idlsave' ; 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+1)*!values.f_nan newfm(*,0:nreg-1)=transpose(fullmxd(*,*)) newfm(*,nreg)=nhts(*) fullmxd=newfm ; ; Get regional temperature series ; restore,filename='datastore/regbest_mxd_fixed1950.idlsave' ; 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)) ) bestmxd=fullmxd(kcomp,*) ; ; Combine hemispheric and regional temperatures ; regname=[regname,nhtit] newrt=fltarr(nyr,nreg+1)*!values.f_nan newrt(*,0:nreg-1)=regtemp(*,*) newrt(*,nreg)=nhtemp(*) regtemp=newrt nreg=nreg+1 ; ; Prepare for plotting ; loadct,39 multi_plot,nrow=3,layout='large',/landscape if !d.name eq 'X' then begin window,xsize=850 !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] ; calregts=fltarr(nfull,nreg)*!values.f_nan ; ; Process each region separately ; for ireg = 0 , nreg-1 do begin print,regname(ireg) ; ; Extract series ; tstemp=reform(regtemp(*,ireg)) tsmxd=reform(bestmxd(*,ireg)) tsfull=reform(fullmxd(*,ireg)) ; ; Identify calibration and verification subsets ; calkl=where( finite(tstemp+tsmxd) and $ (timey ge calper(0)) and (timey le calper(1)) , ncal ) verkl=where( finite(tstemp+tsmxd) and $ (timey ge verper(0)) and (timey le verper(1)) , nver ) print,'Calibration period:'+string(calper,format='(2I5)')+$ ' length:'+string(ncal,format='(I5)') ; print,'Verification period:'+string(verper,format='(2I5)')+$ ; ' length:'+string(nver,format='(I5)') ; ; Calibrate the MXD ; mkcalibrate,tsmxd(calkl),tstemp(calkl),fitcoeff=fitcoeff,autocoeff=autocoeff print,'RAW' print,' r alpha beta SEalpha SEbeta sig rmse' print,fitcoeff,format='(F6.2,4F8.4,F7.2,F8.4)' print,' a_mxd a_tem a_resid ess' print,autocoeff,format='(3F6.2,F8.1)' ; Scale up the standard error on the regression coefficient if there is ; autocorrelation in the residuals aresid=autocoeff(2) serialfac=(1.+abs(aresid))/(1.-abs(aresid)) sebeta=fitcoeff(4)*serialfac print,'SEbeta scaled up=',sebeta ; ; Generate calibrated record and uncertainties etc. ; Use the extended MXD record (after checking it matches the short one!) ; xxx=tsfull(kcomp) yyy=tsmxd tse=total( (xxx-yyy)^2 , /nan ) ; print,tse if tse gt 0.03 then message,'Series do not match' calmxd=tsfull*fitcoeff(2)+fitcoeff(1) calregts(*,ireg)=calmxd(*) ; ; Now compute the one standard error for each reconstructed value ; (optionally follow procedure for smoothed data) ; ; To do so, 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 ; ; SE due to uncertainty in (a) slope and (b) intercept coefficients, ; and (c) due to unexplained variance ; 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] se_a=lowpred*sebeta varfac=serialfac*sqrt(total(allwt^2)) < 1. ; don't let errors get bigger print,'Unexplained variance scaled (due to smoothing) by',varfac se_c=fitcoeff(6)*varfac prets=callow endif else begin se_a=predanom*sebeta se_c=fitcoeff(6) prets=calmxd endelse se_b=fitcoeff(3) serr=sqrt( se_a^2 + se_b^2 + se_c^2 ) ; 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 !p.multi[0]=0 plot,timey,tstemp,thick=2,$ /xstyle,xrange=[1870,1965],xtitle='Year',$ ytitle='Temperature anomaly (!Uo!NC wrt 1961-90)',$ title=regname(ireg) oplot,yrmxd,calmxd,color=11,thick=2 filter_cru,/nan,10.,tsin=tstemp,tslow=t1 ;oplot,timey,t1,thick=6 filter_cru,/nan,10.,tsin=calmxd,tslow=t1 ;oplot,yrmxd,t1,thick=6,color=11 ; 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 ; endfor ; end