<gh_stars>1-10
#!/usr/bin/env python

# Incremental data reduction.

# Suppose you have model output in files file1,file2,...,fileN.
# All files are the same except that they cover different times (in sequence; together they cover
# all time without overlapping); and of course variable values are different.
# Suppose you have reduced the data, e.g. time from daily to monthly averages, maybe also some space
# reductions and maybe limited to a few variables.  The reduced data is in a file reducf.
# Now we get the next file fileN+1.  Our goal here is to update reducf, without re-opening the
# first N files.  Of course, there is the initial case N=0, where reducf doesn't exist.

# This script, when finished, will also be able to compute climatologies.  Just replace each
# model output file's time units from "...since 1850-00-00" to "...since 0000-00-00" which
# would be the time units of a climatology file.  And shift all its time values by an integer
# number of years, so that all time values fall within the same year (time bounds may slip
# into another year.  There may be an issue here, especially if the time isn't centered
# within its bounds.)

# It will be essential that the time axis be present, with a bounds attribute.
# (However, the output file will also have a climatology attribute on the time axis.)
# Time averaging will be weighted by the size of the time interval.  Thus, for the noleap calendar
# with units in days, January would have a weight of 31 relative to the February weight of 28.
# If units were months, each month would have the same weight.

# On machines, such as rhea.ccs.ornl.gov, which impose a very large random cost to opening
# a file, this script will perform better than one which is based on cdms datasets and cdscan.

# The partial time reduction depends, obviously, on what we use for target time intervals.
# The most general way to specify them is to provide a function, so that will be an option.
# The general case is complicated because the intervals need not be equal.
# But usually users will want "annual", "monthly", or "seasonal climatologies".
# I will provide functions for these cases.
# Note that "annual" intervals are unequal in a Gregorian calenar with leap years, and
# "monthly" intervals are unequal except in a 360-day calendar.
# A keyword + calendar + units will be sufficient if it's one of the likely requests.

# >>>> TO DO:
# Create a suitable time axis when reading climo data without one.

import numpy, cdms2, sys, os, math, logging
from numbers import Number
from pprint import pprint
import time
from metrics.packages.acme_regridder.scripts.acme_regrid import addVariable
import logging

logger = logging.getLogger(__name__)

#import debug

# Silence annoying messages about setting the NetCDF file type.  Also, these three lines will
# ensure that we get NetCDF-3 files.  Expansion of a FileAxis may not work on NetCDF-4 files.
cdms2.setNetcdfShuffleFlag(0)
cdms2.setNetcdfDeflateFlag(0)
cdms2.setNetcdfDeflateLevelFlag(0)

def daybounds( season ):
    """Input is a 3-character season, e.g. JAN, JJA, ANN.  Output is bounds for that season
    in units of "days since 0".
    In most cases the bounds will be like [[lo,hi]].  But for DJF it will be like
    [[lo1,hi1],[lo2,hi2]].
    The noleap (365 day) calendar is required."""
    # This is an incredibly crude way to do it, but works for now.
    # If permanent, I would make these dictionaries global for performance.
    dbddic = { 'JAN': [[0,31]],
               'FEB': [[31,59]],
               'MAR': [[59,90]],
               'APR': [[90,120]],
               'MAY': [[120,151]],
               'JUN': [[151,181]],
               'JUL': [[181,212]],
               'AUG': [[212,243]],
               'SEP': [[243,273]],
               'OCT': [[273,304]],
               'NOV': [[304,334]],
               'DEC': [[334,365]] }
    sesdic = { 'MAM': [[ dbddic['MAR'][0][0], dbddic['MAY'][0][1] ]],
               'JJA': [[ dbddic['JUN'][0][0], dbddic['AUG'][0][1] ]],
               'SON': [[ dbddic['SEP'][0][0], dbddic['NOV'][0][1] ]],
               'DJF': [ dbddic['DEC'][0],
                        [dbddic['JAN'][0][0], dbddic['FEB'][0][1]] ],
               'ANN': [[ dbddic['JAN'][0][0], dbddic['DEC'][0][1] ]] }
    dbddic.update( sesdic )
    return dbddic[season]

def setClimoBounds( time, tfile=None, value=None ):
    """time is a TransientAxis or FileAxis.  tfile is the open file it came from, value a new value
    for its bounds.
    The _bounds_ attribute will be set for future use, as long as the axis is a time axis.
    This will be set from the provided value if any, otherwise will come from files: from the
    time.climatology variable if possible; otherwise the regular getBounds() will be called.
    If both a file and a value are provided, then also the climatology or bounds variables in the
    file will be set from the provided value."""
    if not time.isTime():
        return
    if value is None and tfile is not None:
        # Get _bounds_ from tfile.
        if hasattr( time, 'climatology' ):
            bounds = tfile( time.climatology )[:,:]
            if getattr(bounds,'initialized','yes')=='no':
                time._bounds_ = time.getBounds()
            else:
                time._bounds_ = bounds
        else:
            time._bounds = time.getBounds()
    elif value is None:  # tfile is None
        if hasattr(time,'climatology') and hasattr(time,'parent') and\
                time.parent.__class__.__name__=='CdmsFile':
            setClimoBounds( time, time.parent )
        else:
            time._bounds_ = time.getBounds()
    else:   # value is not None
        # Set _bounds_ from provided value.
        time.setBounds(value)
        time._bounds_ = value  # needed to set _bounds_ for FileVariable
        if tfile is not None:
            # Set bounds variables in tfile.
            if hasattr( time, 'climatology' ):
                f( time.climatology )[:,:] = value
            if hasattr( time, 'bounds' ):
                f( time.bounds )[:,:] = value

    return

def getClimoBounds( time ):
    """time is a TransientAxis or FileAxis, and setClimoBounds() should have already been
    called on it.  This function will return its bounds, from the _bounds_ attribute."""
    if not hasattr( time, '_bounds_' ):
        setClimoBounds( time )
    return time._bounds_

def next_tbounds_copyfrom_data( red_time_bnds, data_time_bnds ):
    """This is one of many ways to provide a time interval for the computing a partial time
    reduction of data.  This one works for no real time reduction.  A typical case is that
    we have a sequence of monthly averages and the 'reduced' time is monthly averages.
    The idea is simply that if the data time bounds are greater than the reduced time bounds,
    use the data time bounds as the next interval in the reduced time sequence.
    red_time_bnds and data_time_bnds may be treated as arrays like a typical time_bnds variable,
    e.g. numpy.array([[1234,1244],[1244,1256],[1256,1267]]).  In practice, red_time_bnds is
    normally a FileVariable with such a value.
    You get an initial value for red_time_bnds by calling this function with red_time_bnds==[].
    """
    if len(data_time_bnds)==0:
        return red_time_bnds
    if len(red_time_bnds)==0:
        return data_time_bnds
    if data_time_bnds[0][0]>=red_time_bnds[-1][1]:
        return numpy.concatenate( ( red_time_bnds, data_time_bnds ), axis=0 )
    else:
        return red_time_bnds

def next_tbounds_prescribed_step( red_time_bnds, data_time_bnds, dt ):
    """The next bounds will be a prescribed step above the existing one.
    This uses data_time_bnds only to use its lowest bound in initialization of red_time_bnds.
    This is intended to be called from another function, e.g. a lambda expression which includes
    a fixed step size."""
    # Normally the numbers in the return value are an integers between 0 and 365,
    # representing a time in days
    if len(red_time_bnds)==0:
        t0 = data_time_bnds[0][0]
        return numpy.array([[t0,t0+dt]])
    t0 = red_time_bnds[-1][0]
    t1 = red_time_bnds[-1][1]
    t1t2 = numpy.array([[t0+dt, t1+dt]],dtype=numpy.int32)
    return numpy.concatenate( ( red_time_bnds, t1t2 ), axis=0 )

def initialize_redfile( filen, axisdict, typedict, attdict, varnames, lock=None, dt=-1 ):
    """Initializes the file containing reduced data.  Input is the filename, axes
    variables representing bounds on axes, and names of variables to be initialized.
    Axes are represented as a dictionary, with items like varname:[timeaxis, lataxis, lonaxis].
    Time should always be present as the first axis.
    Note that cdms2 will automatically write the bounds variable if it writes an axis which has bounds.
    """
    if lock is not None:  lock.acquire()
    #t1 = time.time()
    g = cdms2.open( filen, 'w' )
    #t2 = time.time()
    if lock is not None:  lock.release()
    if 'time_bnds' not in varnames:
        addVariable( g, 'time_bnds', typedict['time_bnds'], axisdict['time_bnds'], attdict['time_bnds'] )
    if 'climatology_bnds' not in varnames and dt==0:
        addVariable( g, 'climatology_bnds', typedict['time_bnds'], axisdict['time_bnds'], attdict['time_bnds'] )
    for varn in varnames:
        addVariable( g, varn, typedict[varn], axisdict[varn], attdict[varn] )
    addVariable( g, 'time_weights', typedict['time'], [axisdict['time']], [] )
    g['time_weights'].initialized = 'no'
    return g;

def initialize_redfile_from_datafile( redfilename, varnames, datafilen, dt=-1, init_red_tbounds=None,
                                      force_double=False, lock=None ):
    """Initializes a file containing the partially-time-reduced average data, given the names
    of the variables to reduce and the name of a data file containing those variables with axes.
    The file is created and returned in an open state.  Closing the file is up to the calling function.
    Also returned is a list of names of variables which have been initialized in the file.
    If provided, dt is a fixed time interval for defining the time axis.
    dt=0 is special in that it creates a climatology file.
    If provided, init_red_tbounds is the first interval for time_bnds.
    For example in the 360-day calendar with time units of days, dt=360 and
    init_red_tbounds=[[150,240]] would set up the partially-time-reduced averages for annual
    trends of values averaged over the JJA season.
    """
    # Note: some of the attributes are set appropriately for climatologies, but maybe not for
    # other time reductions.  In particular, c.f. time long_name and varn cell_methods.
    boundless_axes = set([])
    if lock is not None:  lock.acquire()
    #t1 = time.time()
    f = cdms2.open( datafilen )
    #t2 = time.time()
    if lock is not None:  lock.release()
    axisdict = {}
    typedict = {}
    attdict = {}
    timeaxis = None
    if dt>=0:   # time axis should have fixed intervals, size dt
        if init_red_tbounds is None or init_red_tbounds==[] or init_red_tbounds==[[]]:
            init_red_tbounds=[]
            red_time_bnds = next_tbounds_prescribed_step(
                init_red_tbounds, f.getAxis('time').getBounds(), dt )
        else:
            red_time_bnds = init_red_tbounds
        red_time = [ 0.5*(tb[0]+tb[1]) for tb in red_time_bnds ]
        # tb[?] can be expected to be an integer,  red_time 64-bit float.
        timeaxis = cdms2.createAxis( red_time, red_time_bnds, 'time' )
        # ... this will set the _bounds_ attribute to red_time_bnds by calline timeaxis.setBounds().
        timeaxis.bounds = 'time_bnds'
        timeaxis.climatology = 'climatology_bnds'
        axisdict['time'] = timeaxis
        typedict['time'] = f['time'].typecode()
        attdict['time'] = { 'bounds':'time_bnds', 'long_name':'climatological_time' }
        if dt==0:
            attdict['time']['climatology'] = 'climatology_bnds'
        if force_double and typedict['time']=='f':
            # I think that time is always double, but let's make sure:
            typedict['time']='d'
            new_dtype = numpy.float64
            if '_FillValue' in attdict['time']:
                attdict['time']['_FillValue'] = attdict['time']['_FillValue'].astype( new_dtype )
                if 'missing_value' in attdict['time']:
                    attdict['time']['missing_value'] =\
                        attdict['time']['missing_value'].astype( new_dtype )

    out_varnames = []
    for varn in varnames:
        if f[varn] is None:
            logger.warning("%s was not found in %s", varn, datafilen)
            continue
        try:
            dtnom = f[varn].dtype.name  # works for variables
        except:
            # Probably an axis.  Most of what we're doing here makes no sense for an axis.
            # FWIW, this would get the type name: dtnom = f[varn].dtype('spam').name
            logger.warning("WARNING, averager will ignore axis %s.", varn)
            continue
        if dtnom.find('string')==0:
            # We can't handle string variables.
            logger.warning("Ignoring string variable %s", varn)
            continue
        out_varnames.append(varn)
        dataaxes = f[varn].getAxisList()
        axes = []
        for ax in dataaxes:
            if dt>=0 and ax.isTime():
                axes.append( timeaxis )
            else:
                axes.append( ax )
            if not hasattr( ax,'bounds' ): boundless_axes.add(ax)
        axisdict[varn] = axes
        typedict[varn] = f[varn].typecode()

        # attributes:
        vdict = f[varn].__dict__
        attdict[varn] = {a:vdict[a] for a in vdict if (a=='_FillValue' or a[0]!='_')
                         and a!='parent' and a!='autoApiInfo' and a!='domain' and a!='vwgts'}
        # ...It would be worth reconstructing the domain attribute with the new axes.
        if 'cell_methods' not in attdict[varn]:
            attdict[varn]['cell_methods'] = 'time: mean'
        elif 'time: mean' not in attdict[varn]['cell_methods']:
            attdict[varn]['cell_methods'] += ', time: mean'
        attdict[varn]['initialized'] = 'no'

        if force_double and typedict[varn]=='f':
            typedict[varn]='d'
            new_dtype = numpy.float64
            if '_FillValue' in attdict[varn]:
                attdict[varn]['_FillValue'] = attdict[varn]['_FillValue'].astype( new_dtype )
                if 'missing_value' in attdict[varn]:
                    attdict[varn]['missing_value'] =\
                        attdict[varn]['missing_value'].astype( new_dtype )

    for ax in boundless_axes:
        logger.warning("Axis %s has no bounds", ax.id)
    if timeaxis is not None:
        tbndaxis = cdms2.createAxis( [0,1], None, 'tbnd' )

        axisdict['time_bnds'] = [timeaxis,tbndaxis]
        typedict['time_bnds'] = timeaxis.typecode() # always 'd'
        attdict['time_bnds'] = {'initialized':'no'}

        if dt==0:
            out_varnames.append('climatology_bnds')
            axisdict['climatology_bnds'] = [timeaxis,tbndaxis]
            typedict['climatology_bnds'] = timeaxis.typecode() # always 'd'
            attdict['climatology_bnds'] = {'initialized':'no'}

    g = initialize_redfile( redfilename, axisdict, typedict, attdict, out_varnames,
                            lock=lock, dt=dt )

    if timeaxis is not None:
        g['time_bnds'][:] = red_time_bnds
        g['time_bnds'].initialized = 'yes'
        g['time_weights'][:] = numpy.zeros(len(timeaxis))
        g['time_weights'].initialized = 'yes'
    ftime = f.getAxis('time')
    if hasattr(ftime,'climatology'):
        tb = f( ftime.climatology ).data
    else:
        tb = ftime.getBounds()
    if tb is None:
        tmin = f.getAxis('time')[0]
        tmax = f.getAxis('time')[-1]
    else:
        tmin = tb[0][0]
        tmax = tb[-1][-1]
    if lock is not None:  lock.acquire()
    f.close()
    if lock is not None:  lock.release()
    return g, out_varnames, tmin, tmax

