; ; Creates a regional timeseries and adjust it for changing sample size ; JUST DOES ALL, AND DOES TIME-DEPENDENT RBAR!!! ; plot,[0,1] multi_plot,nrow=4,layout='large' if !d.name eq 'X' then begin window,ysize=900 initx endif ; ncid=ncdf_open('tree_dens_nh.nc') ncdf_diminq,ncid,'time',dummy,ntime ncdf_diminq,ncid,'station',dummy,nstat ncdf_varget,ncid,'year',x ncdf_varget,ncid,'density',density ncdf_attget,ncid,'density','missing_value',valmiss ncdf_varget,ncid,'fraction',weight ncdf_close,ncid ; misslist=where(density eq valmiss,nmiss) density(misslist)=!values.f_nan ; ; Get regional tree lists and rbar ; restore,filename='rbar_dens_18811960.idlsave' restore,filename='reglists.idlsave' ; i=0 ; do ALL only ;for i = 0 , nreg-1 do begin ; dens=density(*,treelist(0:ntree(i)-1,i)) wt=weight(*,treelist(0:ntree(i)-1,i)) ; ; Make a time-dependent rbar first ; nt=n_elements(x) rbar_td=fltarr(nt)*!values.f_nan print,nt for it = 0 , nt-1 do begin mask = where( finite(dens(it,*)) , nsamp) print,it,nsamp if nsamp gt 0 then begin if nsamp eq 1 then begin rbar_td(it)=1. ; if n=1 then n'=1 regardless of rbar endif else begin ; get reduced grandr redr=reform(grandr(mask,*)) redr=redr(*,mask) ; Total it and average it, but only do i<>j. totr=0. totn=0. for j = 1 , nsamp-1 do begin for k = 0 , j-1 do begin totr=totr+redr(j,k) ; should have no missing values I hope totn=totn+1. if redr(j,k) eq 1. then message,'Ooops!' endfor endfor rbar_td(it)=totr/totn endelse endif endfor ; n=total(dens*0.+1.,2,/nan) ; compute timeseries of number of chronols ncore=total(dens*0.+wt,2,/nan) ; compute timeseries of summed core fraction totwdens=total(dens*wt,2,/nan) totwvals=total(float(finite(dens))*wt,2) wdens=totwdens/float(totwvals) ; this is the weighted anomaly ; ; Now need to adjust for varying sample size ; ;xadj=mkeffective(wdens,n,rbar=allrbar(i),neff=neff) xadj=mkeffective(wdens,n,rbar=rbar_td,neff=neff) wdens=xadj ; ; Now normalise w.r.t. 1881-1960 ; mknormal,wdens,x,refperiod=[1881,1960],refmean=refmean,refsd=refsd print,refmean,refsd ; ; Now plot them ; pause filter_cru,20,tsin=wdens,tslow=tslow,/nan cpl_barts,x,wdens,title='Weighted mean normalised max density anomaly'+$ ' for region: '+regname(i),xtitle='Year',/xstyle,$ zeroline=tslow,yrange=[-8,4] oplot,x,tslow,thick=2 moms=moment(wdens(where(finite(wdens))),sdev=sdev) oplot,!x.crange,[0.,0.] oplot,!x.crange,[-2.,-2.],linestyle=1 ; plot,x,n,ytitle='Sample size',/xstyle ; plot,x,rbar_td,/xstyle,ytitle='rbar' plot,x,neff,/xstyle,ytitle='Effective sample size' ; ; Now save the weighted mean timeseries, for further analysis ; densadj=wdens save,filename='densadj_'+regname(i)+'.idlsave',x,densadj,n,neff,ncore ; ;endfor ; end