I often find myself needing to compare datasets with different temporal resolution, like 3-hourly and 6-hourly. This is not a trivial thing to do on the fly, and storing redundant datasets for these occasions is wasteful. So I made an NCL function that averages the leftmost dimension into blocks, effectively degrading the time resolution by whatever factor I want.I call this function “block_avg()” because I’m averaging consecutive bundles of data together. For example, I tend to store climate model data as 6-hourly averages, so the most common use for this is to convert 6-hourly data into daily means. Since there are four 6-hour time samples for each day, the actual syntax to use my function would simply be:
X_daily = block_avg(X_6hour, 4)
The great thing about this function is that it’s really fast. The secret to making it fast is to not do any big loops. Instead I simply loop over the size of the block. Then by using the “stride” syntax I can pull all data associated with the n-th part of the block (e.g. x(n: : stride) ).
In the example of converting 6-hourly to daily, I first find all the data associated with hour 00UTC. I then find all the data associated with 06UTC and perform a running average with the 00UTC data that I previously isolated. I repeat this two more times for the 18UTC and 21UTC data.
I haven’t figured out a clean way to do this for an n-dimensional input, so far I only go up to 5 dimensions. If you have an idea how to generalize this for n-dimensions please let me know!
;======================================================================= ;======================================================================= ; This function averages the leftmost dimension into ; blocks of size 'nblock'. the resulting size of the ; leftmost dimensions of the output is odim(0)/nblock ; Useful for reducing the temporal resolution of data. ; Ex. converting 6-hourly data into daily averages ; X_daily = block_avg(X_6hour,4) ; Created by Walter Hannah ; University of Miami Sept. 2014 undef("block_avg") function block_avg(Xin,nblock) local Xin,nblock,Xout,tmp,odim,rank,n,cnt,tcnt begin odim = dimsizes(Xin) ; Output dimension sizes rank = dimsizes(odim) ; # of dimensions ;------------------------------------------- ; check that nblock is evenly divisible if odim(0)%nblock.ne. 0 then print("block_avg - Error: dimension 0 ("+odim(0)+") \ not divisible by nblock ("+nblock+")") exit end if ;------------------------------------------- ; average data in blocks of the left dimension odim(0) = odim(0)/nblock Xout = new(odim,typeof(Xin)) do n = 0,nblock-1 if rank.eq.5 then tmp = Xin(n::nblock,:,:,:,:) end if if rank.eq.4 then tmp = Xin(n::nblock,:,:,:) end if if rank.eq.3 then tmp = Xin(n::nblock,:,:) end if if rank.eq.2 then tmp = Xin(n::nblock,:) end if if rank.eq.1 then tmp = Xin(n::nblock) end if if n.eq.0 then cnt = tmp*0. end if inc = where(.not.ismissing(tmp),1.,0.) tmp = where(.not.ismissing(tmp),tmp,0.) cnt = cnt+inc tcnt = where(cnt.ne.0,cnt,1.) if n.eq.0 then Xout = tmp end if if n.ne.0 then Xout = ( tmp + (cnt-inc)*Xout )/tcnt end if end do ;------------------------------------------- ; Copy the metadata copy_VarAtts(Xin,Xout) dims = getvardims(Xin) do n = 0,dimsizes(dims)-1 if .not.ismissing(dims(n)) then Xout!n = dims(n) if iscoord(Xout,dims(n)) then if n.eq.0 then Xout&$dims(n)$ = Xin&$dims(n)$(::nblock) end if if n.gt.0 then Xout&$dims(n)$ = Xin&$dims(n)$ end if end if end if end do ;------------------------------------------- return Xout end ;======================================================================= ;=======================================================================