; ; Need to provide a list of years in vector yrlist ; nyr=n_elements(yrlist) if nyr eq 0 then message,'Need list of years in yrlist' ; ; Get display ready ; loadct,39 if nyr eq 2 then multi_plot,nrow=2,ncol=1 $ else multi_plot,nrow=3,ncol=2,layout='large' if !d.name eq 'X' then begin window,ysize=950 fac=1. endif else begin fac=2. endelse ; ; N. Hemi. map down to 30N, no labels ; map=def_map(/npolar) map.limit(0)=30. labels=def_labels(/off) labels.gridon=1 & labels.dlon=90. coast=def_coast(/get_device) coast.double=0 ; ; Small amount of smoothing, with filling in around the edges to extend to ; region that is plotted ; sm=def_sm() sm.method=2 sm.thresh=1.0 ; ;print,'Enter year required:' ;read,iyr ;yrlist=[1601,1641,1669,1699,1816,1817,1818,1878,1884,1912,1964] ;yrlist=[1818,1878,1884,1912,1964] for iyy = 0 , nyr-1 do begin iyr=yrlist(iyy) print,iyr densfd=rd1yr(iyr,x=x,y=y,frac=frac) ; ; Remove missing data and data south of 30N ; keeplist=where((densfd ne -9.99) and (y ge 30.),nkeep) if nkeep gt 0 then begin frac=frac(keeplist) densfd=densfd(keeplist) x=x(keeplist) y=y(keeplist) endif print,nkeep ; ; First of all, store each station value into its 0.5 by 0.5 grid box. ; When more than one falls in a box, average them - this is where the ; weighting by the fraction of cores available comes into it! It also ; prevents duplicate points going forward, which can upset spherical ; triangulation. ; dx=1. & dy=1. gnx=360./dx gny=(90.-30.)/dy gx=findgen(gnx)*dx-180. gy=90.-findgen(gny)*dy gridfd=gridit(gnx,gny,gx,gy,x,y,densfd,frac,nstat=chron) ; ; Convert boxes back to a list of stations ; gx2d=gx # (intarr(gny)+1.) gy2d=transpose(gy # (intarr(gnx)+1.)) keeplist=where(finite(gridfd)) gridfd=gridfd(keeplist) gx=gx2d(keeplist) gy=gy2d(keeplist) ; pause ; levels=findgen(21)-10. nlev=n_elements(levels) zeroloc=where(levels eq 0,nzero) c_labels=intarr(nlev)+1 & if nzero gt 0 then c_labels(zeroloc)=0 c_charsize=0.6 ;c_thick=(levels lt 0.)*fac+2. & if nzero gt 0 then c_thick(zeroloc)=1. c_thick=1.5*fac & if nzero gt 0 then c_thick(zeroloc)=1. print,levels print,c_thick ; inter_const,gridfd,gx,gy,map=map,$ maxdist=3000.,$ gs=[2.,2.],$ sm=sm,labels=labels,$ shade=1,sh_thresh=0.,sh_grey=[235,180],miss_grey='white',$ levels=levels,c_labels=c_labels,c_charsize=c_charsize,c_thick=c_thick,$ /follow,$ title=string(iyr,format='(I4)') ; xtitle='Tree ring density anomaly' ; endfor ; end