; ; Reads in the MXD data set, grids it onto a 5 by 5 grid, adjusts each grid ; box for sample size variations (requires rbar for each box, which is also ; computed here), and then normalises over the period 1881-1960 (although ; the period can be altered). When gridding, multiple chronologies are ; weighted by their "fraction of cores" values before averaging. ; Only those with complete data over 1800-1964 are kept. Those with complete ; data over 1800-1976 are used to infill missing values between 1965 and 1976, ; and those complete over 1697-1964 are used to infill missing values between ; 1697 and 1799. Any that are not complete over 1697-1976 after infilling ; are removed, and only the 1697-1976 period is kept and saved. ; *** NOW REMOVES MXD SERIES WITH POOR CORRELATIONS WITH LOCAL SUMMER TEMP *** ; *** NOW IT ALSO DOES LOTS OF INFILLING USING VARIOUS PERIODS *** ; ; *** NOW ALSO WORKS FOR TRW! *** ; ;---------------------------------------------------------------------- ; beglist=[1400,1453,1540,1583,1612,1660,1697,1743,1777] nbeg=n_elements(beglist) endlist=[1994,1988,1980] nend=n_elements(endlist) ; trv=1 ; selects tree-ring-variable: 0=MXD 1=TRW cutr=0.22 ; a 'best' series is based on sites with local r >= cutr (0.22) case trv of 0: begin fnadd='mxd' iseas=18 ; use local r with Apr-Sep temp to select chronologies end 1: begin fnadd='trw' iseas=20 ; use local r with Jun-Aug temp to select chronologies end endcase titadd=strupcase(fnadd) ; ; First read in the MXD data ; print,'Reading '+titadd+' data' restore,filename='../all'+fnadd+'.idlsave' ; nchron,idno,idname,location,country,tree,yrstart,yrend,statlat,statlon,$ ; mxd,fraction,timey,nyr ; ; Next read in the correlations between MXD and local climate ; print,'Reading '+titadd+' local correlations' restore,filename='../'+fnadd+'_moncorr.idlsave' ; allr,ncorr,nvar,nchron,varname,corrname,statlat,statlon,allp,moir,moir2 ; ; Remove any that are almost entirely missing, or that ; aren't well correlated with their local summer temperature ; nkeep=0 for i = 0 , nchron-1 do begin if allr(i,iseas,1) ge cutr then begin dummy=where(finite(mxd(*,i)),ngot) if ngot ge 20 then begin if nkeep eq 0 then kl=i else kl=[kl,i] nkeep=nkeep+1 endif endif endfor print,'Had',nchron,' chronologies, but',nchron-nkeep,' are empty or poor' nstat=nkeep density=mxd(*,kl) weight=fraction(*,kl) statlon=statlon(kl) statlat=statlat(kl) statcty=country(kl) ; ;---------------------------------------------------------------------- ; ; Get grid info (IPCC instrumental dataset grid), NH only ; g=def_grid(/ipcc) ky=where(g.y ge 0.,newny) g = { nx: g.nx, ny: newny, x: g.x, y: g.y(ky) } ; ;---------------------------------------------------------------------- ; ; Create a box by box list of which chronologies fall in each box ; print,'Locating the chronologies into grid boxes' chronlists=ingrid(g.x,g.y,statlon,statlat,nstat=chrondens) maxlist=max(chrondens) keepbox=where(chrondens gt 0.,maxbox) keepi=keepbox mod g.nx keepj=fix(keepbox/g.nx) chronlists=reform(chronlists,maxlist,g.nx*g.ny) chrondens=reform(chrondens,g.nx*g.ny) ; ;---------------------------------------------------------------------- ; ; Define array to contain box by box density timeseries ; boxdens=fltarr(maxbox,nyr) cty=strarr(maxbox) ; ;---------------------------------------------------------------------- ; ; For each box, extract chronologies, compute rbar if more than one, ; compute weighted mean chronology, adjust variance and normalise ; print,'Computing variance adjusted, normalised, box-mean timeseries' ;print,replicate('-',maxbox),format='(80A1)' refper=[1881,1960] for ibox = 0 , maxbox-1 do begin ;print,'.',format='($,A1)' nlist=chrondens(keepbox(ibox)) onelist=chronlists(0:nlist-1,keepbox(ibox)) mxd=density(*,onelist) wt=weight(*,onelist) cty(ibox)=statcty(onelist(0)) if nlist gt 1 then begin ; compute weighted mean of chronologies in the box totwtmxd=total(mxd*wt,2,/nan) totwt=total(float(finite(mxd))*wt,2) wtmxd=totwtmxd/float(totwt) ; compute rbar between chronologies inx=transpose(mxd) onerbar=mkrbar(inx) print,nlist,onerbar ; adjust variance for sample size nsample=total(mxd*0.+1.,2,/nan) adjmxd=mkeffective(wtmxd,nsample,rbar=onerbar) endif else begin ; if only one chronology in box, then no need to do anything to it print,nlist adjmxd=reform(mxd) endelse ; now normalise mknormal,adjmxd,timey,refperiod=refper boxdens(ibox,*)=adjmxd endfor print ; ; Now store the full gridded time series into a northern hemisphere array ; mxdfd=fltarr(g.nx,g.ny,nyr)*!values.f_nan for i = 0 , maxbox-1 do mxdfd(keepi(i),keepj(i),*)=boxdens(i,*) mxdyear=timey mxdnyr=nyr mxdfd2=mxdfd ; mxdfd2 will also contain the infilled data ; ;---------------------------------------------------------------------- ; ; Remove those that are too short (assume that they are continuous, so just ; look for the first non-missing value) ; maxstart=1810 print,'Removing boxes that dont precede',maxstart yrstart=fltarr(maxbox) yrend=fltarr(maxbox) for ibox = 0 , maxbox-1 do begin adjmxd=boxdens(ibox,*) keeplist=where(finite(adjmxd),nkeep) yrstart(ibox)=timey(keeplist(0)) yrend(ibox)=timey(keeplist(nkeep-1)) endfor keeplist=where(yrstart le maxstart,nkeep) print,'Reduce boxes from',maxbox,' to',nkeep maxbox=nkeep keepbox=keepbox(keeplist) keepi=keepi(keeplist) keepj=keepj(keeplist) boxdens=boxdens(keeplist,*) cty=cty(keeplist) yrstart=yrstart(keeplist) yrend=yrend(keeplist) ; ; Make a copy containing un-infilled data, and one to contain the final ; multi-infilled data ; infilldens=boxdens nofilldens=boxdens ; ;---------------------------------------------------------------------- ; ; Locate those that are complete for 1800 to 1985 and use them to infill ; those that stop between 1965 and 1985 ; *** NOW TRY IT FOR 1697 to 1976 (280 years) *** ; for iend = 0 , nend-1 do begin lastyr=endlist(iend) print,'Infilling to',lastyr baselist=where(yrend ge lastyr,nbase) print,'Out of',maxbox,' boxes,',nbase,' reach',lastyr if nbase eq 0 then message,'Cannot infill' filllist=where(yrend lt lastyr,nfill) print,'Fill in',nfill ncant1=0 if nfill gt 0 then begin boxdens=nofilldens infill,boxdens,g.x(keepi),g.y(keepj),timey,boxname=cty,$ baselist=baselist,filllist=filllist,$ overlap=[1810,1970],fillperiod=[1971,lastyr],$ search=[800.,1600.,2400.,3000.],$ cantlist=cantlist1 ncant1=n_elements(cantlist1) print,'Cannot fill',ncant1 newlist=where(finite(boxdens),nnew) if nnew gt 0 then infilldens(newlist)=boxdens(newlist) endif endfor ; ;---------------------------------------------------------------------- ; ; Locate those that are complete for 1700 to 1964 and use them to infill ; those that start between 1700 and 1800 ; for ibeg = 0 , nbeg-1 do begin begyr=beglist(ibeg) print,'Infilling to',begyr baselist=where(yrstart le begyr,nbase) print,'Out of',maxbox,' boxes,',nbase,' reach',begyr if nbase eq 0 then message,'Cannot infill' filllist=where(yrstart gt begyr,nfill) print,'Fill in',nfill if nfill gt 0 then begin boxdens=nofilldens infill,boxdens,g.x(keepi),g.y(keepj),timey,boxname=cty,$ baselist=baselist,filllist=filllist,$ overlap=[1810,1970],fillperiod=[begyr,1809],$ search=[800.,1600.,2400.,3000.],$ cantlist=cantlist2 ncant2=n_elements(cantlist2) print,'Cannot fill',ncant2 newlist=where(finite(boxdens),nnew) if nnew gt 0 then infilldens(newlist)=boxdens(newlist) endif endfor ; ; Store the infilled ones in the full array ; boxdens=infilldens for i = 0 , maxbox-1 do mxdfd2(keepi(i),keepj(i),*)=boxdens(i,*) ; ;---------------------------------------------------------------------- ; ; Get rid of any that couldn't be infilled, and keep only 1700 to 1985 only. ; Any that contain any missing data during that period are removed. ; if ncant1 gt 0 then boxdens(cantlist1,*)=!values.f_nan if ncant2 gt 0 then boxdens(cantlist2,*)=!values.f_nan ; keeplist=where((timey ge begyr) and (timey le lastyr),nkeep) timey=timey(keeplist) boxdens=boxdens(*,keeplist) ntime=nkeep ; totdens=total(boxdens,2) keeplist=where(finite(totdens),nkeep) print,'Ended up with',nkeep,' boxes' keepi=keepi(keeplist) keepj=keepj(keeplist) boxdens=boxdens(keeplist,*) oldn=maxbox maxbox=nkeep ; ; Also work out where the un-infilled (complete) boxes are in the new ; reduced dataset ; ilist=fltarr(oldn) ; define 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ... ilist(baselist)=1. ; mark 0 0 1 1 0 1 0 0 1 1 1 0 0 1 1 0 1 0 0 ... ilist=ilist(keeplist) ; remove 0 0 1 1 0 0 1 1 1 0 1 1 0 0 ... baselist=where(ilist eq 1,nbase) ; locate remaining 1's ; ;---------------------------------------------------------------------- ; save,filename=fnadd+'_infill.idlsave',ntime,maxbox,boxdens,g,keepi,keepj,timey,$ baselist,nbase,mxdfd,mxdyear,mxdnyr,mxdfd2,beglist,endlist ; end