def adjust_time_for_climatology( newtime, redtime ):
    """Changes the time axis newtime for use in a climatology calculation.  That is, the values of
    the time axis are shifted to lie within a single year suitable for averaging with the existing
    reduced (climatology) variables redvars.
    """
    newtime.toRelativeTime( redtime.units, redtime.getCalendar() )
    assert( newtime.calendar=='noleap' )
    # btw, newtime.getCalendar()==4113 means newtime.calendar=="noleap"

    # We need to adjust the time to be in the year 0, the year we use for climatology calculations.
    # It would be natural to compute the adjustment from newtime[0].  But sometimes it's equal to
    # one of its bounds, on the border between two years!  So it's better to use the midpoint
    # between two bounds.
    # Usually we can assume that the time bounds are on monthly borders.  This should be checked, but isn't yet.<<<<
    timebnds = getClimoBounds(newtime)
    midtime = 0.5*(timebnds[0][0]+timebnds[0][1])
    N = math.floor(midtime/365.)   # >>>> assumes noleap calendar, day time units!
    newtime[:] -= N*365               # >>>> assumes noleap calendar, day time units!
    setClimoBounds( newtime, value=(getClimoBounds(newtime) - N*365) ) # >>>> assumes noleap calendar, day time units!
    return newtime

def two_pt_avg( mv1, mv2, i, a2, sw1, sw2, aw2=None, force_scalar_avg=False ):
    """Input:
    FileVariables mv1,mv2;  i, (time) index to subset mv1 as a1=mv1[i], a
    corresponding array a2 from mv2, also formed by fixing a point on the first, time, axis;
    scalar weights sw1,sw2.
    Output:   weighted average a of a1 and a2.
    The weights will be scalar weights unless an array weight must be used because either:
    an array weight is already an attribute of mv1 or mv2; or
    v1 and v2 do not have the same mask.
    If an array weight is used, it will become an attribute of the first variable mv1,
    and take its shape.  Normally mv1 is going to be where averages get accumulated,
    so if there is more than one time, it will appear in mv1 but not the input data mv2.
    The optional argument aw2 is meant to support array weights in the case where we are
    averaging different parts of the same variable, i.e. mv1=mv2 but a1==mv1[i]!=a2.
    If aw2 be supplied, it will be used as the array weight for mv2,a2.
    The optional argument force_scalar_avg argument is for testing.  If set to True,
    all computations involving array weights will be bypassed and simple scalar
    weights used instead.
    """
    a1 = mv1[i]
    if not force_scalar_avg:
        w1 = sw1
        if aw2 is not None:
            w2 = aw2
        else:
            w2 = sw2

        # If array weights already exist for mv1 _or_ mv2, use them, and create array weights
        # for whatever doesn't have them.  If mv1,mv2 have different masks, then we need array
        # weights now and henceforth; create them.  The shape attribute is a convenient way
        # to detect whether w1 or w2 is scalar or array.
        if hasattr( mv1, 'vwgts' ):
            f1 = mv1.parent # the (open) file corresponding to the FileVariable mv1
            w1 = f1(mv1.vwgts)
        if hasattr( mv2, 'vwgts' ) and aw2 is None:
            f2 = mv2.parent # the (open) file corresponding to the FileVariable mv2
            if f2 is None:
                f2 = mv2.from_file # the (open) file corresponding to the TransientVariable mv2
            w2 = f2(mv2.vwgts)
        if (not hasattr(w1,'shape') or len(w1.shape)==0) and hasattr(w2,'shape') and len(w2.shape)>0:
            w1 = numpy.full( mv1.shape, -1 )
            w1[i] = sw1
            w1 = numpy.ma.masked_less(w1,0)
        if (not hasattr(w2,'shape') or len(w2.shape)==0) and hasattr(w1,'shape') and len(w1.shape)>0:
            w2 = numpy.full( mv1.shape, -1 )   # assumes mv1 time axis is >= mv2 time axis
            w2[i] = sw2
            w2 = numpy.ma.masked_less(w2,0)
        if (not hasattr(w1,'shape') or len(w1.shape)==0) and\
                (not hasattr(w2,'shape') or len(w2.shape)==0):
            mask1 = False
            mask2 = False
            if hasattr(mv1,'_mask'):
                mask1 = mv1.mask
            elif getattr(mv1,'missing_value'):  # if no missing_value, then no mask
                valu = mv1.getValue()
                if hasattr( valu, '_mask' ):
                    mask1 = valu._mask
            if hasattr(mv2,'_mask'):
                mask2 = mv2.mask
            elif getattr(mv2,'missing_value'):  # if no missing_value, then no mask
                valu = mv2.getValue()
                if hasattr( valu, '_mask' ):
                    mask2 = valu._mask
            if not numpy.all(mask1==mask2):
                # Note that this test requires reading all the data.  That has to be done anyway to compute
                # the average.  Let's hope that the system caches well enough so that there won't be any
                # file access cost to this.
                #jfp was w1 = numpy.full( a1.shape, w1 )
                #jfp was w2 = numpy.full( a2.shape, w2 )
                w1 = numpy.full( mv1.shape, -1 )
                w1[i] = sw1
                w1 = numpy.ma.masked_less(w1,0)
                w2 = numpy.full( mv1.shape, -1 )
                w2[i] = sw2
                w2 = numpy.ma.masked_less(w2,0)

    if not force_scalar_avg and\
            hasattr(w1,'shape') and len(w1.shape)>0 and hasattr(w2,'shape') and len(w2.shape)>0:
        if w1[i].mask.all():   # if w1[i] is all missing values:
            w1[i] = sw1
        if w2[i].mask.all():   # if w2[i] is all missing values:
            w2[i] = sw2
        # Here's what I think the numpy.ma averager does about weights and missing values:
        # The output weight w(i)=sw1+sw2 if there be no missing value for mv1(i) and mv2(i) (maybe
        # also if both be missing, because the weight doesn't matter then).
        # But if i be missing for mv1, drop sw1, thus w(i)=sw2.  If i be missing for mv2, drop sw2.
        a,w = numpy.ma.average( numpy.ma.array((a1,a2)), axis=0,
                                weights=numpy.ma.array((w1[i],w2[i])), returned=True )
        # Avoid the occasional surprise about float32/float64 data types:
        a = a.astype( a1.dtype )
        w = w.astype( a.dtype )

        f1 = mv1.parent # the (open) file corresponding to the FileVariable mv1
        w1id = mv1.id+'_vwgts'
        if not hasattr(mv1,'vwgts'):
            w1axes = mv1.getAxisList()
            w1attributes = {}
            addVariable( f1, w1id, 'd', w1axes, w1attributes )
        f1w = f1[w1id]
        f1w[:] = w1
        f1w[i] = w
        # TypeError: 'CdmsFile' object does not support item assignment    f1[w1id] = w
        mv1.vwgts = w1id
    else:
        # This is what happens most of the time.  It's a simple average (of two compatible numpy
        # arrays), weighted by scalars.  These scalars are the length of time represented by
        # by mv1, mv2.
        # Note that a1,a2 are type TransientVariable which inherits from numpy.ma.MaskedArray
        sw12 = sw1+sw2
        a = ( a1*(sw1/sw12) + a2*(sw2/sw12) )
        try:
            if a.dtype != a1.dtype:
                a = a.astype( a1.dtype )
        except Exception as e:
            # happens if a1 isn't a numpy array, e.g. a float.  Then it's ok to just go on.
            #print "In arithmetic average of",mv1.id,"in two_pt_avg, encountered exception:",e
            #print "a1=",a1,type(a1),a1.dtype if hasattr(a1,'dtype') else None
            pass
    return a

