; ; Reads in density data, puts a low-frequency filter through each retained ; grid-box, and plots anomaly (normalised w.r.t. 1901-40) map ; ; Get display ready ; loadct,39 multi_plot,nrow=1,layout='large' if !d.name eq 'X' then begin window,ysize=700,xsize=700 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 ;sm.thresh=0.1 ; ; Read in data ; 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',frac ncdf_varget,ncid,'longitude',statlon ncdf_varget,ncid,'latitude',statlat ncdf_close,ncid ; misslist=where(density eq valmiss,nmiss) density(misslist)=!values.f_nan ; thalf=10. ; fns='dens1978.idlsave' dummy=findfile(fns,count=nfile) if nfile gt 0 then begin restore,filename=fns endif else begin ; ; For each available box, normalise w.r.t. 1901-40, and filter with thalf ; statvals=fltarr(nstat) iwant=where(x eq 1978) fracvals=frac(iwant,*) print,'nstat=',nstat for i = 0 , nstat-1 do begin ts1=reform(density(*,i)) dummy=where(finite(ts1) and (x ge 1901.) and (x le 1940.),nkkkk) print,i,nkkkk if nkkkk gt 10 then begin mknormal,ts1,x,refperiod=[1901,1940] filter_cru,thalf,tsin=ts1,tslow=ts3,/nan statvals(i)=ts3(iwant) endif else begin statvals(i)=!values.f_nan endelse endfor save,filename=fns,statvals,fracvals endelse ; ; Remove missing data and data south of 30N ; keeplist=where(finite(statvals) and (statlat ge 30.),nkeep) if nkeep gt 0 then begin statvals=statvals(keeplist) statlat=statlat(keeplist) statlon=statlon(keeplist) ; boxlat=boxlat(keeplist) ; boxlon=boxlon(keeplist) fracvals=fracvals(keeplist) ; densfd=densfd(keeplist) ; x=x(keeplist) ; y=y(keeplist) endif print,'No. of boxes with data:',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=2.5 & dy=2.5 gnx=360./dx gny=(90.-30.)/dy gx=findgen(gnx)*dx-180. gy=90.-findgen(gny)*dy gridfd=gridit(gnx,gny,gx,gy,statlon,statlat,statvals,fracvals,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) ; levels=findgen(16)*0.4-3. 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 def_1color,r,g,b,10,color='purple' def_1color,r,g,b,16,color='lgreen' def_smearcolor,r,g,b,fromto=[10,16] def_1color,r,g,b,17,color='lgrey' def_1color,r,g,b,18,color='lyellow' def_1color,r,g,b,24,color='red' def_smearcolor,r,g,b,fromto=[18,24] levcol=indgen(15)+10 ; scale_vert,levels=levels,c_colors=levcol ; inter_const,gridfd,gx,gy,map=map,$ ;inter_const,statvals,boxlon,boxlat,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,$ miss_grey='white',$ levels=levels,c_colors=levcol,$ /cell_fill,fdsm=fdsm,$ title='1978 (20 yr) normalised density anomaly' ; contour,fdsm.fd,fdsm.x,fdsm.y,levels=levels,/overplot ; scale_vert,/off ; end