; ; Need to provide a list of years in vector yrlist ; yrlist=[1453,1601,1641,1695,1816] md=[2300,2300,2300,2300,2300] ; ; Get display ready ; loadct,39 def_1color,r,g,b,90,color='black' def_1color,r,g,b,91,color='purple' def_1color,r,g,b,92,color='blue' def_1color,r,g,b,93,color='green' def_1color,r,g,b,94,color='red' ; def_1color,r,g,b,40,color='purple' def_1color,r,g,b,42,color='dblue' def_1color,r,g,b,43,color='blue' def_1color,r,g,b,46,color='lgreen' def_1color,r,g,b,47,color='lyellow' def_1color,r,g,b,50,color='red' def_smearcolor,r,g,b,fromto=[40,42] def_smearcolor,r,g,b,fromto=[43,46] def_smearcolor,r,g,b,fromto=[47,50] multi_plot,nrow=3,ncol=2,layout='large' if !d.name eq 'X' then begin window, ysize=950 !p.font=-1 endif else begin !p.font=0 device,/helvetica,/bold,font_size=9,$ xsize=17.5,xoffset=1.5 endelse !x.omargin=[0,15] ; ; N. Hemi. map down to 30N, no labels ; map=def_map(/npolar) map.limit(0)=30. labels=def_labels(/off) ; ; 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=0.05 ; ; First plot a map of the regions and the chronology locations (with different ; colours according to chronology length ; inter_boxfd,/nodata,map=map,labels=labels ; !p.ticklen=0.03 lon_polar,map=map,[-180,-150.,-120.,-90.,-60.,-30.,30.,60.,120.,150.],/dotick ; x1=[180.,60.,-40.,-105.,-105.,-107.] y1=[30., 30., 30., 30., 55., 55.] x2=[180.,60.,-40.,-105.,-107.,-107.] y2=[90., 90., 90., 55., 55., 90.] n=n_elements(x1) for i = 0 , n-1 do begin map_plots,[x1(i),x2(i)],[y1(i),y2(i)],thick=6 endfor ; xs=findgen(99)-40. ys=replicate(54.8,99) map_plots,xs,ys,thick=6 ; ; Plots maps of chronology location ; ncid=ncdf_open('tree_dens_nh.nc') ncdf_varget,ncid,'latitude',ylat ncdf_varget,ncid,'longitude',xlon ncdf_varget,ncid,'year',timey ncdf_varget,ncid,'density',density ncdf_attget,ncid,'density','missing_value',valmiss ncdf_close,ncid misslist=where(density eq valmiss,nmiss) density(misslist)=!values.f_nan ; restore,filename='reglists.idlsave' iall=where(regname eq 'ALL') i=iall(0) mxdall=density(*,treelist(0:ntree(i)-1,i)) x=xlon(treelist(0:ntree(i)-1,i)) y=ylat(treelist(0:ntree(i)-1,i)) ; allcut=reverse([1500,1600,1700,1800,1910]) ncut=n_elements(allcut) ; cpl_usersym,/circle,/fill for icut = 0 , ncut-1 do begin ct=allcut(icut) yw=where(timey eq ct) yw=yw(0) kl=where(finite(mxdall(yw,*)),nkeep) print,icut,ct,nkeep,90+icut if ct eq 1500 then print,x(kl) plots,x(kl),y(kl),psym=8,color=90+icut,symsize=!p.charsize*0.9 endfor ; ; Now plot the yearly maps (read from file containing mxd normalised ; w.r.t. 1881-1960). ; nyr=n_elements(yrlist) for iyy = 0 , nyr-1 do begin iyr=yrlist(iyy) print,'Year:',iyr densfd=rd1yr_norm(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,'Number of chronology values:',nkeep ; ; First of all, store each station value into its 1. by 1. 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=2. & dy=2. 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) print,'Min/Max values:',min(gridfd,/nan,max=mmm) & print,mmm ; levels=[-12,-6,-4,-3,-2,-1,-0.5,0,0.5,1,2,3] nlev=n_elements(levels) zeroloc=where(levels eq 0,nzero) c_labels=intarr(nlev)+0 & if nzero gt 0 then c_labels(zeroloc)=0 c_charsize=0.6 c_colors=findgen(11)+40. ; inter_const,gridfd,gx,gy,map=map,$ maxdist=md(iyy),$ gs=[2.,2.],$ sm=sm,labels=labels,fdsm=fdsm,$ miss_grey='white',$ /cell_fill,levels=levels,c_colors=c_colors contour,fdsm.fd,fdsm.x,fdsm.y,/overplot,levels=levels,c_labels=c_labels,$ c_charsize=c_charsize,/follow ; xyouts,-150.,27.,string(iyr,format='(I4)'),align=1. ; endfor ; ; Now plot a colour scale ; device,font_size=12 !p.multi(0)=1 !p.region=[0.9,0.35,1.,0.65] !y.margin=map.ymargin scale_vert,levels=levels,$ c_colors=c_colors,/noextremes,$ title='Normalised!Cdensity!Canomaly' device,font_size=9 scale_vert,/off !p.region=[0,0,0,0] ; end