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