def update_time_avg( redvars, redtime_bnds, redtime_wts, newvars, next_tbounds, dt=None,
                     new_time_weights=None, force_scalar_avg=False ):
    """Updates the time-reduced data for a list of variables.  The reduced-time and averaged
    variables are listed in redvars.  Its weights (for time averaging) are another variable,
    redtime_wts.
    (Each member of redvars should have the same the time axis.  Each member of newvars should have
    the same time axis.  Normally it has an attribute wgts which is the same as redtime_wts.id.)
    The new data is listed in newvars, and this list should correspond to redvars, e.g.
    newvars[i].id==redvars[i].id.  The length of both should be equal and at least one.
    Each variable is an MV.  Normally redvar and redtime_wts will be a FileVariable (required if
    they might change) and newvar a TransientVariable.
    If newvar needs any spatial reductions to match redvar, they should have been performed before
    calling this function.
    next_tbounds is the next time interval, used if newvar is defined on a time beyond redvar's
    present time axis.  If next_tbounds==[], newvar will be ignored on such times.  Normally
    next_tbounds will be set to [] when updating a climatology file which has been initialized.
    The penultimate argument dt is used only in that dt=0 means that we are computing climatologies - hence
    the new data's time axis must be adjusted before averaging the data into redvars.
    The last argument is the time_weights global attribute of the data file, if any; it corresponds
    to newvars.  This is expected to occur iff the data file is a climatology file written by
    an earlier use of this module.
    The optional argument force_scalar_avg argument is for testing and is passed on to two_pt_avg.
    """

    # >>>> TO DO <<<< Ensure that each redvar, redtime_wts, newvar have consistent units
    # >>>> for the variable and for time.  Assert that they have the same shape, axes, etc.

    if redvars is None or len(redvars)==0:  # formerly redvar was initialized here
        raise Exception("update_time_avg requires a reduced variable list")
    nvars = len(redvars)
    # The following two asserts express my assumption that the first index of the variable is time.
    # This is almost universal, but in the future I should generalize this code.  That would make
    # slicing more verbose, e.g. if time were changed from the first index to second then
    # v[j] would become v[:,j,:] (for a 2-D variable v).

    for var in redvars:
        redtime = var.getTime()  # partially-reduced time axis
        if redtime is not None:  # some variables have no time axis
            try:
                assert( var.getDomain()[0][0].isTime() )
            except Exception as e:
                logging.exception("redvars=%s, var=%s, var.getDomain()=%s", redvars, var, var.getDomain())
                raise e
            break
    assert( redtime is not None )
    redtime_len = redtime.shape[0]
    for var in newvars:
        newtime = var.getTime()  # original time axis, from the new variable
        if newtime is not None:
            try:
                assert( var.getDomain()[0][0].isTime() )
            except Exception as e:
                logging.exception("redvars=%s, var=%s, var.getDomain()=%s", redvars, var, var.getDomain())
                raise e
            break
    assert( newtime is not None ) # The input data should have a time axis!
    if dt==0:
        newtime = adjust_time_for_climatology( newtime, redtime )
    newtime_bnds = getClimoBounds(newtime)
    # newtime_wts[j,i] is the weight applied to the data at time newtime[j] in computing
    # an average, reduced time for time newtime[ newtime_rti[i] ], 0<=i<2.
    # If newtime_rti[i]<0, that means the weight is 0.
    maxolaps = 3  # Maximum number of data time intervals which could overlap with a single
    #               reduced-time interval.  We're unlikely to see more than 2.
    newtime_wts = numpy.zeros(( newtime.shape[0], maxolaps ))
    newtime_rti = numpy.zeros(( newtime.shape[0], maxolaps ), numpy.int32) - 1
    for j in range( newtime.shape[0] ):
        # First, extend redtime and redtime_bnds if necessary:
        # This should be moved to a separate function.
        if newtime_bnds[j][1] > redtime_bnds[-1][1]:
            bndmin = max( newtime_bnds[j][0], next_tbounds[0] )
            bndmax = min( newtime_bnds[j][1], next_tbounds[1] )
            weight = bndmax-bndmin
            if weight>0:
                # Extend the time axis to add a new time, time bounds, and weight.  With one more
                # silly step (see below), this will also have the effect of extending redvars along
                # the time axis.
                redtime_bnds[redtime_len] = next_tbounds
                redtime[redtime_len] = 0.5*(
                    redtime_bnds[redtime_len][1] + redtime_bnds[redtime_len][0] )
                redtime_wts[redtime_len] = 0.0
                redtime_len +=1
        for iv in range(nvars):
            # Without this silly step, the data in redvars[iv] won't be expanded to match the
            # newly expanded time axis...
            dummy = redvars[iv].shape
            # This also will do the job, but it looks like a lot of i/o:
            #   redvars[iv].parent.write(redvars[iv])
            # It doesn't help to write redtime_wts or redtime_bnds.  You need to write a variable
            # with the same axes as redvars.
            # This doesn't do the job:  redvars[iv].parent.sync()

        # The weight of time newtime[j] is the part of its bounds which lie within some reduced-
        # time bounds.  We'll also need to remember the indices of the reduced-times for
        # which this is nonzero (there will be few of those, many reduced times total)
        k = -1
        for i,red_bnds in enumerate( redtime_bnds ):
            bndmin = max( newtime_bnds[j][0], red_bnds[0] )
            bndmax = min( newtime_bnds[j][1], red_bnds[1] )
            weight = bndmax-bndmin
            if weight<0:
                continue
            else:
                k += 1
                newtime_wts[j,k] = weight
                newtime_rti[j,k] = i
        #  This much simpler expression works if there is no overlap:
        #newtime_wts[j] = newtime_bnds[j][1] - newtime_bnds[j][0]
        kmax = k
        assert( kmax<maxolaps ) # If k be unlimited, coding is more complicated.
        # This is the first point at which we decide whether the input file covers any times of interest,
        # e.g. for climatology of DJF, here is where we decide whether the file is a D,J,or F file.
        # Here kmax<0 if the file has no interesting time.

    if new_time_weights is not None:
        # The weights from a climatology file make sense only when time is simply structured.
        # Otherwise, we don't know what to do.
        for j,nt in enumerate(newtime):
            for k in range(kmax+1)[1:]:
                assert( newtime_wts[j,k] ==0 )
        newtime_wts = numpy.array([new_time_weights.data])
    for j,nt in enumerate(newtime):
        for k in range(kmax+1):
            i = int( newtime_rti[j][k] )
            # This is the second point at which we decide whether the input file covers any times of interest,
            # e.g. for climatology of DJF, here is where we decide whether the file is a D,J,or F file.
            # Here i<0 if the file has no interesting time.
            if i<0: continue
            for iv in range(nvars):
                redvar = redvars[iv]
                newvar = newvars[iv]
                if redvar.id=='time_bnds' or redvar.id=='time_weights':
                    continue
                if redvar.dtype.kind=='i' and newvar.dtype.kind=='i' or\
                        redvar.dtype.kind=='S' and newvar.dtype.kind=='S' :
                    # integer, any length, or string.  Time average makes no sense.
                    data = newvar
                    if redvar.shape==newvar.shape:
                        redvar.assignValue(newvar)
                    else:
                        redvar[i] = newvar[j]
                    continue
                if 'time' not in redvar.getAxisIds() and redvar.initialized=='yes':
                    # ...This assumes time axis is named 'time', but a full-blown check for a time
                    # axis, e.g. via redvar.getTime(), is surprisingly expensive because of a
                    # udunits inefficiency.
                    # Although there's no time axis (hence index i is irrelevant), values may
                    # differ from one file to the next, so we still have to do a time average.
                    if newtime_wts[j,k]>0:
                        # Maybe the assignValue call will work for other shapes.  To do: try it,
                        # this code will simplify if it works
                        if len(redvar.shape)==0:
                            redvar.assignValue(
                                ( redvar.subSlice()*redtime_wts[i] + newvar*newtime_wts[j,k] ) /\
                                    ( redtime_wts[i] + newtime_wts[j,k] ) )
                        elif len(redvar.shape)==1:
                            redvar[:] =\
                                ( redvar[:]*redtime_wts[i] + newvar[:]*newtime_wts[j,k] ) /\
                                ( redtime_wts[i] + newtime_wts[j,k] )
                        elif len(redvar.shape)==2:
                            redvar[:,:] =\
                                ( redvar[:,:]*redtime_wts[i] + newvar[:,:]*newtime_wts[j,k] ) /\
                                ( redtime_wts[i] + newtime_wts[j,k] )
                        elif len(redvar.shape)==3:
                            redvar[:,:,:] =\
                                ( redvar[:,:,:]*redtime_wts[i] + newvar[:,:,:]*newtime_wts[j,k] ) /\
                                ( redtime_wts[i] + newtime_wts[j,k] )
                        else:
                            #won't work because redvar is a FileVariable
                            logger.warning("Probably miscomputing average of %s", redvar.id)
                            redvar =\
                                ( redvar*redtime_wts[i] + newvar*newtime_wts[j,k] ) /\
                                ( redtime_wts[i] + newtime_wts[j,k] )
                elif 'time' not in redvar.getAxisIds() and redvar.initialized=='no':
                    redvar.assignValue(newvar)
                elif (not hasattr(redvar[i],'mask') and redvar.initialized=='yes') or\
                        (hasattr(redvar[i],'mask') and not redvar[i].mask.all()):
                    # ... i.e., redvar[i] is initialized and is not entirely missing data
                    if newtime_wts[j,k]>0:

                        # Don't miss this line, it's where the average is computed:
                        redvar[i] = two_pt_avg( redvar, newvar, i, newvar[j],
                                                redtime_wts[i], newtime_wts[j,k],
                                                force_scalar_avg=force_scalar_avg )

                else:
                    # For averaging, unitialized is same as value=0 because redtime_wts is
                    # initialized to 0
                    redvar[i] = redvar.dtype.type( newvar[j] )
                    if hasattr( newvar, 'vwgts' ):
                        redvar.vwgts = newvar.vwgts
                        wid = redvar.vwgts
                        frv = redvar.parent
                        if wid not in frv.variables:
                            waxes = redvar.getAxisList()
                            addVariable( frv, wid, 'd', waxes, {} )
                            frvw = frv[wid]
                            frvw[:] = newvar.from_file( wid )[:]
            if redtime_wts.initialized=='yes':
                redtime_wts[i] += newtime_wts[j,k]
            else:      # uninitialized is same as value=0
                redtime_wts[i]  = newtime_wts[j,k]
    for iv in range(nvars):
        redvars[iv].initialized = 'yes'
    redtime_wts.initialized = 'yes'

    return redvars,redtime_wts,redtime

def update_time_avg_from_files( redvars0, redtime_bnds, redtime_wts, filenames,
                                fun_next_tbounds=next_tbounds_copyfrom_data,
                                redfiles=[], dt=None, force_scalar_avg=False, lock=None ):
    """Updates the time-reduced data for a several variables.  The reduced-time and averaged
    variables are the list redvars.  Its weights (for time averaging) are another variable, redtime_wts.
    (Each variable redvar of redvars has the same time axis, and it normally has an attribute wgts
    which is the same as redtime_wts.id.)  redtime_bnds is a list of reduced time bounds.
    Corresponging to it is the weights redtime_wts.
    The new data will be read from the specified files, which should cover a time sequence in order.
    We expect each member of redvar; and redtime_bnds[i] and redtime_wts[i], to be a FileVariable.
    This function doesn't perform any spatial reductions.
    The fun_next_tbounds argument is a function which will compute the next_tbounds argument of
    update_time_avg, i.e. the next time interval (if any) to be appended to the bounds of the time
    axis of redvars.
    redfiles is a list of reduced-time files for output.  They should already be open as 'r+'.
    The optional argument dt is used only in that dt=0 means that we are computing climatologies.
    The optional argument force_scalar_avg argument is for testing and is passed on to two_pt_avg.
    """
    tnewvar = 0
    ttotal = time.time()
    tmin = 1.0e10
    tmax = -1.0e10
    for filen in filenames:
        if lock is not None:  lock.acquire()
        #t1 = time.time()
        f = cdms2.open(filen)
        #t2 = time.time()
        if lock is not None:  lock.release()
        ftime = f.getAxis('time').clone()  # clone saves it in memory - faster
        data_tbounds = getClimoBounds(ftime)
        # Compute tmin, tmax which represent the range of times for the data we computed with.
        # For regular model output files, these should come from the time bounds.
        # For a climatology file, they should come from the climatology attribute of time.
        tmin = min( tmin, data_tbounds[0][0] )
        tmax = max( tmax, data_tbounds[-1][-1] )
        newvard = {}
        redvard = {}
        varids = []
        if 'time_weights' in f.variables:
            new_time_weights=f('time_weights')
        else:
            new_time_weights=None
        for redvar in redvars0:
            try:
                varid = redvar.id
                if varid not in f.variables.keys(): continue #jfp testing
                t1=time.time()
                newvar = f(varid)
                tnewvar += time.time()-t1
                # For testing done for <NAME>...
                #if hasattr(newvar,'dtype') and newvar.dtype=='float32':
                #    newvar = newvar.astype('float64')
                if hasattr(newvar,'id'):  # excludes an ordinary number
                    # Usually newvar is a TransientVariable, hence doesn't have the file as :parent.
                    # We need the file to get associated variables, so:
                    newvar.from_file = f
                testarray = numpy.array([0],newvar.dtype)
                if isinstance( testarray[0], Number):
                    varids.append( varid )
                    if newvar.__class__.__name__.find('Variable')<0:
                        # newvar should be a TransientVariable.  But cdms2 doesn't work very well
                        # for scalar variables.  It reads them as numpy numbers.
                        newvar = cdms2.createVariable(newvar,id=varid)
                    newvard[varid] = newvar
                    redvard[varid] = redvar
                else:
                    logger.info( "skipping %s",redvar.id )
            except Exception as e:
                if varid!='climatology_bnds':  # I know about this one.
                    logging.exception("skipping %s due to exception", redvar.id)
                    logging.exception(e)
                pass
        if len(varids)==0:
            continue
        redvars = [ redvard[varid] for varid in varids ]
        newvars = [ newvard[varid] for varid in varids ]
        if len(redfiles)==0:
            #redvars = redvars0
            tbnds = apply( fun_next_tbounds, ( redtime_bnds, data_tbounds ) )
            update_time_avg( redvars, redtime_bnds, redtime_wts, newvars, tbnds[-1], dt=dt,
                             new_time_weights=new_time_weights, force_scalar_avg=force_scalar_avg )
        else:
            for g in redfiles:
                redtime = g.getAxis('time')
                redtime_wts = g['time_weights']
                redtime_bnds = getClimoBounds(redtime)
                redvars = [ g[varn] for varn in varids ]
                tbnds = apply( fun_next_tbounds, ( redtime_bnds, data_tbounds ) )
                update_time_avg( redvars, redtime_bnds, redtime_wts, newvars, tbnds[-1], dt=dt,
                                 new_time_weights=new_time_weights,
                                 force_scalar_avg=force_scalar_avg )
        if lock is not None:  lock.acquire()
        f.close()
        if lock is not None:  lock.release()
    ttotal = time.time() - ttotal
    return tmin, tmax

