function gridit,nx,ny,xlon,ylat,statlon,statlat,statval,statweight,$ nstat=nstat,statflag=statflag,minv=minv,maxv=maxv ; ; **** THIS VERSION HAS A WEIGHTING FOR EACH STATION, THAT IS USED WHEN ; AVERAGING A NUMBER OF STATIONS INTO A SINGLE BOX **** ; ### AND GIVES MIN AND MAX IN EACH CELL ### ; Returns a gridded result from a list of station values ; where the grid is defined by xlon and ylat (centre of boxes) ; and where the stations are defined by statlon and statlat. ; If statval is 2D, its first dimension represents a sequence of fields to ; create. The number of stations contributing to each grid-box can also ; be returned if required. If more than one station contributes to a box, ; then the values are simply averaged. If no stations in a box, NaN is ; returned. Missing data must be stored as NaN to be accounted for. ; If statflag is given, then any station with a non-zero flag is ignored. ; ;------------------------------------------------------------------- ; ; Check dimensions, and ensure a 2D input array ; n1=n_elements(statlon) n2=n_elements(statlat) statsize=size(statval) ndim=statsize(0) if (ndim lt 1) and (ndim gt 2) then message,'statval must be 1D or 2D' if ndim eq 1 then begin ; print,'1D input' nt=1 ns=statsize(1) stat2d=reform(statval,nt,ns) ; print,nt,ns ; help,stat2d endif else begin ; print,'2D input' nt=statsize(1) ns=statsize(2) stat2d=statval ; print,nt,ns ; help,stat2d endelse if (ns ne n1) or (ns ne n2) then message,'Incompatible station arrays' if n_elements(statweight) eq 0 then statweight=fltarr(ns)+1. ;print,statweight ; ; Define arrays ; fdnew=dblarr(nt,nx,ny) fdmin=fltarr(nt,nx,ny) & fdmin(*,*,*)=1.E10 fdmax=fltarr(nt,nx,ny) & fdmax(*,*,*)=-1.E10 nstat=fltarr(nt,nx,ny) wstat=dblarr(nt,nx,ny) ; ; Compute grid-box edges ; dx=xlon(1)-xlon(0) xlonedge=[xlon(0)-0.5*dx,xlon+0.5*dx] ;print,'xlonedge' ;print,xlonedge dy=ylat(1)-ylat(0) ylatedge=[ylat(0)-0.5*dy,ylat+0.5*dy] ;print,'ylatedge' ;print,ylatedge ; ; Start filling the arrays ; for istat = 0 , ns-1 do begin ; print,'Station:',istat ; ; Ignore it if requested to do so ; iwant=1 if n_elements(statflag) gt 0 then begin if statflag(istat) ne 0 then iwant=0 endif if iwant eq 1 then begin ; print,'Using it' ; ; Locate it's grid-box ; xlist=where( (xlonedge(0:nx-1) lt statlon(istat)) and $ (xlonedge(1:nx) ge statlon(istat)),nbox1) ylist=where( (ylatedge(0:ny-1) ge statlat(istat)) and $ (ylatedge(1:ny) lt statlat(istat)),nbox2) if (nbox1 gt 1) or (nbox2 gt 1) then message,'Oops!' if (nbox1 eq 1) and (nbox2 eq 1) then begin ; print,'Location:',xlon(xlist),ylat(ylist) ; ; Check each time value separately ; for itime = 0 , nt-1 do begin if finite(stat2d(itime,istat)) then begin fdnew(itime,xlist,ylist)=fdnew(itime,xlist,ylist)+$ stat2d(itime,istat)*statweight(istat) nstat(itime,xlist,ylist)=nstat(itime,xlist,ylist)+1. wstat(itime,xlist,ylist)=wstat(itime,xlist,ylist)+statweight(istat) fdmin(itime,xlist,ylist)=min([fdmin(itime,xlist,ylist),$ stat2d(itime,istat)]) fdmax(itime,xlist,ylist)=max([fdmax(itime,xlist,ylist),$ stat2d(itime,istat)]) endif endfor ; endif ; endif ; endfor ; ; Average the totals ; for i = 0 , nx-1 do begin for j = 0 , ny-1 do begin for itime = 0 , nt-1 do begin if nstat(itime,i,j) gt 0. then begin fdnew(itime,i,j)=fdnew(itime,i,j)/wstat(itime,i,j) if nstat(itime,i,j) eq 1. then begin fdmin(itime,i,j)=!values.f_nan fdmax(itime,i,j)=!values.f_nan endif endif else begin fdnew(itime,i,j)=!values.f_nan fdmin(itime,i,j)=!values.f_nan fdmax(itime,i,j)=!values.f_nan endelse endfor endfor endfor ; nstat=reform(nstat) minv=reform(fdmin) maxv=reform(fdmax) return,float(reform(fdnew)) ; end