NCL: Averaging in blocks (i.e. converting to longer time steps)

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
;=======================================================================
;=======================================================================

 

 

Leave a Reply

Your email address will not be published. Required fields are marked *