def test_time_avg( redfilename, varnames, datafilenames ):
    #dt = None   # if None, the "reduced" time bounds are exactly the same as in the input data
    dt = 365   # if >0, a fixed time step for partially time-reduced data
    #init_red_tbounds = [[]] # if [[]] or  [] or None, a reduced time interval is dt, the time step from
    #                       one interval to the next; thus all time is covered
    f = cdms2.open(datafilenames[0])
    data_time = f.getAxis('time')  # a FileAxis
    init_data_tbounds = data_time.getBounds()[0]
    # N is used to start the intervals off in the right year.  Note that this works only if calendar has a fixed-length year
    N = math.floor(data_time[0]/365.)
    logger.info("N= %s",N)
    logger.info("data_time= %s %s %s",data_time, (data_time[0]-N*365), (data_time[-1]-N*365))
    logger.info("init_data_tbounds= %s %s %s ",init_data_tbounds,[init_data_tbounds[0]-N*365],[init_data_tbounds[1]-N*365])
    init_red_tbounds = numpy.array([[150,240]], dtype=numpy.int32) # example time interval for a single season - not all times
    init_red_tbounds = init_red_tbounds + N*365
    logger.info("init_red_tbounds= %s",init_red_tbounds)

    #                              go into time-reduced data if dt=365
    initialize_redfile_from_datafile( redfilename, varnames, datafilenames[0], dt,
                                      init_red_tbounds )

    g = cdms2.open( redfilename, 'r+' )
    redtime_wts = g['time_weights']
    redtime_bnds = g[ g.getAxis('time').bounds ]
    redvars = [ g[varn] for varn in varnames ]

    if dt is None:
        update_time_avg_from_files( redvars, [redtime_bnds], [redtime_wts], datafilenames )
    else:
        update_time_avg_from_files(
            redvars, [redtime_bnds], [redtime_wts], datafilenames,
            fun_next_tbounds = (lambda rtb,dtb,dt=dt: next_tbounds_prescribed_step(rtb,dtb,dt) ) )

    g.close()

    # For testing, print results...
    g = cdms2.open( redfilename )
    redtime = g.getAxis('time')
    redtime_bnds = g( redtime.bounds )
    redtime_wts = g('time_weights')
    TS = g('TS')
    PS = g('PS')
    logger.info("redtime= %s",redtime)
    logger.info("redtime_bnds= %s",redtime_bnds)
    logger.info("redtime_wts= %s",redtime_wts)
    logger.info("TS= %s",TS)
    logger.info("PS= %s",PS)

if __name__ == '__main__':
    if len( sys.argv )>=2:
        datafilenames = sys.argv[1:]
    else:
        datafilenames = ['b30.009.cam2.h0.0600-01.nc','b30.009.cam2.h0.0600-02.nc']
    redfilename = 'redtest.nc'
    varnames = ['TS', 'PS']
    #test_time_avg( redfilename, varnames, datafilenames )
    test_climos( redfilename, varnames, datafilenames )

#!/usr/bin/env python

# Incremental data reduction.

# Suppose you have model output in files file1,file2,...,fileN.
# All files are the same except that they cover different times (in sequence; together they cover
# all time without overlapping); and of course variable values are different.
# Suppose you have reduced the data, e.g. time from daily to monthly averages, maybe also some space
# reductions and maybe limited to a few variables.  The reduced data is in a file reducf.
# Now we get the next file fileN+1.  Our goal here is to update reducf, without re-opening the
# first N files.  Of course, there is the initial case N=0, where reducf doesn't exist.

# This script, when finished, will also be able to compute climatologies.  Just replace each
# model output file's time units from "...since 1850-00-00" to "...since 0000-00-00" which
# would be the time units of a climatology file.  And shift all its time values by an integer
# number of years, so that all time values fall within the same year (time bounds may slip
# into another year.  There may be an issue here, especially if the time isn't centered
# within its bounds.)

# It will be essential that the time axis be present, with a bounds attribute.
# (However, the output file will also have a climatology attribute on the time axis.)
# Time averaging will be weighted by the size of the time interval.  Thus, for the noleap calendar
# with units in days, January would have a weight of 31 relative to the February weight of 28.
# If units were months, each month would have the same weight.

# On machines, such as rhea.ccs.ornl.gov, which impose a very large random cost to opening
# a file, this script will perform better than one which is based on cdms datasets and cdscan.

# The partial time reduction depends, obviously, on what we use for target time intervals.
# The most general way to specify them is to provide a function, so that will be an option.
# The general case is complicated because the intervals need not be equal.
# But usually users will want "annual", "monthly", or "seasonal climatologies".
# I will provide functions for these cases.
# Note that "annual" intervals are unequal in a Gregorian calenar with leap years, and
# "monthly" intervals are unequal except in a 360-day calendar.
# A keyword + calendar + units will be sufficient if it's one of the likely requests.

# >>>> TO DO:
# Create a suitable time axis when reading climo data without one.

import numpy, cdms2, sys, os, math, logging
from numbers import Number
from pprint import pprint
import time
from metrics.packages.acme_regridder.scripts.acme_regrid import addVariable
import logging

logger = logging.getLogger(__name__)

#import debug

# Silence annoying messages about setting the NetCDF file type.  Also, these three lines will
# ensure that we get NetCDF-3 files.  Expansion of a FileAxis may not work on NetCDF-4 files.
cdms2.setNetcdfShuffleFlag(0)
cdms2.setNetcdfDeflateFlag(0)
cdms2.setNetcdfDeflateLevelFlag(0)

def daybounds( season ):
    """Input is a 3-character season, e.g. JAN, JJA, ANN.  Output is bounds for that season
    in units of "days since 0".
    In most cases the bounds will be like [[lo,hi]].  But for DJF it will be like
    [[lo1,hi1],[lo2,hi2]].
    The noleap (365 day) calendar is required."""
    # This is an incredibly crude way to do it, but works for now.
    # If permanent, I would make these dictionaries global for performance.
    dbddic = { 'JAN': [[0,31]],
               'FEB': [[31,59]],
               'MAR': [[59,90]],
               'APR': [[90,120]],
               'MAY': [[120,151]],
               'JUN': [[151,181]],
               'JUL': [[181,212]],
               'AUG': [[212,243]],
               'SEP': [[243,273]],
               'OCT': [[273,304]],
               'NOV': [[304,334]],
               'DEC': [[334,365]] }
    sesdic = { 'MAM': [[ dbddic['MAR'][0][0], dbddic['MAY'][0][1] ]],
               'JJA': [[ dbddic['JUN'][0][0], dbddic['AUG'][0][1] ]],
               'SON': [[ dbddic['SEP'][0][0], dbddic['NOV'][0][1] ]],
               'DJF': [ dbddic['DEC'][0],
                        [dbddic['JAN'][0][0], dbddic['FEB'][0][1]] ],
               'ANN': [[ dbddic['JAN'][0][0], dbddic['DEC'][0][1] ]] }
    dbddic.update( sesdic )
    return dbddic[season]

def setClimoBounds( time, tfile=None, value=None ):
    """time is a TransientAxis or FileAxis.  tfile is the open file it came from, value a new value
    for its bounds.
    The _bounds_ attribute will be set for future use, as long as the axis is a time axis.
    This will be set from the provided value if any, otherwise will come from files: from the
    time.climatology variable if possible; otherwise the regular getBounds() will be called.
    If both a file and a value are provided, then also the climatology or bounds variables in the
    file will be set from the provided value."""
    if not time.isTime():
        return
    if value is None and tfile is not None:
        # Get _bounds_ from tfile.
        if hasattr( time, 'climatology' ):
            bounds = tfile( time.climatology )[:,:]
            if getattr(bounds,'initialized','yes')=='no':
                time._bounds_ = time.getBounds()
            else:
                time._bounds_ = bounds
        else:
            time._bounds = time.getBounds()
    elif value is None:  # tfile is None
        if hasattr(time,'climatology') and hasattr(time,'parent') and\
                time.parent.__class__.__name__=='CdmsFile':
            setClimoBounds( time, time.parent )
        else:
            time._bounds_ = time.getBounds()
    else:   # value is not None
        # Set _bounds_ from provided value.
        time.setBounds(value)
        time._bounds_ = value  # needed to set _bounds_ for FileVariable
        if tfile is not None:
            # Set bounds variables in tfile.
            if hasattr( time, 'climatology' ):
                f( time.climatology )[:,:] = value
            if hasattr( time, 'bounds' ):
                f( time.bounds )[:,:] = value

    return

def getClimoBounds( time ):
    """time is a TransientAxis or FileAxis, and setClimoBounds() should have already been
    called on it.  This function will return its bounds, from the _bounds_ attribute."""
    if not hasattr( time, '_bounds_' ):
        setClimoBounds( time )
    return time._bounds_

def next_tbounds_copyfrom_data( red_time_bnds, data_time_bnds ):
    """This is one of many ways to provide a time interval for the computing a partial time
    reduction of data.  This one works for no real time reduction.  A typical case is that
    we have a sequence of monthly averages and the 'reduced' time is monthly averages.
    The idea is simply that if the data time bounds are greater than the reduced time bounds,
    use the data time bounds as the next interval in the reduced time sequence.
    red_time_bnds and data_time_bnds may be treated as arrays like a typical time_bnds variable,
    e.g. numpy.array([[1234,1244],[1244,1256],[1256,1267]]).  In practice, red_time_bnds is
    normally a FileVariable with such a value.
    You get an initial value for red_time_bnds by calling this function with red_time_bnds==[].
    """
    if len(data_time_bnds)==0:
        return red_time_bnds
    if len(red_time_bnds)==0:
        return data_time_bnds
    if data_time_bnds[0][0]>=red_time_bnds[-1][1]:
        return numpy.concatenate( ( red_time_bnds, data_time_bnds ), axis=0 )
    else:
        return red_time_bnds

def next_tbounds_prescribed_step( red_time_bnds, data_time_bnds, dt ):
    """The next bounds will be a prescribed step above the existing one.
    This uses data_time_bnds only to use its lowest bound in initialization of red_time_bnds.
    This is intended to be called from another function, e.g. a lambda expression which includes
    a fixed step size."""
    # Normally the numbers in the return value are an integers between 0 and 365,
    # representing a time in days
    if len(red_time_bnds)==0:
        t0 = data_time_bnds[0][0]
        return numpy.array([[t0,t0+dt]])
    t0 = red_time_bnds[-1][0]
    t1 = red_time_bnds[-1][1]
    t1t2 = numpy.array([[t0+dt, t1+dt]],dtype=numpy.int32)
    return numpy.concatenate( ( red_time_bnds, t1t2 ), axis=0 )

