Re: Computing monthly-mean using daily data

From: Christian Pagé <page.christian_at_nyahnyahspammersnyahnyah>
Date: Fri, 8 Feb 2008 17:27:20 +0100

Dennis,

I did the following code which works very well and is quite fast,
once I extracted the time period I wanted before re-ordering the
dimensions....

=================================================
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

;
; Local functions and procedures
;
procedure usage()
begin
  quote = inttochar(34)
  print("")
  print("generate_mean_variance:: USAGE")
  print("Program to compute monthly mean and variance of CF-1.0 convention
NetCDF file.")
  print("Needed variables: lat (2D), lon (2D), x (1D), y (1D), time (1D)")
  print("Fields must be 3D (time,x,y)")
  print("")
  print("Example usage: ncl infile=\" + quote + "T.DAT_france_1985_1986.nc\"
+ quote + " outfile=\" + quote + "out.nc\" + quote + "
generate_mean_variance.ncl")
  print("")
  print("infile=input NetCDF file")
  print("outfile=output NetCDF file")
  print("")
end

;
; Command-line arguments
;
if (.not. isvar("infile")) then ; is infile on command line?
  usage()
  exit
end if
if (.not. isvar("outfile")) then ; is outfile on command line?
  usage()
  exit
end if

; Open files
filein = addfile(infile,"r")

; remove any pre-existing file if necessary
system ("rm -f " + outfile)
fileout = addfile(outfile, "c")

; Get input global attributes
att_names = getvaratts(filein)
if(.not.all(ismissing(att_names))) then
; Assign input global attributes to output
  do i = 0,dimsizes(att_names)-1
    print("Copying Global Attributes: " + att_names(i) + " = " +
filein@$att_names(i)$)
    fileout@$att_names(i)$ = filein@$att_names(i)$ ; copy inputfile
global attributes
  end do
end if

; Get variable names
var_names = getfilevarnames(filein)
numvar = dimsizes(var_names)

time = filein->time ; time:units = "hours
since"
TIME = ut_calendar(time, 0) ; type float
year = floattointeger( TIME(:,0) )
month = floattointeger( TIME(:,1) )
day = floattointeger( TIME(:,2) )
ymd = ut_calendar(time, -2) ; type integer in YYYYMMDD
ymdw = ymd ; mask array when parsing year/month
couple
ymdw@_FillValue = -99
nt = dimsizes(year)
; Number of months in input file
nmonths = nt / 30
; New time vector for monthly data
newtime = new((/nmonths/), integer)
; Copy attributes from old to new time vector
copy_VarAtts(filein->time, newtime)

;
; loop over variables of inputfile and compute mean and variance, then
output to file
;
do i=0,numvar-1

 print("Writing NetCDF output file: i=" + i + " name=" + var_names(i))

; Get number of dimensions of input variable
 rank = dimsizes(dimsizes(filein->$var_names(i)$))

 done = False
 if (var_names(i).ne."x" .and. var_names(i).ne."y" .and.
var_names(i).ne."time") then
; x and y are already generated automatically when writing other dependent
variables, so skip them
; Also skip time because we recreate a new time variable

; Other 2D variables (lat, lon, etc)
   if (rank .eq. 2 .and. done .eq. False) then
     fileout->$var_names(i)$ = filein->$var_names(i)$
     done = True
   end if

; 3D variables : we compute mean and variance over time
   if (rank .eq. 3 .and. done .eq. False) then

; Get dimensions
     dimsiz = dimsizes(filein->$var_names(i)$)
     ny = dimsiz(1)
     nx = dimsiz(2)

; Loop over time and compute monthly mean and variance of field

     t = 0
     mm = 0

; Allocate new variance and mean variables
     ncvar = new((/nmonths,ny,nx/), float)
     ncmean = new((/nmonths,ny,nx/), float)

     do while (t .lt. nt)

; Compute current year and month. Compute also next year and month.
       cury = ymdw(t) / 10000
       curm = (ymdw(t) / 100) - (cury * 100)
       curym = (cury * 100) + curm
       curymn = curym + 1

; Find all indices for the current month/year couple
       curindex = ind( (ymdw .ge. curym*100) .and. (ymdw .lt. curymn*100) )
       if (.not.all(ismissing(curindex)))

; First select only times we need (for MUCH faster reordering for next
step!!!)
         tmpvar = filein->$var_names(i)$(curindex,:,:)
; Then reorder dimensions
         tmpvarr = tmpvar(y|:,x|:,time|:)

; Mask values which will be processed in ymdw so that we don't
; do it more than one time for one month/year couple
         ymdw(curindex) = ymdw@_FillValue

; Compute mean and variance of 2D field
; and store for current month
         ncmean(mm,:,:) = dim_avg_Wrap( tmpvarr ) ; ==> xmean(lat,lon)
         ncvar(mm,:,:) = dim_variance_Wrap( tmpvarr ) ; ==>
xvar(lat,lon)

; Some mandatory cleanup
         delete(tmpvar)
         delete(tmpvarr)

; Set new time vector to number of hours
         newtime(mm) = t*24

; Increment month number
         mm = mm + 1

         done = True

       end if
       delete(curindex)

       t = t + 1

     end do

; Write new mean and variance data to NetCDF
; Set time dimension to new vector
     ncmean!0 = "time"
     ncmean&time = newtime
     ncmean_at_long_name = filein->$var_names(i)$@long_name + "_mean"
     namevar = var_names(i) + "_mean"
     fileout->$namevar$ = ncmean

     ncvar!0 = "time"
     ncvar&time = newtime
     ncvar_at_long_name = filein->$var_names(i)$@long_name + "_variance"
     namevar = var_names(i) + "_variance"
     fileout->$namevar$ = ncvar

; Cleanup
     delete(ncvar)
     delete(ncmean)

   end if

; Other cases (scalar, etc.)
   if (done .eq. False) then
     fileout->$var_names(i)$ = filein->$var_names(i)$
   end if

 end if

end do

print(fileout) ; print overview of contents of output file
=================================================

2008/2/8, Christian Pagé <page.christian_at_gmail.com>:
>
> Thanks Dennis,
>
> I will give it a try and let you know how it goes.
>
> Regards
>
> Christian
>
> 2008/2/7, Dennis Shea <shea_at_ucar.edu>:
> >
> > This should be:
> > i = ind(year(nyr).eq.year .and. month(nmo).eq.(nmo+1) )
> >
> > Dennis Shea wrote:
> > > C/fortran ... any compiled language is generally faster than
> > > interpreted languages. f90 is an array language so the following
> > > is more like f90 than C.
> > >
> > > Completely untested .... You could make a function out of this
> > > if you want.
> > >
> > > nyrs = max(year)-min(year)+1
> > > ntim = nyrs*12
> > >
> > > MONAVG = ((/nyrs,12, nlat,mlon/), typeof(x), getFillValue(x))
> > > or
> > > monavg = ( (/ntim,nlat,mlon/) typeof(x), getFillValue(x))
> > >
> > > nt = -1
> > > do nyr=0,nyrs-1
> > > do nmo=0,11
> > > i = ind(year(nyr).eq.year .and. month(nmo).eq.(nmo)+1)
> > >
> > > MONAVG(nyr,nmo,:,:) = dim_avg(lat|:,lon|:, time|i)
> > > or
> > > nt = nt+1
> > > monavg(nt,:,:) = dim_avg(lat|:,lon|:, time|i)
> > >
> > > delete(i)
> > > end do
> > > end to
> > >
> > >
> > >
> > > Christian Pagé wrote:
> > >> Hello all,
> > >>
> > >> I am still a newbie regarding NCL advanced array functions and
> > >> subscripting (I always think more in a C-like programming fashion).
> > >> It is probably why I cannot find a proper way to do the following
> > >> efficiently using NCL :
> > >>
> > >> Goal : compute monthly-mean using daily data...
> > >>
> > >> I have a CF-1.0 compliant NetCDF file with some daily-data 2D fields
> > >> (time,x,y).
> > >> I also have the time in a 1D vector in hours.
> > >> I want to compute the monthly mean and variance and generate two
> > >> monthly data fields.
> > >>
> > >> I already dealt with time :
> > >> time = filein->time ; time:units = "hours
> > >> since" TIME = ut_calendar(time,
> > >> 0) ; type float
> > >> year = floattointeger( TIME(:,0) )
> > >> month = floattointeger( TIME(:,1) )
> > >> day = floattointeger( TIME(:,2) )
> > >>
> > >> and I have my field, accessing it with : filein->dailyfield
> > >>
> > >> Now, the question is :
> > >> Is there an efficient way to loop over time, and :
> > >> - Find indices which belong to a specific month/year couple
> > >> - Mean my field on all these indices
> > >> - Mark this month/year couple as "done"
> > >>
> > >> as I don't know in advance :
> > >> - the number of years in the input data
> > >> - the number of months for each year in the input data (whole years
> > >> are not in the same files)
> > >>
> > >> This can be done easily in C using standard loops and extra vectors
> > >> to keep track of what has been done.
> > >>
> > >> Many thanks for any hints or solutions...! :)
> > >>
> > >> --
> > >> Christian Pagé
> > >>
> > ------------------------------------------------------------------------
> > >>
> > >> _______________________________________________
> > >> ncl-talk mailing list
> > >> ncl-talk_at_ucar.edu
> > >> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> > >>
> > >
> > >
>
>
>

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Feb 08 2008 - 09:27:20 MST

This archive was generated by hypermail 2.2.0 : Tue Feb 12 2008 - 14:45:27 MST