Running Mean

참고: here

 

A Running average filter The similarity between NCL and fortran-90 syntax can readily be seen in this example. Note that an NCL built-in function (runave) exists that can perform this calculation. The example is still useful however in demonstrating function construction. The following built-in NCL functions are used: dimsizes(to check the number of elements in a vector), isatt (to check to for an attribute), new (dynamic memory allocation), all (test all elements of a vector/array; true or false), ismissing (test for missing value in a vector or array; true or false), sum (arithmetic sum of elements), delete (deallocate memory), and return (to return a scalar or vector/array). ;========================NCL Function: runave=========================== function runave(x[*]:float, nptav:integer, kopt:integer)   ; note the 'typing' ; this routine will perform a running average on x                               ; .   x may contain missing data ; ; usage: ;                 y = runave(x,nptav,kopt)      ; y returned ; nomenclature:                                                                 ; .   x         - on input the series to be smoothed                             ; .   nptav     - no. of pts to be utilized in the running average               ; .               nptav must be odd  and nptav >= 3                             ; .   kopt      - end-point option                                               ; .               kopt < 0 : utilize cyclic conditions                           ; .                          e.g.,  nptav=3 at n=1                               ; .                          x(1) = (x(nptx)+x(1)+x(2)) / 3.                     ; .               kopt = 0 : set un-smoothed beginning and end pts to spvl       ; .                          e.g.,  nptav=3 , x(1) = spvl and x(nptx) = spvl     ; .               kopt > 0 : utilize reflective (symmetric) conditions           ; .                          e.g.,  nptav=3 at n=1                               ; .                          x(1) = (x(2)+x(1)+x(2)) / 3.                       local nptx,w,nptw,n,na,rnptav,nxb,nxu,nwb,nwu begin       nptx = dimsizes(x)                      ; number of elements in vector "x"         if (nptav.gt.nptx) then                 return(-16)         else                 if ((nptav%2).eq.0) then      ; nptav must be odd;  0s "mod"                         return(-17)                 else                         if ((isatt(x,"_FillValue")).and.(all(ismissing(x)))) then                                 return(-18)                         end if                 end if         end if       nptw = nptx + (nptav-1)               ; work vector: extra values for end pts       if (isatt(x,"_FillValue")) then           w = new ( (/nptw/), float, x@_FillValue)           w@_FillValue = x@_FillValue       ; use x FillValue       else           wmsg = -999.0                     ; local FillValue           w = new ( (/nptw/), float, wmsg)           w@_FillValue = wmsg                       end if                                                    ; remember NCL is zero based       na  = nptav/2                          ; subscript offset       nxb = 0                                ; base  dim of x       nxu = nptx-1                           ; upper dim of x       nwb = na                               ; w(nwb) = x(0)       nwu = nxu+na                           ; w(nwu) = x(nxu)       w(nwb:nwu) = x(:)                      ; fill middle of w                         if (kopt.lt.0) then           w(0:na-1)       = x(nptx-na:nxu)           w(nwu+1:nwu+na) = x(0:na-1)       else         if (kopt.eq.0) then           w(0:na-1)       = w@_FillValue           w(nwu+1:nwu+na) = w@_FillValue               else              ; must be kopt>0           w(0:na-1:-1)    = x(1:na)           w(nwu+na:nwu+1:-1) = x(nxu-na:nxu-1:-1)               end if       end if       y = new ( (/nptx/),float,x@_FillValue); dynamic allocation return vector       do n=0,nptx-1          if (all(.not.ismissing(w(n:n+nptav-1) ))) then              y(n) = sum(w(n:n+nptav-1))/nptav ; float returned          end if       end do       delete (w)                            ; delete work vector; no longer needed       return (y)                            ; return vector end The Fortran-90 subroutine code which performs the same task is: ;========================F90 Function: runave===========================       subroutine runave_90 (x,npts,nptav,xmsg,kopt,ier)       implicit  none       integer , intent(in)     :: npts, nptav, kopt       real    , intent(in)     :: xmsg       integer , intent(inout)  :: ier       real    , intent(inout)  :: x(npts)                                                                                 ! this routine will perform a running average on x                               ! .   x may containg missing data ! .   the end points may be treated in three different ways                                                                                                     ! nomenclature:                                                                 ! .   x         - on input the series to be smoothed                             ! .               on output the smoothed series                                 ! .   npts      - no of pts in x                                                 ! .   nptav     - no. of pts to be utilized in the runaing average               ! .               nptav must be odd.                                             ! .   kopt      - end-point option                                               ! .               kopt < 0 : utilize cyclic conditions                           ! .                          e.g.,  nptav=3 at n=1                               ! .                          x(1) = (x(npts)+x(1)+x(2)) / 3.                     ! .               kopt = 0 : set unsmoothed beginaing and end pts to spvl       ! .                          e.g.,  nptav=3 , x(1) = spvl and x(npts) = spvl     ! .               kopt > 0 : utilize reflective (symmetric) conditions           ! .                          e.g.,  nptav=3 at n=1                               ! .                          x(1) = (x(2)+x(1)+x(2)) / 3.                       ! .   xmsg      - if kopt=0 set end points to this special value                 ! .               xmsg = 0. or 1.e+36 or the series mean are common             ! .   ier       - error flag                                                     ! .               ier = 0 : no errors                                           ! .               ier = -11 : npts @ 0                                           ! .               ier = -12 : npts < nptav                                       ! .               ier = -13 : nptav is not odd                                         real,   allocatable, dimension(:) :: work       integer n, na, np1, nm1       ier = 0                                                                         if (npts.lt.1) ier = -11                                                         if (mod(nptav,2).eq.0) ier = -13                                                 if (nptav.gt.npts) ier = -12                                                     if (ier.lt.0) return                                                             allocate (work(1-nptav/2:npts+nptav/2), stat=ier)       if (ier.ne.0) then           return       endif       work(1:npts) = x(1:npts)    ! vector copy                                         ! handle end points       na  = nptav/2-1             ! subscript offset       np1 = npts+1                ! temporary       nm1 = npts-1                ! temporary       if (kopt.lt.0) then           work(-na:0)      = x(npts-na:npts)           work(np1:np1+na) = x(1:1+na)       elseif (kopt.eq.0) then           work(-na:0)      = xmsg           work(np1:np1+na) = xmsg       else           work(0:-na:-1)   = x(2:2+na)           work(np1:np1+na) = x(nm1:nm1-na:-1)       endif       x  = xmsg                   ! preset to xmsg       na = nptav/2       do n=1,npts          if (all(work(n-na:n+na) /= xmsg) ) then              x(n) = sum( work(n-na:n+na) )/real(nptav)          endif       enddo       deallocate (work,  stat=ker)       return       end