def initialize_redfile( filen, axisdict, typedict, attdict, varnames, lock=None, dt=-1 ):
    """Initializes the file containing reduced data.  Input is the filename, axes
    variables representing bounds on axes, and names of variables to be initialized.
    Axes are represented as a dictionary, with items like varname:[timeaxis, lataxis, lonaxis].
    Time should always be present as the first axis.
    Note that cdms2 will automatically write the bounds variable if it writes an axis which has bounds.
    """
    if lock is not None:  lock.acquire()
    #t1 = time.time()
    g = cdms2.open( filen, 'w' )
    #t2 = time.time()
    if lock is not None:  lock.release()
    if 'time_bnds' not in varnames:
        addVariable( g, 'time_bnds', typedict['time_bnds'], axisdict['time_bnds'], attdict['time_bnds'] )
    if 'climatology_bnds' not in varnames and dt==0:
        addVariable( g, 'climatology_bnds', typedict['time_bnds'], axisdict['time_bnds'], attdict['time_bnds'] )
    for varn in varnames:
        addVariable( g, varn, typedict[varn], axisdict[varn], attdict[varn] )
    addVariable( g, 'time_weights', typedict['time'], [axisdict['time']], [] )
    g['time_weights'].initialized = 'no'
    return g;

def initialize_redfile_from_datafile( redfilename, varnames, datafilen, dt=-1, init_red_tbounds=None,
                                      force_double=False, lock=None ):
    """Initializes a file containing the partially-time-reduced average data, given the names
    of the variables to reduce and the name of a data file containing those variables with axes.
    The file is created and returned in an open state.  Closing the file is up to the calling function.
    Also returned is a list of names of variables which have been initialized in the file.
    If provided, dt is a fixed time interval for defining the time axis.
    dt=0 is special in that it creates a climatology file.
    If provided, init_red_tbounds is the first interval for time_bnds.
    For example in the 360-day calendar with time units of days, dt=360 and
    init_red_tbounds=[[150,240]] would set up the partially-time-reduced averages for annual
    trends of values averaged over the JJA season.
    """
    # Note: some of the attributes are set appropriately for climatologies, but maybe not for
    # other time reductions.  In particular, c.f. time long_name and varn cell_methods.
    boundless_axes = set([])
    if lock is not None:  lock.acquire()
    #t1 = time.time()
    f = cdms2.open( datafilen )
    #t2 = time.time()
    if lock is not None:  lock.release()
    axisdict = {}
    typedict = {}
    attdict = {}
    timeaxis = None
    if dt>=0:   # time axis should have fixed intervals, size dt
        if init_red_tbounds is None or init_red_tbounds==[] or init_red_tbounds==[[]]:
            init_red_tbounds=[]
            red_time_bnds = next_tbounds_prescribed_step(
                init_red_tbounds, f.getAxis('time').getBounds(), dt )
        else:
            red_time_bnds = init_red_tbounds
        red_time = [ 0.5*(tb[0]+tb[1]) for tb in red_time_bnds ]
        # tb[?] can be expected to be an integer,  red_time 64-bit float.
        timeaxis = cdms2.createAxis( red_time, red_time_bnds, 'time' )
        # ... this will set the _bounds_ attribute to red_time_bnds by calline timeaxis.setBounds().
        timeaxis.bounds = 'time_bnds'
        timeaxis.climatology = 'climatology_bnds'
        axisdict['time'] = timeaxis
        typedict['time'] = f['time'].typecode()
        attdict['time'] = { 'bounds':'time_bnds', 'long_name':'climatological_time' }
        if dt==0:
            attdict['time']['climatology'] = 'climatology_bnds'
        if force_double and typedict['time']=='f':
            # I think that time is always double, but let's make sure:
            typedict['time']='d'
            new_dtype = numpy.float64
            if '_FillValue' in attdict['time']:
                attdict['time']['_FillValue'] = attdict['time']['_FillValue'].astype( new_dtype )
                if 'missing_value' in attdict['time']:
                    attdict['time']['missing_value'] =\
                        attdict['time']['missing_value'].astype( new_dtype )

    out_varnames = []
    for varn in varnames:
        if f[varn] is None:
            logger.warning("%s was not found in %s", varn, datafilen)
            continue
        try:
            dtnom = f[varn].dtype.name  # works for variables
        except:
            # Probably an axis.  Most of what we're doing here makes no sense for an axis.
            # FWIW, this would get the type name: dtnom = f[varn].dtype('spam').name
            logger.warning("WARNING, averager will ignore axis %s.", varn)
            continue
        if dtnom.find('string')==0:
            # We can't handle string variables.
            logger.warning("Ignoring string variable %s", varn)
            continue
        out_varnames.append(varn)
        dataaxes = f[varn].getAxisList()
        axes = []
        for ax in dataaxes:
            if dt>=0 and ax.isTime():
                axes.append( timeaxis )
            else:
                axes.append( ax )
            if not hasattr( ax,'bounds' ): boundless_axes.add(ax)
        axisdict[varn] = axes
        typedict[varn] = f[varn].typecode()

        # attributes:
        vdict = f[varn].__dict__
        attdict[varn] = {a:vdict[a] for a in vdict if (a=='_FillValue' or a[0]!='_')
                         and a!='parent' and a!='autoApiInfo' and a!='domain' and a!='vwgts'}
        # ...It would be worth reconstructing the domain attribute with the new axes.
        if 'cell_methods' not in attdict[varn]:
            attdict[varn]['cell_methods'] = 'time: mean'
        elif 'time: mean' not in attdict[varn]['cell_methods']:
            attdict[varn]['cell_methods'] += ', time: mean'
        attdict[varn]['initialized'] = 'no'

        if force_double and typedict[varn]=='f':
            typedict[varn]='d'
            new_dtype = numpy.float64
            if '_FillValue' in attdict[varn]:
                attdict[varn]['_FillValue'] = attdict[varn]['_FillValue'].astype( new_dtype )
                if 'missing_value' in attdict[varn]:
                    attdict[varn]['missing_value'] =\
                        attdict[varn]['missing_value'].astype( new_dtype )

    for ax in boundless_axes:
        logger.warning("Axis %s has no bounds", ax.id)
    if timeaxis is not None:
        tbndaxis = cdms2.createAxis( [0,1], None, 'tbnd' )

        axisdict['time_bnds'] = [timeaxis,tbndaxis]
        typedict['time_bnds'] = timeaxis.typecode() # always 'd'
        attdict['time_bnds'] = {'initialized':'no'}

        if dt==0:
            out_varnames.append('climatology_bnds')
            axisdict['climatology_bnds'] = [timeaxis,tbndaxis]
            typedict['climatology_bnds'] = timeaxis.typecode() # always 'd'
            attdict['climatology_bnds'] = {'initialized':'no'}

    g = initialize_redfile( redfilename, axisdict, typedict, attdict, out_varnames,
                            lock=lock, dt=dt )

    if timeaxis is not None:
        g['time_bnds'][:] = red_time_bnds
        g['time_bnds'].initialized = 'yes'
        g['time_weights'][:] = numpy.zeros(len(timeaxis))
        g['time_weights'].initialized = 'yes'
    ftime = f.getAxis('time')
    if hasattr(ftime,'climatology'):
        tb = f( ftime.climatology ).data
    else:
        tb = ftime.getBounds()
    if tb is None:
        tmin = f.getAxis('time')[0]
        tmax = f.getAxis('time')[-1]
    else:
        tmin = tb[0][0]
        tmax = tb[-1][-1]
    if lock is not None:  lock.acquire()
    f.close()
    if lock is not None:  lock.release()
    return g, out_varnames, tmin, tmax

def adjust_time_for_climatology( newtime, redtime ):
    """Changes the time axis newtime for use in a climatology calculation.  That is, the values of
    the time axis are shifted to lie within a single year suitable for averaging with the existing
    reduced (climatology) variables redvars.
    """
    newtime.toRelativeTime( redtime.units, redtime.getCalendar() )
    assert( newtime.calendar=='noleap' )
    # btw, newtime.getCalendar()==4113 means newtime.calendar=="noleap"

    # We need to adjust the time to be in the year 0, the year we use for climatology calculations.
    # It would be natural to compute the adjustment from newtime[0].  But sometimes it's equal to
    # one of its bounds, on the border between two years!  So it's better to use the midpoint
    # between two bounds.
    # Usually we can assume that the time bounds are on monthly borders.  This should be checked, but isn't yet.<<<<
    timebnds = getClimoBounds(newtime)
    midtime = 0.5*(timebnds[0][0]+timebnds[0][1])
    N = math.floor(midtime/365.)   # >>>> assumes noleap calendar, day time units!
    newtime[:] -= N*365               # >>>> assumes noleap calendar, day time units!
    setClimoBounds( newtime, value=(getClimoBounds(newtime) - N*365) ) # >>>> assumes noleap calendar, day time units!
    return newtime

def two_pt_avg( mv1, mv2, i, a2, sw1, sw2, aw2=None, force_scalar_avg=False ):
    """Input:
    FileVariables mv1,mv2;  i, (time) index to subset mv1 as a1=mv1[i], a
    corresponding array a2 from mv2, also formed by fixing a point on the first, time, axis;
    scalar weights sw1,sw2.
    Output:   weighted average a of a1 and a2.
    The weights will be scalar weights unless an array weight must be used because either:
    an array weight is already an attribute of mv1 or mv2; or
    v1 and v2 do not have the same mask.
    If an array weight is used, it will become an attribute of the first variable mv1,
    and take its shape.  Normally mv1 is going to be where averages get accumulated,
    so if there is more than one time, it will appear in mv1 but not the input data mv2.
    The optional argument aw2 is meant to support array weights in the case where we are
    averaging different parts of the same variable, i.e. mv1=mv2 but a1==mv1[i]!=a2.
    If aw2 be supplied, it will be used as the array weight for mv2,a2.
    The optional argument force_scalar_avg argument is for testing.  If set to True,
    all computations involving array weights will be bypassed and simple scalar
    weights used instead.
    """
    a1 = mv1[i]
    if not force_scalar_avg:
        w1 = sw1
        if aw2 is not None:
            w2 = aw2
        else:
            w2 = sw2

        # If array weights already exist for mv1 _or_ mv2, use them, and create array weights
        # for whatever doesn't have them.  If mv1,mv2 have different masks, then we need array
        # weights now and henceforth; create them.  The shape attribute is a convenient way
        # to detect whether w1 or w2 is scalar or array.
        if hasattr( mv1, 'vwgts' ):
            f1 = mv1.parent # the (open) file corresponding to the FileVariable mv1
            w1 = f1(mv1.vwgts)
        if hasattr( mv2, 'vwgts' ) and aw2 is None:
            f2 = mv2.parent # the (open) file corresponding to the FileVariable mv2
            if f2 is None:
                f2 = mv2.from_file # the (open) file corresponding to the TransientVariable mv2
            w2 = f2(mv2.vwgts)
        if (not hasattr(w1,'shape') or len(w1.shape)==0) and hasattr(w2,'shape') and len(w2.shape)>0:
            w1 = numpy.full( mv1.shape, -1 )
            w1[i] = sw1
            w1 = numpy.ma.masked_less(w1,0)
        if (not hasattr(w2,'shape') or len(w2.shape)==0) and hasattr(w1,'shape') and len(w1.shape)>0:
            w2 = numpy.full( mv1.shape, -1 )   # assumes mv1 time axis is >= mv2 time axis
            w2[i] = sw2
            w2 = numpy.ma.masked_less(w2,0)
        if (not hasattr(w1,'shape') or len(w1.shape)==0) and\
                (not hasattr(w2,'shape') or len(w2.shape)==0):
            mask1 = False
            mask2 = False
            if hasattr(mv1,'_mask'):
                mask1 = mv1.mask
            elif getattr(mv1,'missing_value'):  # if no missing_value, then no mask
                valu = mv1.getValue()
                if hasattr( valu, '_mask' ):
                    mask1 = valu._mask
            if hasattr(mv2,'_mask'):
                mask2 = mv2.mask
            elif getattr(mv2,'missing_value'):  # if no missing_value, then no mask
                valu = mv2.getValue()
                if hasattr( valu, '_mask' ):
                    mask2 = valu._mask
            if not numpy.all(mask1==mask2):
                # Note that this test requires reading all the data.  That has to be done anyway to compute
                # the average.  Let's hope that the system caches well enough so that there won't be any
                # file access cost to this.
                #jfp was w1 = numpy.full( a1.shape, w1 )
                #jfp was w2 = numpy.full( a2.shape, w2 )
                w1 = numpy.full( mv1.shape, -1 )
                w1[i] = sw1
                w1 = numpy.ma.masked_less(w1,0)
                w2 = numpy.full( mv1.shape, -1 )
                w2[i] = sw2
                w2 = numpy.ma.masked_less(w2,0)

    if not force_scalar_avg and\
            hasattr(w1,'shape') and len(w1.shape)>0 and hasattr(w2,'shape') and len(w2.shape)>0:
        if w1[i].mask.all():   # if w1[i] is all missing values:
            w1[i] = sw1
        if w2[i].mask.all():   # if w2[i] is all missing values:
            w2[i] = sw2
        # Here's what I think the numpy.ma averager does about weights and missing values:
        # The output weight w(i)=sw1+sw2 if there be no missing value for mv1(i) and mv2(i) (maybe
        # also if both be missing, because the weight doesn't matter then).
        # But if i be missing for mv1, drop sw1, thus w(i)=sw2.  If i be missing for mv2, drop sw2.
        a,w = numpy.ma.average( numpy.ma.array((a1,a2)), axis=0,
                                weights=numpy.ma.array((w1[i],w2[i])), returned=True )
        # Avoid the occasional surprise about float32/float64 data types:
        a = a.astype( a1.dtype )
        w = w.astype( a.dtype )

        f1 = mv1.parent # the (open) file corresponding to the FileVariable mv1
        w1id = mv1.id+'_vwgts'
        if not hasattr(mv1,'vwgts'):
            w1axes = mv1.getAxisList()
            w1attributes = {}
            addVariable( f1, w1id, 'd', w1axes, w1attributes )
        f1w = f1[w1id]
        f1w[:] = w1
        f1w[i] = w
        # TypeError: 'CdmsFile' object does not support item assignment    f1[w1id] = w
        mv1.vwgts = w1id
    else:
        # This is what happens most of the time.  It's a simple average (of two compatible numpy
        # arrays), weighted by scalars.  These scalars are the length of time represented by
        # by mv1, mv2.
        # Note that a1,a2 are type TransientVariable which inherits from numpy.ma.MaskedArray
        sw12 = sw1+sw2
        a = ( a1*(sw1/sw12) + a2*(sw2/sw12) )
        try:
            if a.dtype != a1.dtype:
                a = a.astype( a1.dtype )
        except Exception as e:
            # happens if a1 isn't a numpy array, e.g. a float.  Then it's ok to just go on.
            #print "In arithmetic average of",mv1.id,"in two_pt_avg, encountered exception:",e
            #print "a1=",a1,type(a1),a1.dtype if hasattr(a1,'dtype') else None
            pass
    return a

