; ; Need to provide a list of years in vector yrlist ; In fact, yrlist is 2D array (n,m), with n different maps each made up ; of the mean of m different years. If m equal NaN then it is ignored. ; titlist=list of plot titles ; if n_elements(yrlist) eq 0 then message,'Need list of years in yrlist' ; ; Get display ready ; loadct,39 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 def_1color,r,g,b,40,color='brightblue' def_1color,r,g,b,49,color='paleblue' def_1color,r,g,b,50,color='palered' def_1color,r,g,b,59,color='brightred' def_smearcolor,r,g,b,fromto=[40,49] def_smearcolor,r,g,b,fromto=[50,59] ; ; 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 ; lsz=size(yrlist) nmap=lsz(1) print,lsz print,nmap for i = 0 , nmap-1 do begin alist=reform(yrlist(i,*)) onelist=alist( where( finite(alist) ) ) nyr=n_elements(onelist) for iyy = 0 , nyr-1 do begin iyr=onelist(iyy) print,iyy,iyr onefd=rd1yr(iyr,x=x,y=y,frac=frac) if iyy eq 0 then begin nchron=n_elements(onefd) allfd=fltarr(nchron,nyr) endif allfd(*,iyy)=onefd endfor misslist=where(allfd eq -9.99,nmiss) if nmiss gt 0 then allfd(misslist)=!values.f_nan ; make sure the fields to be totalled are 2D even if only one year is used! summap=total(reform(allfd,nchron,nyr),/nan,2) coumap=total(reform(finite(allfd),nchron,nyr),2) densfd=summap/float(coumap) ; ; Remove missing data and data south of 30N ; keeplist=where(finite(densfd) 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,'No. of points 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=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=[-20.,findgen(19)-9.,20.] levels=levels*0.5 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=1.5*fac & if nzero gt 0 then c_thick(zeroloc)=1. c_colors=findgen(20)+40. ; inter_const,gridfd,gx,gy,map=map,$ maxdist=2500.,$ gs=[2.,2.],$ sm=sm,labels=labels,$ miss_grey='white',$ c_colors=c_colors,levels=levels,$ /cell_fill,$ title=titlist(i),fdsm=fdsm contour,fdsm.fd,fdsm.x,fdsm.y,levels=levels,c_labels=c_labels,$ c_charsize=c_charsize,/follow,c_thick=c_thick,/overplot ; endfor ; end