def update_time_avg( redvars, redtime_bnds, redtime_wts, newvars, next_tbounds, dt=None,
                     new_time_weights=None, force_scalar_avg=False ):
    """Updates the time-reduced data for a list of variables.  The reduced-time and averaged
    variables are listed in redvars.  Its weights (for time averaging) are another variable,
    redtime_wts.
    (Each member of redvars should have the same the time axis.  Each member of newvars should have
    the same time axis.  Normally it has an attribute wgts which is the same as redtime_wts.id.)
    The new data is listed in newvars, and this list should correspond to redvars, e.g.
    newvars[i].id==redvars[i].id.  The length of both should be equal and at least one.
    Each variable is an MV.  Normally redvar and redtime_wts will be a FileVariable (required if
    they might change) and newvar a TransientVariable.
    If newvar needs any spatial reductions to match redvar, they should have been performed before
    calling this function.
    next_tbounds is the next time interval, used if newvar is defined on a time beyond redvar's
    present time axis.  If next_tbounds==[], newvar will be ignored on such times.  Normally
    next_tbounds will be set to [] when updating a climatology file which has been initialized.
    The penultimate argument dt is used only in that dt=0 means that we are computing climatologies - hence
    the new data's time axis must be adjusted before averaging the data into redvars.
    The last argument is the time_weights global attribute of the data file, if any; it corresponds
    to newvars.  This is expected to occur iff the data file is a climatology file written by
    an earlier use of this module.
    The optional argument force_scalar_avg argument is for testing and is passed on to two_pt_avg.
    """

    # >>>> TO DO <<<< Ensure that each redvar, redtime_wts, newvar have consistent units
    # >>>> for the variable and for time.  Assert that they have the same shape, axes, etc.

    if redvars is None or len(redvars)==0:  # formerly redvar was initialized here
        raise Exception("update_time_avg requires a reduced variable list")
    nvars = len(redvars)
    # The following two asserts express my assumption that the first index of the variable is time.
    # This is almost universal, but in the future I should generalize this code.  That would make
    # slicing more verbose, e.g. if time were changed from the first index to second then
    # v[j] would become v[:,j,:] (for a 2-D variable v).

    for var in redvars:
        redtime = var.getTime()  # partially-reduced time axis
        if redtime is not None:  # some variables have no time axis
            try:
                assert( var.getDomain()[0][0].isTime() )
            except Exception as e:
                logging.exception("redvars=%s, var=%s, var.getDomain()=%s", redvars, var, var.getDomain())
                raise e
            break
    assert( redtime is not None )
    redtime_len = redtime.shape[0]
    for var in newvars:
        newtime = var.getTime()  # original time axis, from the new variable
        if newtime is not None:
            try:
                assert( var.getDomain()[0][0].isTime() )
            except Exception as e:
                logging.exception("redvars=%s, var=%s, var.getDomain()=%s", redvars, var, var.getDomain())
                raise e
            break
    assert( newtime is not None ) # The input data should have a time axis!
    if dt==0:
        newtime = adjust_time_for_climatology( newtime, redtime )
    newtime_bnds = getClimoBounds(newtime)
    # newtime_wts[j,i] is the weight applied to the data at time newtime[j] in computing
    # an average, reduced time for time newtime[ newtime_rti[i] ], 0<=i<2.
    # If newtime_rti[i]<0, that means the weight is 0.
    maxolaps = 3  # Maximum number of data time intervals which could overlap with a single
    #               reduced-time interval.  We're unlikely to see more than 2.
    newtime_wts = numpy.zeros(( newtime.shape[0], maxolaps ))
    newtime_rti = numpy.zeros(( newtime.shape[0], maxolaps ), numpy.int32) - 1
    for j in range( newtime.shape[0] ):
        # First, extend redtime and redtime_bnds if necessary:
        # This should be moved to a separate function.
        if newtime_bnds[j][1] > redtime_bnds[-1][1]:
            bndmin = max( newtime_bnds[j][0], next_tbounds[0] )
            bndmax = min( newtime_bnds[j][1], next_tbounds[1] )
            weight = bndmax-bndmin
            if weight>0:
                # Extend the time axis to add a new time, time bounds, and weight.  With one more
                # silly step (see below), this will also have the effect of extending redvars along
                # the time axis.
                redtime_bnds[redtime_len] = next_tbounds
                redtime[redtime_len] = 0.5*(
                    redtime_bnds[redtime_len][1] + redtime_bnds[redtime_len][0] )
                redtime_wts[redtime_len] = 0.0
                redtime_len +=1
        for iv in range(nvars):
            # Without this silly step, the data in redvars[iv] won't be expanded to match the
            # newly expanded time axis...
            dummy = redvars[iv].shape
            # This also will do the job, but it looks like a lot of i/o:
            #   redvars[iv].parent.write(redvars[iv])
            # It doesn't help to write redtime_wts or redtime_bnds.  You need to write a variable
            # with the same axes as redvars.
            # This doesn't do the job:  redvars[iv].parent.sync()

        # The weight of time newtime[j] is the part of its bounds which lie within some reduced-
        # time bounds.  We'll also need to remember the indices of the reduced-times for
        # which this is nonzero (there will be few of those, many reduced times total)
        k = -1
        for i,red_bnds in enumerate( redtime_bnds ):
            bndmin = max( newtime_bnds[j][0], red_bnds[0] )
            bndmax = min( newtime_bnds[j][1], red_bnds[1] )
            weight = bndmax-bndmin
            if weight<0:
                continue
            else:
                k += 1
                newtime_wts[j,k] = weight
                newtime_rti[j,k] = i
        #  This much simpler expression works if there is no overlap:
        #newtime_wts[j] = newtime_bnds[j][1] - newtime_bnds[j][0]
        kmax = k
        assert( kmax<maxolaps ) # If k be unlimited, coding is more complicated.
        # This is the first point at which we decide whether the input file covers any times of interest,
        # e.g. for climatology of DJF, here is where we decide whether the file is a D,J,or F file.
        # Here kmax<0 if the file has no interesting time.

    if new_time_weights is not None:
        # The weights from a climatology file make sense only when time is simply structured.
        # Otherwise, we don't know what to do.
        for j,nt in enumerate(newtime):
            for k in range(kmax+1)[1:]:
                assert( newtime_wts[j,k] ==0 )
        newtime_wts = numpy.array([new_time_weights.data])
    for j,nt in enumerate(newtime):
        for k in range(kmax+1):
            i = int( newtime_rti[j][k] )
            # This is the second point at which we decide whether the input file covers any times of interest,
            # e.g. for climatology of DJF, here is where we decide whether the file is a D,J,or F file.
            # Here i<0 if the file has no interesting time.
            if i<0: continue
            for iv in range(nvars):
                redvar = redvars[iv]
                newvar = newvars[iv]
                if redvar.id=='time_bnds' or redvar.id=='time_weights':
                    continue
                if redvar.dtype.kind=='i' and newvar.dtype.kind=='i' or\
                        redvar.dtype.kind=='S' and newvar.dtype.kind=='S' :
                    # integer, any length, or string.  Time average makes no sense.
                    data = newvar
                    if redvar.shape==newvar.shape:
                        redvar.assignValue(newvar)
                    else:
                        redvar[i] = newvar[j]
                    continue
                if 'time' not in redvar.getAxisIds() and redvar.initialized=='yes':
                    # ...This assumes time axis is named 'time', but a full-blown check for a time
                    # axis, e.g. via redvar.getTime(), is surprisingly expensive because of a
                    # udunits inefficiency.
                    # Although there's no time axis (hence index i is irrelevant), values may
                    # differ from one file to the next, so we still have to do a time average.
                    if newtime_wts[j,k]>0:
                        # Maybe the assignValue call will work for other shapes.  To do: try it,
                        # this code will simplify if it works
                        if len(redvar.shape)==0:
                            redvar.assignValue(
                                ( redvar.subSlice()*redtime_wts[i] + newvar*newtime_wts[j,k] ) /\
                                    ( redtime_wts[i] + newtime_wts[j,k] ) )
                        elif len(redvar.shape)==1:
                            redvar[:] =\
                                ( redvar[:]*redtime_wts[i] + newvar[:]*newtime_wts[j,k] ) /\
                                ( redtime_wts[i] + newtime_wts[j,k] )
                        elif len(redvar.shape)==2:
                            redvar[:,:] =\
                                ( redvar[:,:]*redtime_wts[i] + newvar[:,:]*newtime_wts[j,k] ) /\
                                ( redtime_wts[i] + newtime_wts[j,k] )
                        elif len(redvar.shape)==3:
                            redvar[:,:,:] =\
                                ( redvar[:,:,:]*redtime_wts[i] + newvar[:,:,:]*newtime_wts[j,k] ) /\
                                ( redtime_wts[i] + newtime_wts[j,k] )
                        else:
                            #won't work because redvar is a FileVariable
                            logger.warning("Probably miscomputing average of %s", redvar.id)
                            redvar =\
                                ( redvar*redtime_wts[i] + newvar*newtime_wts[j,k] ) /\
                                ( redtime_wts[i] + newtime_wts[j,k] )
                elif 'time' not in redvar.getAxisIds() and redvar.initialized=='no':
                    redvar.assignValue(newvar)
                elif (not hasattr(redvar[i],'mask') and redvar.initialized=='yes') or\
                        (hasattr(redvar[i],'mask') and not redvar[i].mask.all()):
                    # ... i.e., redvar[i] is initialized and is not entirely missing data
                    if newtime_wts[j,k]>0:

                        # Don't miss this line, it's where the average is computed:
                        redvar[i] = two_pt_avg( redvar, newvar, i, newvar[j],
                                                redtime_wts[i], newtime_wts[j,k],
                                                force_scalar_avg=force_scalar_avg )

                else:
                    # For averaging, unitialized is same as value=0 because redtime_wts is
                    # initialized to 0
                    redvar[i] = redvar.dtype.type( newvar[j] )
                    if hasattr( newvar, 'vwgts' ):
                        redvar.vwgts = newvar.vwgts
                        wid = redvar.vwgts
                        frv = redvar.parent
                        if wid not in frv.variables:
                            waxes = redvar.getAxisList()
                            addVariable( frv, wid, 'd', waxes, {} )
                            frvw = frv[wid]
                            frvw[:] = newvar.from_file( wid )[:]
            if redtime_wts.initialized=='yes':
                redtime_wts[i] += newtime_wts[j,k]
            else:      # uninitialized is same as value=0
                redtime_wts[i]  = newtime_wts[j,k]
    for iv in range(nvars):
        redvars[iv].initialized = 'yes'
    redtime_wts.initialized = 'yes'

    return redvars,redtime_wts,redtime

def update_time_avg_from_files( redvars0, redtime_bnds, redtime_wts, filenames,
                                fun_next_tbounds=next_tbounds_copyfrom_data,
                                redfiles=[], dt=None, force_scalar_avg=False, lock=None ):
    """Updates the time-reduced data for a several variables.  The reduced-time and averaged
    variables are the list redvars.  Its weights (for time averaging) are another variable, redtime_wts.
    (Each variable redvar of redvars has the same time axis, and it normally has an attribute wgts
    which is the same as redtime_wts.id.)  redtime_bnds is a list of reduced time bounds.
    Corresponging to it is the weights redtime_wts.
    The new data will be read from the specified files, which should cover a time sequence in order.
    We expect each member of redvar; and redtime_bnds[i] and redtime_wts[i], to be a FileVariable.
    This function doesn't perform any spatial reductions.
    The fun_next_tbounds argument is a function which will compute the next_tbounds argument of
    update_time_avg, i.e. the next time interval (if any) to be appended to the bounds of the time
    axis of redvars.
    redfiles is a list of reduced-time files for output.  They should already be open as 'r+'.
    The optional argument dt is used only in that dt=0 means that we are computing climatologies.
    The optional argument force_scalar_avg argument is for testing and is passed on to two_pt_avg.
    """
    tnewvar = 0
    ttotal = time.time()
    tmin = 1.0e10
    tmax = -1.0e10
    for filen in filenames:
        if lock is not None:  lock.acquire()
        #t1 = time.time()
        f = cdms2.open(filen)
        #t2 = time.time()
        if lock is not None:  lock.release()
        ftime = f.getAxis('time').clone()  # clone saves it in memory - faster
        data_tbounds = getClimoBounds(ftime)
        # Compute tmin, tmax which represent the range of times for the data we computed with.
        # For regular model output files, these should come from the time bounds.
        # For a climatology file, they should come from the climatology attribute of time.
        tmin = min( tmin, data_tbounds[0][0] )
        tmax = max( tmax, data_tbounds[-1][-1] )
        newvard = {}
        redvard = {}
        varids = []
        if 'time_weights' in f.variables:
            new_time_weights=f('time_weights')
        else:
            new_time_weights=None
        for redvar in redvars0:
            try:
                varid = redvar.id
                if varid not in f.variables.keys(): continue #jfp testing
                t1=time.time()
                newvar = f(varid)
                tnewvar += time.time()-t1
                # For testing done for <NAME>...
                #if hasattr(newvar,'dtype') and newvar.dtype=='float32':
                #    newvar = newvar.astype('float64')
                if hasattr(newvar,'id'):  # excludes an ordinary number
                    # Usually newvar is a TransientVariable, hence doesn't have the file as :parent.
                    # We need the file to get associated variables, so:
                    newvar.from_file = f
                testarray = numpy.array([0],newvar.dtype)
                if isinstance( testarray[0], Number):
                    varids.append( varid )
                    if newvar.__class__.__name__.find('Variable')<0:
                        # newvar should be a TransientVariable.  But cdms2 doesn't work very well
                        # for scalar variables.  It reads them as numpy numbers.
                        newvar = cdms2.createVariable(newvar,id=varid)
                    newvard[varid] = newvar
                    redvard[varid] = redvar
                else:
                    logger.info( "skipping %s",redvar.id )
            except Exception as e:
                if varid!='climatology_bnds':  # I know about this one.
                    logging.exception("skipping %s due to exception", redvar.id)
                    logging.exception(e)
                pass
        if len(varids)==0:
            continue
        redvars = [ redvard[varid] for varid in varids ]
        newvars = [ newvard[varid] for varid in varids ]
        if len(redfiles)==0:
            #redvars = redvars0
            tbnds = apply( fun_next_tbounds, ( redtime_bnds, data_tbounds ) )
            update_time_avg( redvars, redtime_bnds, redtime_wts, newvars, tbnds[-1], dt=dt,
                             new_time_weights=new_time_weights, force_scalar_avg=force_scalar_avg )
        else:
            for g in redfiles:
                redtime = g.getAxis('time')
                redtime_wts = g['time_weights']
                redtime_bnds = getClimoBounds(redtime)
                redvars = [ g[varn] for varn in varids ]
                tbnds = apply( fun_next_tbounds, ( redtime_bnds, data_tbounds ) )
                update_time_avg( redvars, redtime_bnds, redtime_wts, newvars, tbnds[-1], dt=dt,
                                 new_time_weights=new_time_weights,
                                 force_scalar_avg=force_scalar_avg )
        if lock is not None:  lock.acquire()
        f.close()
        if lock is not None:  lock.release()
    ttotal = time.time() - ttotal
    return tmin, tmax

def test_time_avg( redfilename, varnames, datafilenames ):
    #dt = None   # if None, the "reduced" time bounds are exactly the same as in the input data
    dt = 365   # if >0, a fixed time step for partially time-reduced data
    #init_red_tbounds = [[]] # if [[]] or  [] or None, a reduced time interval is dt, the time step from
    #                       one interval to the next; thus all time is covered
    f = cdms2.open(datafilenames[0])
    data_time = f.getAxis('time')  # a FileAxis
    init_data_tbounds = data_time.getBounds()[0]
    # N is used to start the intervals off in the right year.  Note that this works only if calendar has a fixed-length year
    N = math.floor(data_time[0]/365.)
    logger.info("N= %s",N)
    logger.info("data_time= %s %s %s",data_time, (data_time[0]-N*365), (data_time[-1]-N*365))
    logger.info("init_data_tbounds= %s %s %s ",init_data_tbounds,[init_data_tbounds[0]-N*365],[init_data_tbounds[1]-N*365])
    init_red_tbounds = numpy.array([[150,240]], dtype=numpy.int32) # example time interval for a single season - not all times
    init_red_tbounds = init_red_tbounds + N*365
    logger.info("init_red_tbounds= %s",init_red_tbounds)

    #                              go into time-reduced data if dt=365
    initialize_redfile_from_datafile( redfilename, varnames, datafilenames[0], dt,
                                      init_red_tbounds )

    g = cdms2.open( redfilename, 'r+' )
    redtime_wts = g['time_weights']
    redtime_bnds = g[ g.getAxis('time').bounds ]
    redvars = [ g[varn] for varn in varnames ]

    if dt is None:
        update_time_avg_from_files( redvars, [redtime_bnds], [redtime_wts], datafilenames )
    else:
        update_time_avg_from_files(
            redvars, [redtime_bnds], [redtime_wts], datafilenames,
            fun_next_tbounds = (lambda rtb,dtb,dt=dt: next_tbounds_prescribed_step(rtb,dtb,dt) ) )

    g.close()

    # For testing, print results...
    g = cdms2.open( redfilename )
    redtime = g.getAxis('time')
    redtime_bnds = g( redtime.bounds )
    redtime_wts = g('time_weights')
    TS = g('TS')
    PS = g('PS')
    logger.info("redtime= %s",redtime)
    logger.info("redtime_bnds= %s",redtime_bnds)
    logger.info("redtime_wts= %s",redtime_wts)
    logger.info("TS= %s",TS)
    logger.info("PS= %s",PS)

if __name__ == '__main__':
    if len( sys.argv )>=2:
        datafilenames = sys.argv[1:]
    else:
        datafilenames = ['b30.009.cam2.h0.0600-01.nc','b30.009.cam2.h0.0600-02.nc']
    redfilename = 'redtest.nc'
    varnames = ['TS', 'PS']
    #test_time_avg( redfilename, varnames, datafilenames )
    test_climos( redfilename, varnames, datafilenames )
#!/usr/bin/env python # Incremental data reduction. # Suppose you have model output in files file1,file2,...,fileN. # All files are the same except that they cover different times (in sequence; together they cover # all time without overlapping); and of course variable values are different. # Suppose you have reduced the data, e.g. time from daily to monthly averages, maybe also some space # reductions and maybe limited to a few variables.  The reduced data is in a file reducf. # Now we get the next file fileN+1.  Our goal here is to update reducf, without re-opening the # first N files.  Of course, there is the initial case N=0, where reducf doesn't exist. # This script, when finished, will also be able to compute climatologies.  Just replace each # model output file's time units from "...since 1850-00-00" to "...since 0000-00-00" which # would be the time units of a climatology file.  And shift all its time values by an integer # number of years, so that all time values fall within the same year (time bounds may slip # into another year.  There may be an issue here, especially if the time isn't centered # within its bounds.) # It will be essential that the time axis be present, with a bounds attribute. # (However, the output file will also have a climatology attribute on the time axis.) # Time averaging will be weighted by the size of the time interval.  Thus, for the noleap calendar # with units in days, January would have a weight of 31 relative to the February weight of 28. # If units were months, each month would have the same weight. # On machines, such as rhea.ccs.ornl.gov, which impose a very large random cost to opening # a file, this script will perform better than one which is based on cdms datasets and cdscan. # The partial time reduction depends, obviously, on what we use for target time intervals. # The most general way to specify them is to provide a function, so that will be an option. # The general case is complicated because the intervals need not be equal. # But usually users will want "annual", "monthly", or "seasonal climatologies". # I will provide functions for these cases. # Note that "annual" intervals are unequal in a Gregorian calenar with leap years, and # "monthly" intervals are unequal except in a 360-day calendar. # A keyword + calendar + units will be sufficient if it's one of the likely requests. # >>>> TO DO: # Create a suitable time axis when reading climo data without one. #import debug # Silence annoying messages about setting the NetCDF file type.  Also, these three lines will # ensure that we get NetCDF-3 files.  Expansion of a FileAxis may not work on NetCDF-4 files. Input is a 3-character season, e.g. JAN, JJA, ANN.  Output is bounds for that season     in units of "days since 0".     In most cases the bounds will be like [[lo,hi]].  But for DJF it will be like     [[lo1,hi1],[lo2,hi2]].     The noleap (365 day) calendar is required. # This is an incredibly crude way to do it, but works for now. # If permanent, I would make these dictionaries global for performance. time is a TransientAxis or FileAxis.  tfile is the open file it came from, value a new value     for its bounds.     The _bounds_ attribute will be set for future use, as long as the axis is a time axis.     This will be set from the provided value if any, otherwise will come from files: from the     time.climatology variable if possible; otherwise the regular getBounds() will be called.     If both a file and a value are provided, then also the climatology or bounds variables in the     file will be set from the provided value. # Get _bounds_ from tfile. # tfile is None # value is not None # Set _bounds_ from provided value. # needed to set _bounds_ for FileVariable # Set bounds variables in tfile. time is a TransientAxis or FileAxis, and setClimoBounds() should have already been     called on it.  This function will return its bounds, from the _bounds_ attribute. This is one of many ways to provide a time interval for the computing a partial time     reduction of data.  This one works for no real time reduction.  A typical case is that     we have a sequence of monthly averages and the 'reduced' time is monthly averages.     The idea is simply that if the data time bounds are greater than the reduced time bounds,     use the data time bounds as the next interval in the reduced time sequence.     red_time_bnds and data_time_bnds may be treated as arrays like a typical time_bnds variable,     e.g. numpy.array([[1234,1244],[1244,1256],[1256,1267]]).  In practice, red_time_bnds is     normally a FileVariable with such a value.     You get an initial value for red_time_bnds by calling this function with red_time_bnds==[]. The next bounds will be a prescribed step above the existing one.     This uses data_time_bnds only to use its lowest bound in initialization of red_time_bnds.     This is intended to be called from another function, e.g. a lambda expression which includes     a fixed step size. # Normally the numbers in the return value are an integers between 0 and 365, # representing a time in days Initializes the file containing reduced data.  Input is the filename, axes     variables representing bounds on axes, and names of variables to be initialized.     Axes are represented as a dictionary, with items like varname:[timeaxis, lataxis, lonaxis].     Time should always be present as the first axis.     Note that cdms2 will automatically write the bounds variable if it writes an axis which has bounds. #t1 = time.time() #t2 = time.time() Initializes a file containing the partially-time-reduced average data, given the names     of the variables to reduce and the name of a data file containing those variables with axes.     The file is created and returned in an open state.  Closing the file is up to the calling function.     Also returned is a list of names of variables which have been initialized in the file.     If provided, dt is a fixed time interval for defining the time axis.     dt=0 is special in that it creates a climatology file.     If provided, init_red_tbounds is the first interval for time_bnds.     For example in the 360-day calendar with time units of days, dt=360 and     init_red_tbounds=[[150,240]] would set up the partially-time-reduced averages for annual     trends of values averaged over the JJA season. # Note: some of the attributes are set appropriately for climatologies, but maybe not for # other time reductions.  In particular, c.f. time long_name and varn cell_methods. #t1 = time.time() #t2 = time.time() # time axis should have fixed intervals, size dt # tb[?] can be expected to be an integer,  red_time 64-bit float. # ... this will set the _bounds_ attribute to red_time_bnds by calline timeaxis.setBounds(). # I think that time is always double, but let's make sure: # works for variables # Probably an axis.  Most of what we're doing here makes no sense for an axis. # FWIW, this would get the type name: dtnom = f[varn].dtype('spam').name # We can't handle string variables. # attributes: # ...It would be worth reconstructing the domain attribute with the new axes. # always 'd' # always 'd' Changes the time axis newtime for use in a climatology calculation.  That is, the values of     the time axis are shifted to lie within a single year suitable for averaging with the existing     reduced (climatology) variables redvars. # btw, newtime.getCalendar()==4113 means newtime.calendar=="noleap" # We need to adjust the time to be in the year 0, the year we use for climatology calculations. # It would be natural to compute the adjustment from newtime[0].  But sometimes it's equal to # one of its bounds, on the border between two years!  So it's better to use the midpoint # between two bounds. # Usually we can assume that the time bounds are on monthly borders.  This should be checked, but isn't yet.<<<< # >>>> assumes noleap calendar, day time units! # >>>> assumes noleap calendar, day time units! # >>>> assumes noleap calendar, day time units! Input:     FileVariables mv1,mv2;  i, (time) index to subset mv1 as a1=mv1[i], a     corresponding array a2 from mv2, also formed by fixing a point on the first, time, axis;     scalar weights sw1,sw2.     Output:   weighted average a of a1 and a2.     The weights will be scalar weights unless an array weight must be used because either:     an array weight is already an attribute of mv1 or mv2; or     v1 and v2 do not have the same mask.     If an array weight is used, it will become an attribute of the first variable mv1,     and take its shape.  Normally mv1 is going to be where averages get accumulated,     so if there is more than one time, it will appear in mv1 but not the input data mv2.     The optional argument aw2 is meant to support array weights in the case where we are     averaging different parts of the same variable, i.e. mv1=mv2 but a1==mv1[i]!=a2.     If aw2 be supplied, it will be used as the array weight for mv2,a2.     The optional argument force_scalar_avg argument is for testing.  If set to True,     all computations involving array weights will be bypassed and simple scalar     weights used instead. # If array weights already exist for mv1 _or_ mv2, use them, and create array weights # for whatever doesn't have them.  If mv1,mv2 have different masks, then we need array # weights now and henceforth; create them.  The shape attribute is a convenient way # to detect whether w1 or w2 is scalar or array. # the (open) file corresponding to the FileVariable mv1 # the (open) file corresponding to the FileVariable mv2 # the (open) file corresponding to the TransientVariable mv2 # assumes mv1 time axis is >= mv2 time axis # if no missing_value, then no mask # if no missing_value, then no mask # Note that this test requires reading all the data.  That has to be done anyway to compute # the average.  Let's hope that the system caches well enough so that there won't be any # file access cost to this. #jfp was w1 = numpy.full( a1.shape, w1 ) #jfp was w2 = numpy.full( a2.shape, w2 ) # if w1[i] is all missing values: # if w2[i] is all missing values: # Here's what I think the numpy.ma averager does about weights and missing values: # The output weight w(i)=sw1+sw2 if there be no missing value for mv1(i) and mv2(i) (maybe # also if both be missing, because the weight doesn't matter then). # But if i be missing for mv1, drop sw1, thus w(i)=sw2.  If i be missing for mv2, drop sw2. # Avoid the occasional surprise about float32/float64 data types: # the (open) file corresponding to the FileVariable mv1 # TypeError: 'CdmsFile' object does not support item assignment    f1[w1id] = w # This is what happens most of the time.  It's a simple average (of two compatible numpy # arrays), weighted by scalars.  These scalars are the length of time represented by # by mv1, mv2. # Note that a1,a2 are type TransientVariable which inherits from numpy.ma.MaskedArray # happens if a1 isn't a numpy array, e.g. a float.  Then it's ok to just go on. #print "In arithmetic average of",mv1.id,"in two_pt_avg, encountered exception:",e #print "a1=",a1,type(a1),a1.dtype if hasattr(a1,'dtype') else None Updates the time-reduced data for a list of variables.  The reduced-time and averaged     variables are listed in redvars.  Its weights (for time averaging) are another variable,     redtime_wts.     (Each member of redvars should have the same the time axis.  Each member of newvars should have     the same time axis.  Normally it has an attribute wgts which is the same as redtime_wts.id.)     The new data is listed in newvars, and this list should correspond to redvars, e.g.     newvars[i].id==redvars[i].id.  The length of both should be equal and at least one.     Each variable is an MV.  Normally redvar and redtime_wts will be a FileVariable (required if     they might change) and newvar a TransientVariable.     If newvar needs any spatial reductions to match redvar, they should have been performed before     calling this function.     next_tbounds is the next time interval, used if newvar is defined on a time beyond redvar's     present time axis.  If next_tbounds==[], newvar will be ignored on such times.  Normally     next_tbounds will be set to [] when updating a climatology file which has been initialized.     The penultimate argument dt is used only in that dt=0 means that we are computing climatologies - hence     the new data's time axis must be adjusted before averaging the data into redvars.     The last argument is the time_weights global attribute of the data file, if any; it corresponds     to newvars.  This is expected to occur iff the data file is a climatology file written by     an earlier use of this module.     The optional argument force_scalar_avg argument is for testing and is passed on to two_pt_avg. # >>>> TO DO <<<< Ensure that each redvar, redtime_wts, newvar have consistent units # >>>> for the variable and for time.  Assert that they have the same shape, axes, etc. # formerly redvar was initialized here # The following two asserts express my assumption that the first index of the variable is time. # This is almost universal, but in the future I should generalize this code.  That would make # slicing more verbose, e.g. if time were changed from the first index to second then # v[j] would become v[:,j,:] (for a 2-D variable v). # partially-reduced time axis # some variables have no time axis # original time axis, from the new variable # The input data should have a time axis! # newtime_wts[j,i] is the weight applied to the data at time newtime[j] in computing # an average, reduced time for time newtime[ newtime_rti[i] ], 0<=i<2. # If newtime_rti[i]<0, that means the weight is 0. # Maximum number of data time intervals which could overlap with a single #               reduced-time interval.  We're unlikely to see more than 2. # First, extend redtime and redtime_bnds if necessary: # This should be moved to a separate function. # Extend the time axis to add a new time, time bounds, and weight.  With one more # silly step (see below), this will also have the effect of extending redvars along # the time axis. # Without this silly step, the data in redvars[iv] won't be expanded to match the # newly expanded time axis... # This also will do the job, but it looks like a lot of i/o: #   redvars[iv].parent.write(redvars[iv]) # It doesn't help to write redtime_wts or redtime_bnds.  You need to write a variable # with the same axes as redvars. # This doesn't do the job:  redvars[iv].parent.sync() # The weight of time newtime[j] is the part of its bounds which lie within some reduced- # time bounds.  We'll also need to remember the indices of the reduced-times for # which this is nonzero (there will be few of those, many reduced times total) #  This much simpler expression works if there is no overlap: #newtime_wts[j] = newtime_bnds[j][1] - newtime_bnds[j][0] # If k be unlimited, coding is more complicated. # This is the first point at which we decide whether the input file covers any times of interest, # e.g. for climatology of DJF, here is where we decide whether the file is a D,J,or F file. # Here kmax<0 if the file has no interesting time. # The weights from a climatology file make sense only when time is simply structured. # Otherwise, we don't know what to do. # This is the second point at which we decide whether the input file covers any times of interest, # e.g. for climatology of DJF, here is where we decide whether the file is a D,J,or F file. # Here i<0 if the file has no interesting time. # integer, any length, or string.  Time average makes no sense. # ...This assumes time axis is named 'time', but a full-blown check for a time # axis, e.g. via redvar.getTime(), is surprisingly expensive because of a # udunits inefficiency. # Although there's no time axis (hence index i is irrelevant), values may # differ from one file to the next, so we still have to do a time average. # Maybe the assignValue call will work for other shapes.  To do: try it, # this code will simplify if it works #won't work because redvar is a FileVariable # ... i.e., redvar[i] is initialized and is not entirely missing data # Don't miss this line, it's where the average is computed: # For averaging, unitialized is same as value=0 because redtime_wts is # initialized to 0 # uninitialized is same as value=0 Updates the time-reduced data for a several variables.  The reduced-time and averaged     variables are the list redvars.  Its weights (for time averaging) are another variable, redtime_wts.     (Each variable redvar of redvars has the same time axis, and it normally has an attribute wgts     which is the same as redtime_wts.id.)  redtime_bnds is a list of reduced time bounds.     Corresponging to it is the weights redtime_wts.     The new data will be read from the specified files, which should cover a time sequence in order.     We expect each member of redvar; and redtime_bnds[i] and redtime_wts[i], to be a FileVariable.     This function doesn't perform any spatial reductions.     The fun_next_tbounds argument is a function which will compute the next_tbounds argument of     update_time_avg, i.e. the next time interval (if any) to be appended to the bounds of the time     axis of redvars.     redfiles is a list of reduced-time files for output.  They should already be open as 'r+'.     The optional argument dt is used only in that dt=0 means that we are computing climatologies.     The optional argument force_scalar_avg argument is for testing and is passed on to two_pt_avg. #t1 = time.time() #t2 = time.time() # clone saves it in memory - faster # Compute tmin, tmax which represent the range of times for the data we computed with. # For regular model output files, these should come from the time bounds. # For a climatology file, they should come from the climatology attribute of time. #jfp testing # For testing done for <NAME>... #if hasattr(newvar,'dtype') and newvar.dtype=='float32': #    newvar = newvar.astype('float64') # excludes an ordinary number # Usually newvar is a TransientVariable, hence doesn't have the file as :parent. # We need the file to get associated variables, so: # newvar should be a TransientVariable.  But cdms2 doesn't work very well # for scalar variables.  It reads them as numpy numbers. # I know about this one. #redvars = redvars0 #dt = None   # if None, the "reduced" time bounds are exactly the same as in the input data # if >0, a fixed time step for partially time-reduced data #init_red_tbounds = [[]] # if [[]] or  [] or None, a reduced time interval is dt, the time step from #                       one interval to the next; thus all time is covered # a FileAxis # N is used to start the intervals off in the right year.  Note that this works only if calendar has a fixed-length year # example time interval for a single season - not all times #                              go into time-reduced data if dt=365 # For testing, print results... #test_time_avg( redfilename, varnames, datafilenames )