cdscan.py

#!/usr/bin/env python
from __future__ import print_function

import sys
import getopt
import cdms2
from cdms2.grid import lookupArray
from cdms2.axis import calendarToTag, tagToCalendar
from cdms2.cdmsobj import CdFromObject, CdString, CdScalar, CdFloat, CdDouble, CdShort, CdInt, CdLong
import numpy
import string
import cdtime
import os.path
import copy
from cdms2 import cdmsNode
import re
from functools import reduce
from cdms2.error import CDMSError
from collections import OrderedDict
from six import string_types
from cdms2.util import getenv_bool
usage = """Usage:
    cdscan [options] <files>

    Scan a list of files producing a CDMS dataset in XML representation. See Notes below
    for a more complete explanation.

Arguments:

    <files> is a list of file paths to scan. The files can be listed in any order, and may
    be in multiple directories.  A file may also be a CDML dataset (.xml or .cdml), in
    which case the dataset(s) and files are combined into a new dataset.

Options:

    -a alias_file: change variable names to the aliases defined in an alias file.
                   Each line of the alias file consists of two blank separated
                   fields: variable_id alias. 'variable_id' is the ID of the variable
                   in the file, and 'alias' is the name that will be substituted for
                   it in the output dataset. Only variables with entries in the alias_file
                   are renamed.

    -c calendar:   either "gregorian", "proleptic_gregorian", "julian", "noleap", or "360_day". Default:
                   "gregorian". This option should be used with caution, as it will
                   override any calendar information in the files.

    -d dataset_id: dataset identifier. Default: "none"

    -e newattr:    Add or modify attributes of a file, variable, or
           axis. The form of 'newattr' is either:

           'var.attr = value' to modify a variable or attribute, or
           '.attr = value' to modify a global (file) attribute.

           In either case, 'value' may be quoted to preserve spaces
           or force the attribute to be treated as a string. If
           'value' is not quoted and the first character is a
           digit, it is converted to integer or
           floating-point. This option does not modify the input
           datafiles. See notes and examples below.

    --exclude var,var,...
                   Exclude specified variables. The argument
                   is a comma-separated list of variables containing no blanks.
                   In contrast to --exclude-file, this skips the variables regardless
                   of the file(s) in which they are contained, but processes other
                   variables in the files.
                   Also see --include.

    --exclude-file pattern
                   Exclude files with a basename matching the regular expression pattern.
                   In contrast to --exclude, this skips the file entirely. Multiple patterns
                   may be listed by separating with vertical bars (e.g. abc|def ). Note
                   that the match is to the initial part of the basename. For example, the
                   pattern 'st' matches any basename starting with 'st'.

    -f file_list:  file containing a list of absolute data file names, one per
                   line. <files> arguments are ignored.

    --forecast     generate a description of a forecast dataset.
                   This is not compatible with the -i, -r, -t, or -l options.
                   A file can contain data for exactly one forecast; its
                   forecast_reference_time (aka run time, analysis time, starting time,
                   generating time, tau=0 time) is specified by the nbdate,nbsec variables.
                   Each file's time axis will be interpreted as the forecast_period (aka
                   tau, the interval from the forecast_reference_time to the current time)
                   regardless of its units, standard_name, or other attributes.

    -h:            print a help message.

    -i time_delta: scan time as a 'linear' dimension. This is useful if the time dimension
                   is very long. The argument is the time delta, a float or integer.  For
                   example, if the time delta is 6 hours, and the reference units are
                   "hours since xxxx", set the interval delta to 6.  The default value is
                   the difference of the first two timepoints.

    --ignore-open-error:
                   Ignore open errors. Print a warning and continue.

    --include var,var,...
                   Only include specified variables in the output. The argument
                   is a comma-separated list of variables containing no blanks.
                   Also see --exclude.

    --include-file pattern
                   Only include files with a basename matching the regular expression pattern.
                   In contrast to --include, this skips files entirely if they do not
                   match the pattern. Multiple patterns
                   may be listed by separating with vertical bars (e.g. abc|def ). Note
                   that the match is to the initial part of the basename. For example, the
                   pattern 'st' matches any basename starting with 'st'.

    -j:        scan time as a vector dimension. Time values are listed
           individually. Turns off the -i option.

    -l levels:     list of levels, comma-separated. Only specify if files are partitioned by
                   levels.

    -m levelid:    name of the vertical level dimension. The default is the name of the
                   vertical level dimension

    --notrim-lat:  Don't trim latitude values (in degrees) to the range [-90..90]. By default
           latitude values are trimmed.

    -p template:   Compatibility with pre-V3.0 datasets. 'cdimport -h' describes template strings.

    -q:            quiet mode

    -r time_units: time units of the form "<units> since yyyy-mm-dd hh:mi:ss", where
                   <units> is one of "year", "month", "day", "hour", "minute", "second".
                   Trailing fields may be omitted. The default is the units of the first
                   time dimension found.

    -s suffix_file: Append a suffix to variable names, depending on the directory
                   containing the data file.  This can be used to distinguish variables
                   having the same name but generated by different models or ensemble
                   runs. 'suffix_file' is the name of a file describing a mapping between
                   directories and suffixes.  Each line consists of two blank-separated
                   fields: 'directory' 'suffix'. Each file path is compared to the
                   directories in the suffix file. If the file path is in that directory
                   or a subdirectory, the corresponding suffix is appended to the variable
                   IDs in the file. If more than one such directory is found, the first
                   directory found is used. If no match is made, the variable ids are not
                   altered.  Regular expressions can be used: see the example in the Notes
                   section.

    -t timeid:     id of the partitioned time dimension. The default is the name of the time
                   dimension.

    --time-linear tzero,delta,units[,calendar]
                   Override the time dimensions(s) with a linear time dimension. The arguments are
                   a comma-separated list:

                   tzero is the initial time point, a floating-point value.
                   delta is the time delta, floating-point.
                   units are time units as specified in the [-r] option.
                   calendar is optional, and is specified as in the [-c] option. If omitted, it
                     defaults to the value specified by [-c], otherwise as specified in the file.

                   Example: --time-linear '0,1,months since 1980,noleap'

                   Note (6) compares this option with [-i] and [-r]

    --var-locate 'var,file_pattern':
                   Only scan a variable if the basename of the file matches the pattern. This
                   may be used to resolve duplicate variable errors. var and file_pattern are
                   separated by a comma, with no blanks.

                   var is the name of the variable
                   file_pattern is a regular expression following the Python re module syntax.e

                   Example: to scan variable ps from files starting with the string 'ps_':
                     --var-locate 'ps,ps_.*'

    -x xmlfile:    XML filename. By default, output is written to standard output.

Example:

    cdscan -c noleap -d test -x test.xml [uv]*.nc
    cdscan -d pcmdi_6h -i 0.25 -r 'days since 1979-1-1' *6h*.ctl

Notes:

    (1) The files can be in netCDF, GrADS/GRIB, HDF, or DRS format, and can be listed in
    any order. Most commonly, the files are the result of a single experiment, and the
    'partitioned' dimension is time. The time dimension of a variable is the coordinate
    variable having a name that starts with 'time' or having an attribute "axis='T'". If
    this is not the case, specify the time dimension with the -t option. The time
    dimension should be in the form supported by cdtime. If this is not the case (or to
    override them) use the -r option.

    (2) The basic form of the command is 'cdscan <files>'. By default, the time values are
    listed explicitly in the output XML. This can cause a problem if the time dimension is
    very long, say for 6-hourly data. To handle this the form 'cdscan -i delta <files>'
    may be used. This generates a compact time representation of the form <start, length,
    delta>. An exception is raised if the time dimension for a given file is not linear.

    (3) Another form of the command is 'cdscan -l lev1,lev2,..,levn <files>'. This asserts
    that the dataset is partitioned in both time and vertical level dimensions. The level
    dimension of a variable is the dimension having a name that starts with "lev", or
    having an attribute "axis=Z". If this is not the case, set the level name with the -m
    option.

    (4) An example of a suffix file:

    /exp/pr/ncar-a  _ncar-a
    /exp/pr/ecm-a   _ecm-a
    /exp/ta/ncar-a  _ncar-a
    /exp/ta/ecm-a   _ecm-a

    For all files in directory /exp/pr/ncar-a or a subdirectory, the corresponding
    variable ids will be appended with the suffix '_ncar-a'.  Regular expressions can be
    used, as defined in the Python 're' module. For example, The previous example can be
    replaced with the single line:

    /exp/[^/]*/([^/]*) _\g<1>

    Note the use of parentheses to delimit a group. The syntax \g<n> refers to the n-th
    group matched in the regular expression, with the first group being n=1. The string
    [^/]* matches any sequence of characters other than a forward slash.

    (5) Adding or modifying attributes with the -e option:

    time.units = "days since 1979-1-1"

    sets the units of all variables/axes to "Days since 1979-1-1". Note
    that since this is done before any other processing is done, it allows
    overriding of non-COARDS time units.

    .newattr=newvalue

    Set the global file attribute 'newattr' to 'newvalue'.

    (6) The [--time-linear] option overrides the time values in the file(s). The resulting
    dimension does not have any gaps. In contrast, the [-i], [-r] options use the specified
    time units (from [-r]), and calendar from [-c] if specified, to convert the file times
    to the new units. The resulting linear dimension may have gaps.

    In either case, the files are ordered by the time values in the files.

    The [--time-linear] option should be used with caution, as it is applied to all the time
    dimensions found.
"""

# Ensure that arrays are fully printed to XML files
numpy.set_printoptions(threshold=numpy.inf)

calendarMap = tagToCalendar

reverseCalendarMap = calendarToTag

attrPattern = re.compile(r'\s*(\w*)\.(\w+)\s*=\s*(.*)$')
cdms2.setNetcdfUseParallelFlag(0)


def timestamp():
    "Generate a timestamp."
    import time
    y, m, d, h, mi, s, w, dy, ds = time.gmtime(time.time())
    return "%d-%d-%d %d:%d:%d" % (y, m, d, h, mi, s)


def timeindex(value, units, basetime, delta, calendar):
    """ Calculate (t - basetime)/delu
    where t = reltime(value, units)
    and delu is the time interval (delta, delunits) (e.g., 1 month).
    """
    if units.find(" as ") == -1:
        tval = cdtime.reltime(value, units)
    else:
        tval = cdtime.abstime(value, units)
    newval = tval.torel(basetime, calendar)
    if delta is None:
        return newval.value
    else:
        return newval.value / delta


def combineKeys(mydict, typedict, timeIsLinear=0,
                referenceDelta=None, forecast=None):
    """ Combine dictionary keys into an axis.
    dict: (i,j) => (path, axisname)
    typedict is either timedict or levdict or fcdict.
    timeIsLinear is true iff time has a linear representation.
    referenceDelta is only used for error checks if timeIsLinear is true.
    """

    global verbose

    # Sort the projected time, level indices
    keys = list(OrderedDict(sorted(mydict.items())).keys())

    axislist = []
    prevend = None
    prevpath = None
    name0 = None
    compressPart = []
    partition = []
    previ = 0
    firstunits = ""
    prevvals = None
    coordToInd = {(None, None): (None, None)}
    linCoordToInd = {(None, None): (None, None)}
    iadj = None
    errorOccurred = 0
    for i0, i1 in keys:
        path, name = mydict[(i0, i1)]
        if name0 is None:
            name0 = name
        values, units, dummy = typedict[(path, name)]
        if firstunits == "":
            firstunits = units
        if prevend is not None and prevend >= i0:
            if prevend >= i1:
                if verbose:
                    print('Warning, file %s, dimension %s contains values in file %s' % (
                        prevpath, name, path), file=sys.stderr)
                if timeIsLinear:
                    iind = lookupArray(prevvals, values[0])
                    jind = lookupArray(prevvals, values[-1])
                else:
                    iind = lookupArray(prevvals, i0)
                    jind = lookupArray(prevvals, i1)
                if len(values) != (jind - iind + 1):
                    raise RuntimeError(
                        'Dimension %s in files %s [len(%s)=%d], %s [len(%s)=%d], is inconsistent' %
                        (name, prevpath, name, (jind - iind + 1), path, name, len(values)))
                coordToInd[(i0, i1)] = (iind, jind)

                prevspart, prevepart = partition[-1]
                linCoordToInd[(i0, i1)] = (
                    prevspart + iind, prevspart + jind + 1)
                continue
            else:                       # Fix partial overlap
                if timeIsLinear:
                    jind = lookupArray(prevvals, values[0])
                else:
                    jind = lookupArray(prevvals, i0)
                if verbose:
                    print('Warning, file %s, dimension %s overlaps file %s, value=%f' % (
                        prevpath, name, path, prevvals[jind]), file=sys.stderr)
                previ, prevj = compressPart[-1]
                prevj = previ + jind
                axislist[-1] = prevvals[0:jind]
                compressPart[-1] = (previ, prevj)
                coordToInd[(prevvals[0], prevvals[-1])] = (previ, prevj)
                previ = prevj

                prevspart, prevepart = partition[-1]
                prevepart = prevspart + jind
                partition[-1] = (prevspart, prevepart)
                linCoordToInd[(prevvals[0], prevvals[-1])
                              ] = (prevspart, prevepart)

        axislist.append(values)
        prevend = i1
        prevpath = path
        prevj = previ + len(values)
        compressPart.append((previ, prevj))
        coordToInd[(i0, i1)] = (previ, prevj)

        if iadj is None:                # partition has to start with 0
            iadj = int(i0)
        spart = int(i0) - iadj
        epart = int(i1) + 1 - iadj
        partition.append((spart, epart))
        linCoordToInd[(i0, i1)] = (spart, epart)
        if timeIsLinear and len(values) != (epart - spart):
            # Find the bad values
            diffs = values[1:] - values[:-1]
            badindices = numpy.compress(
                numpy.not_equal(diffs, referenceDelta), list(range(len(values))))
            badvalues = numpy.take(values, badindices)
            if verbose:
                print("Error: Missing values in %s after times: %s." % (path, str(badvalues)) +
                      " Set delta with the -i option or turn off linear mode with the -j option.",
                      file=sys.stderr)
            errorOccurred = 1

        prevvals = values
        previ = prevj

    fullaxis = numpy.ma.concatenate(axislist)
    return fullaxis, name0, compressPart, coordToInd, firstunits, partition, linCoordToInd, errorOccurred


def useKeys(mydict, typedict, timeIsLinear=0,
            referenceDelta=None, forecast=None):
    """ Use dictionary keys for an axis.  This is like combineKeys (same arguments, same return values,
    was written by simplifying combineKeys), but this doesn't do nearly so much because this is
    for an axis where there is no splitting across files, hence partitions are not needed.
    dict: (i,j) => (path, axisname)
    typedict is either timedict or levdict or fcdict.
    timeIsLinear is true iff time has a linear representation.
    referenceDelta is only used for error checks if timeIsLinear is true.
    """
    global verbose

    # Sort the projected time, level indices
    keys = list(OrderedDict(sorted(mydict.items())).keys())

    axislist = []
    name0 = None
#    compressPart = []
    compressPart = None
#    partition = []
    partition = None
#    previ = 0
    firstunits = None
#    coordToInd = {(None,None):(None,None)}
#    linCoordToInd = {(None,None):(None,None)}
    coordToInd = None
    linCoordToInd = None
    errorOccurred = 0
    for i0, i1 in keys:
        path, name = mydict[(i0, i1)]
        if name0 is None:
            name0 = name
        values, units, dummy = typedict[(path, name)]
        if firstunits is None:
            firstunits = units
        axislist.append(values)
#        prevj = previ+len(values)
#        coordToInd[(i0,i1)] = (previ, prevj)

    fullaxis = numpy.ma.concatenate(axislist)
    return fullaxis, name0, compressPart, coordToInd, firstunits, partition, linCoordToInd, errorOccurred


def copyDict(mydict):
    """Copy a dictionary-like object dict to a true dictionary"""
    result = {}
    for key in mydict.keys():
        result[key] = mydict[key]

    return result


def disambig(name, mydict, num, comparator, value):
    """ Make an unique name from name, wrt to the keys in dictionary dict.
    Try using num first. comparator(value,dict[name]) returns 0 if equal, 1 if not.
    """
    if name not in mydict or not comparator(value, mydict[name]):
        uniqname = name
    else:
        uniqname = '%s_%d' % (name, num)
        if uniqname in mydict and comparator(value, mydict[uniqname]):
            trial_name = uniqname
            for letter in string.lowercase:
                uniqname = '%s_%s' % (trial_name, letter)
                if uniqname not in mydict or not comparator(
                        value, mydict[uniqname]):
                    break
            else:
                raise BaseException(
                    'Cannot make axis name unique: ' + str(name))

    return uniqname


def compareaxes(axis1, axis2):
    """Return 0 if equal, 1 if not"""
    return ((len(axis1) != len(axis2)) or
            not numpy.ma.allclose(axis1[:], axis2[:]))


def comparedomains(domain1, domain2):
    """Return 0 if equal, 1 if not"""
    if len(domain1) != len(domain2):
        return 1
    for i in range(len(domain1)):
        item1 = domain1[i]
        item2 = domain2[i]
        if not isinstance(item1, type(item2)):
            return 1
        if isinstance(item1, string_types):
            return item1 != item2
        elif compareaxes(item1, item2):
            return 1
    return 0


def compareVarDictValues(val1, val2):
    return comparedomains(val1[0], val2[0])


def cleanupAttrs(attrs):
    for attname in attrs.keys():
        attval = attrs[attname]
        if isinstance(attval, numpy.ndarray):
            if attval.ndim == 0:
                continue
            elif len(attval) == 1:
                attrs[attname] = attval[0]
            else:
                attrs[attname] = str(attval)
    if 'missing_value' in attrs and attrs['missing_value'] is None:
        del attrs['missing_value']
    if 'ndim' in attrs:
        del attrs['ndim']


def validateAttrs(node):
    """Compare attributes against DTD."""

    global verbose

    if hasattr(node, 'datatype'):
        parenttype = node.datatype
    else:
        parenttype = None
    atts = node.getExternalDict()
    for attname in atts.keys():
        (attval, datatype) = atts[attname]  # (XML value, datatype)
        constraint = node.extra.get(attname)
        if constraint is not None:
            # (CdScalar|CdArray, required type)
            (scaletype, reqtype) = constraint
            if reqtype == CdFromObject:
                reqtype = parenttype
            if reqtype != datatype and datatype == CdString and scaletype == CdScalar:
                if reqtype in (CdFloat, CdDouble) and not isinstance(
                        attval, float):
                    try:
                        attval = float(attval)
                    except BaseException:
                        if verbose:
                            print("Warning: %s=%s should be a float, id=%s" % (
                                attname, attval, node.id), end=' ', file=sys.stderr)
                        try:
                            attval = int(attval)
                            attval = float(attval)
                            if verbose:
                                print("(Recasting)")
                            node.setExternalAttr(attname, attval)
                        except BaseException:
                            if attname in [
                                    'modulo', 'add_offset', 'scale_factor']:
                                if verbose:
                                    print("(Removing)")
                                attdict = node.getExternalDict()
                                del attdict[attname]
                            else:
                                if verbose:
                                    print("")
                elif reqtype in (CdShort, CdInt, CdLong) and not isinstance(attval, int):
                    try:
                        attval = int(attval)
                    except BaseException:
                        if verbose:
                            print("Warning: %s=%s should be an integer, id=%s" % (
                                attname, attval, node.id), end=' ', file=sys.stderr)
                        try:
                            attval = float(attval)
                            attval = int(attval)
                            if verbose:
                                print("(Recasting)")
                            node.setExternalAttr(attname, attval)
                        except BaseException:
                            if verbose:
                                print("")


def cloneWithLatCheck(axis):
    """Clone an axis, ensuring that latitudes (in degrees) are in the range [-90:90]"""

    global verbose
    global notrimlat

    axisvals = origvals = axis[:]
    if axis.isLatitude() and hasattr(
            axis, "units") and axis.units[0:6].lower() == "degree":
        if notrimlat == 0:
            axisvals = numpy.maximum(-90.0, numpy.minimum(90.0, axisvals))
        if not numpy.ma.allclose(axisvals, origvals) and verbose:
            print(
                "Warning: resetting latitude values: ",
                origvals,
                " to: ",
                axisvals,
                file=sys.stderr)

    b = axis.getBounds()
    mycopy = cdms2.createAxis(copy.copy(axisvals))
    mycopy.id = axis.id
    try:
        mycopy.setBounds(b)
    except CDMSError:
        b = mycopy.genGenericBounds()
        mycopy.setBounds(b)
    for k, v in axis.attributes.items():
        setattr(mycopy, k, v)
    return mycopy


def addAttrs(fobj, eattrs):
    """Add extra attributes to file/dataset fobj.
    eattrs has the form [(varid,attr,value), (varid,attr,value), ...]
    where if varid is '', set the global attribute."""
    for evar, eattr, evalue in eattrs:
        if evar == '':
            fobj.__dict__[eattr] = evalue
        else:
            varobj = fobj[evar]
            if varobj is not None:
                varobj.__dict__[eattr] = evalue


def setNodeDict(node, mydict):
    for key in mydict.keys():
        value = mydict[key]
        if (isinstance(value, numpy.integer) or
                isinstance(value, int)):
            datatype = CdLong
        elif (isinstance(value, numpy.floating) or isinstance(value, float)):
            datatype = CdDouble
        else:
            datatype = CdString
        node.attribute[key] = (value, datatype)


def initialize_filemap(filemap, timedict, levdict, timeid, extendDset, splitOnTime,
                       referenceTime, timeIsLinear, referenceDelta, splitOnLevel,
                       dirlen, overrideCalendar):
    # This function was formerly part of the body of "main".
    # Initialize filemap : varid => (tc0, tc1, lc0, lc1, path, timeid, levid)
    # where tc0 is the first time index relative to the reference time, tc1 the last,
    # lc0 is the first level, lc1 the last, path is the filename, timeid is the id
    # of the time dimension of the variable, levid is the id of the level dimension
    #
    # timedict : (path, timeid) => (timearray, timeunits, calendar)
    #
    # levdict : (path, levelid) => (levelarray, levelunits, None)
    #
    initfilemap = cdms2.dataset.parseFileMap(extendDset.cdms_filemap)
    dsetdirec = extendDset.directory
    for namelist, slicelist in initfilemap:
        for name in namelist:
            var = extendDset[name]
            timeaxis = var.getTime()
            if timeaxis is not None and not overrideCalendar:
                calendar = timeaxis.getCalendar()
            if splitOnTime and timeaxis is not None:
                if hasattr(timeaxis, 'name_in_file'):
                    timeid = timeaxis.name_in_file
                else:
                    timeid = timeaxis.id
                if referenceTime is None:
                    referenceTime = timeaxis.units
                if timeIsLinear in [None, 1]:
                    timeIsLinear = timeaxis.isLinear()
                    if timeIsLinear:
                        if len(timeaxis) > 1:
                            referenceDelta = timeaxis[1] - timeaxis[0]
                        else:
                            referenceDelta = 1.0
                    else:
                        referenceDelta = None
            else:
                timeid = None
            levelaxis = var.getLevel()
            if splitOnLevel and levelaxis is not None:
                if hasattr(levelaxis, 'name_in_file'):
                    levid = levelaxis.name_in_file
                else:
                    levid = levelaxis.id
            else:
                levid = None

            varmaplist = []
            for t0, t1, lev0, lev1, path in slicelist:
                fullpath = os.path.join(dsetdirec, path)
                basepath = fullpath[dirlen:]
                if t0 is not None:
                    tc0 = timeindex(
                        timeaxis[t0],
                        timeaxis.units,
                        referenceTime,
                        referenceDelta,
                        calendar)
                    tc1 = timeindex(timeaxis[t1 - 1],
                                    timeaxis.units,
                                    referenceTime,
                                    referenceDelta,
                                    calendar)
                    if (basepath, timeid, calendar) not in timedict:
                        values = timeaxis[t0:t1]
                        timedict[(basepath, timeid)] = (
                            values, timeaxis.units, calendar)
                else:
                    tc0 = tc1 = None
                if lev0 is not None:
                    lc0 = levelaxis[lev0]
                    lc1 = levelaxis[lev1 - 1]
                    if (basepath, levid, None) not in levdict:
                        values = levelaxis[lev0:lev1]
                        levdict[(basepath, levid)] = (
                            values, levelaxis.units, None)
                else:
                    lc0 = lc1 = None
                varmaplist.append(
                    (tc0, tc1, lc0, lc1, basepath, timeid, levid, calendar))
            if name in filemap:
                filemap[name].extend(varmaplist)
            else:
                filemap[name] = varmaplist

# ------------------------------------------------------------------------


verbose = 1


def main(argv):

    global verbose
    global notrimlat

    try:
        args, lastargs = getopt.getopt(
            argv[1:], "a:c:d:e:f:hi:jl:m:p:qr:s:t:x:",
            ["include=", "include-file=", "exclude=", "exclude-file=", "forecast", "time-linear=",
             "notrim-lat", "var-locate=", "ignore-open-error"])
    except getopt.error:
        print(sys.exc_info()[1], file=sys.stderr)
        print(usage, file=sys.stderr)
        sys.exit(0)

    calendar = None
    calenkey = None
    timeid = None
    levelid = None
    notrimlat = 0
    referenceTime = None
    referenceDelta = None
    readFromFile = 0
    splitOnTime = 1
    splitOnLevel = 0
    datasetid = "none"
    timeIsLinear = None
    writeToStdout = 1
    templatestr = None
    timeIsVector = None
    modelMapFile = None
    aliasMapFile = None
    overrideCalendar = 0
    extraAttrs = []
    extraDict = {}
    includeList = None
    excludeList = None
    overrideTimeLinear = None
    varLocate = None
    ignoreOpenError = False
    excludePattern = None
    includePattern = None
    forecast = False
    for flag, arg in args:
        if flag == '-a':
            aliasMapFile = arg
        elif flag == '-c':
            calenkey = arg.lower()
            calendar = calendarMap[calenkey]
            overrideCalendar = 1
        elif flag == '-d':
            datasetid = arg
        elif flag == '-e':
            matchObj = attrPattern.match(arg)
            if matchObj is None:
                raise RuntimeError(
                    "Expression must have form '[var].attr=value': %s" % arg)
            matchGroups = matchObj.groups()
            if len(matchGroups) != 3:
                raise RuntimeError(
                    "Expression must have form '[var].attr=value': %s" % arg)
            matchValue = matchGroups[2]
            if len(matchValue) > 0 and (matchValue[0].isdigit(
            ) or matchValue[0] in ['"', "'", "-", "+"]):  # "
                matcheval = eval(matchValue)
            else:
                matcheval = str(matchValue)
            extraAttrs.append((matchGroups[0], matchGroups[1], matcheval))
        elif flag == '--exclude':
            if arg[0] == '-':
                raise RuntimeError("--exclude option requires an argument")
            excludeList = arg.split(',')
        elif flag == '--exclude-file':
            excludePattern = arg
        elif flag == '-f':
            readFromFile = 1
            filelistpath = arg
        elif flag == '--forecast':  # experimental forecast mode
            forecast = True
            splitOnTime = 0
            splitOnLevel = 0
        elif flag == '-h':
            print(usage)
            sys.exit(0)
        elif flag == '-i':
            splitOnTime = 1
            referenceDelta = float(arg)
            timeIsLinear = 1
            timeIsVector = None
        elif flag == '--ignore-open-error':
            ignoreOpenError = True
        elif flag == '--include':
            if arg[0] == '-':
                raise RuntimeError("--include option requires an argument")
            includeList = arg.split(',')
        elif flag == '--include-file':
            includePattern = arg
        elif flag == '-j':
            timeIsVector = 1
            timeIsLinear = None
        elif flag == '-l':
            splitOnLevel = 1
            levelstr = arg.split(',')
            levellist = list(map(float, levelstr))
            levels = numpy.array(levellist)
            levels = numpy.sort(levels)
        elif flag == '-m':
            levelid = arg
            args.append(('-e', '%s.axis=Z' % levelid))  # Add axis=Z attribute
        elif flag == '--notrim-lat':
            notrimlat = 1
        elif flag == '-p':
            templatestr = arg
        elif flag == '-q':
            verbose = 0
        elif flag == '-r':
            splitOnTime = 1
            referenceTime = arg
        elif flag == '-s':
            modelMapFile = arg
        elif flag == '-t':
            splitOnTime = 1
            timeid = arg
            args.append(('-e', '%s.axis=T' % timeid))  # Add axis=T attribute
        elif flag == '--time-linear':
            targlist = arg.split(',')
            ttzero = float(targlist[0])
            tdelta = float(targlist[1])
            tunits = targlist[2].strip()
            if len(targlist) == 4:
                tcalendar = targlist[3].strip()

            else:
                tcalendar = None
            overrideTimeLinear = [ttzero, tdelta, tunits, tcalendar]
        elif flag == '--var-locate':
            if varLocate is None:
                varLocate = {}
            vname, pattern = arg.split(',')
            varLocate[vname] = pattern
        elif flag == '-x':
            writeToStdout = 0
            xmlpath = arg

    # If overriding time, process time as vector so that no gaps result
    if overrideTimeLinear is not None:
        timeIsVector = 1
        timeIsLinear = None
        if overrideCalendar == 1:
            overrideTimeLinear[3] = calenkey

    if verbose:
        print('Finding common directory ...')
    if readFromFile:
        f = open(filelistpath)
        lastargs = f.readlines()
        f.close()

    # Ignore blank paths
    realargs = []
    for arg in lastargs:
        sarg = arg.strip()
        if len(sarg) > 0:
            realargs.append(sarg)
    lastargs = realargs

    # Split lastargs into files and datasets
    fileargs = []
    dsetargs = []
    for arg in lastargs:
        base, suffix = os.path.splitext(arg)
        if suffix.lower() in ['.xml', '.cdml']:
            dsetargs.append(arg)
        else:
            fileargs.append(arg)

    # Generate a list of pathnames for datasets
    dsetfiles = []
    for path in dsetargs:
        dset = cdms2.open(path)
        if not hasattr(dset, 'cdms_filemap'):
            raise RuntimeError(
                'Dataset must have a cdms_filemap attribute: ' + path)
        if not hasattr(dset, 'directory'):
            raise RuntimeError(
                'Dataset must have a directory attribute: ' + path)
        dsetdirec = dset.directory
        initfilemap = cdms2.dataset.parseFileMap(dset.cdms_filemap)
        for namelist, slicelist in initfilemap:
            for t0, t1, lev0, lev1, path in slicelist:
                dsetfiles.append(os.path.join(dsetdirec, path))
    augmentedArgs = fileargs + dsetfiles

    # Find the common directory
    directory = os.path.commonprefix(augmentedArgs)
    firstpath = augmentedArgs[0][len(directory):]
    if not os.path.isfile(os.path.join(directory, firstpath)):
        dnew = os.path.dirname(directory)
        if len(dnew) > 0 and directory[len(dnew)] == '/':
            directory = dnew + '/'
        else:
            directory = dnew
    if verbose:
        print('Common directory:', directory)

    dirlen = len(directory)

    if templatestr is not None:
        if os.path.isabs(templatestr):
            templatestr = templatestr[dirlen:]
        templatere, ignore = cdms2.cdmsobj.templateToRegex(templatestr)
        template = re.compile(templatere + '$')
    else:
        template = None

    axisdict = {}
    vardict = {}
    filemap = {}
    timedict = {}
    levdict = {}
    fcdict = {}
    global_attrs = None
    fctau0 = None

    if modelMapFile is not None:
        mfile = open(modelMapFile)
        modelMap = {}
        modelDirs = []
        for line in mfile.readlines():
            mdirec, model = line.split()
            modelMap[mdirec] = model
            modelDirs.append(mdirec)
        mfile.close()

    if aliasMapFile is not None:
        afile = open(aliasMapFile)
        aliasMap = {}
        for line in afile.readlines():
            if line[0] not in ["'", '"']:  # "
                varid, alias = line.split()
            else:
                dummy, varid, alias = line.split(line[0])
                alias = alias.strip()
            aliasMap[varid] = alias
        afile.close()

    # Save extra attribute information for new axes
    for evar, eattr, evalue in extraAttrs:
        if evar == '':
            continue
        if evar in extraDict:
            curval = extraDict[evar]
            curval.append((eattr, evalue))
        else:
            extraDict[evar] = [(eattr, evalue)]

    # ------------------------------------------------------------------------
    # Initialize dictionaries if adding to an existing dataset
    if verbose and len(dsetargs) > 0:
        print('Scanning datasets ...')
    for extendPath in dsetargs:
        if verbose:
            print(extendPath)
        extendDset = cdms2.open(extendPath)

        # Add/modify attributes
        addAttrs(extendDset, extraAttrs)

        # Copy the global attribute dictionary if necessary. Note that copy.copy
        # can't be used here, since .attributes is now a 'fake' dictionary.
        if global_attrs is None:
            global_attrs = copyDict(extendDset.attributes)

        # Initialize filemap : varid => (tc0, tc1, lc0, lc1, path, timeid, levid)
        # where tc0 is the first time index relative to the reference time, tc1 the last,
        # lc0 is the first level, lc1 the last, path is the filename, timeid is the id
        # of the time dimension of the variable, levid is the id of the level dimension
        #
        # timedict : (path, timeid) => (timearray, timeunits, calendar)
        #
        # levdict : (path, levelid) => (levelarray, levelunits, None)
        #
        initialize_filemap(filemap, timedict, levdict, timeid, extendDset, splitOnTime,
                           referenceTime, timeIsLinear, referenceDelta, splitOnLevel,
                           dirlen, overrideCalendar)

        # axisdict : id => transient_axis
        #   for non-partitioned axes only
        #
        tempmap = {}
        for axis in list(extendDset.axes.values()):
            if not ((splitOnTime and (axis.isTime() or axis.id == timeid)) or
                    (splitOnLevel and (axis.isLevel() or axis.id == levelid))):
                axis = cloneWithLatCheck(axis)
                if axis.id in axisdict:
                    currentaxis = axisdict[axis.id]

                    # Check that the axis has the same length and values as the saved value. If not,
                    # create an unambiguous name in the axis dictionary.
                    if compareaxes(axis, currentaxis):
                        sepname = disambig(
                            axis.id, axisdict, len(axis), compareaxes, axis)
                        axis.name_in_file = axis.id
                        oldid = axis.id
                        axis.id = sepname
                        axisdict[sepname] = axis
                        tempmap[oldid] = sepname
                else:
                    axisdict[axis.id] = axis

        # vardict : varid => [domain, attributeDict, typecode]
        #   where domain = [axis_or_id, axis_or_id,...]
        #   and axis_or_id is the id of a partitioned dimension, or
        #   the transient axis object associated with a non-partitioned dimension
        #

        for var in list(extendDset.variables.values()):
            tempdomain = []
            for id in var.getAxisIds():
                if id in tempmap:
                    id = tempmap[id]
                if id in axisdict:
                    tempdomain.append(axisdict[id])
                else:
                    axis = extendDset[id]
                    if hasattr(axis, 'name_in_file'):
                        id = axis.name_in_file
                    tempdomain.append(id)
            varattrs = copyDict(var.attributes)
            vardict[var.id] = [tempdomain, varattrs, var.typecode()]

        extendDset.close()

        # end of loop "for extendPath in dsetargs"

    # ------------------------------------------------------------------------
    if verbose:
        print('Scanning files ...')

    boundsmap = {}                      # boundsmap : varid => timebounds_id
    boundsdict = {}                     # Same as vardict for time bounds
    for path in fileargs:
        path = path.strip()

        # Check if the path is included
        if includePattern is not None:
            base = os.path.basename(path)
            mobj = re.match(includePattern, base)
            if mobj is None:
                continue

        # Check if the path is excluded
        if excludePattern is not None:
            base = os.path.basename(path)
            mobj = re.match(excludePattern, base)
            if mobj is not None:
                continue

        if verbose:
            print(path)
        try:
            f = cdms2.open(path)
        except BaseException:
            if not ignoreOpenError:
                raise RuntimeError('Error opening file ' + path)
            else:
                print(
                    'Warning: cannot open file, skipping: %s' %
                    path, file=sys.stderr)
                continue

        # Add/modify attributes
        addAttrs(f, extraAttrs)

        # Determine the variable ID suffix, if any
        varsuffix = None
        if modelMapFile is not None:
            for direc in modelDirs:
                mo = re.match(direc, path)
                if mo is not None:
                    suffixPattern = modelMap[direc]

                    def gensuffix(m, mo=mo):
                        i = int(m.group(1))
                        return mo.group(i)
                    varsuffix = re.sub(r'\\g<(\d)>', gensuffix, suffixPattern)
                    break

        # Copy the global attribute dictionary if necessary. Note that copy.copy
        # can't be used here, since .attributes is now a 'fake' dictionary.
        if global_attrs is None:
            global_attrs = copyDict(f.attributes)

        basepath = path[dirlen:]
        if template is not None and template.match(basepath) is None:
            if verbose:
                print('Warning: path %s does not match template %s' % (
                    basepath, templatestr), file=sys.stderr)

        # Find time boundary variables
        boundsids = []
        if splitOnTime:
            tmpdict = {}
            for axisname in f.axes.keys():
                axis = f[axisname]
                # was if axis.isTime() and hasattr(axis, 'bounds'):
                if axis.isTime() and (axis.getBounds() is not None):
                    tmpdict[axis.bounds] = 1
            boundsids = list(tmpdict.keys())

        # For forecasts, get the time at which the forecast begins (tau=0) which
        # is nbdate,nbsec
        if forecast:
            # f('nbdate') is numpy.int32 which gets truncated
            nbdate = numpy.int(f('nbdate'))
            nbsec = f('nbsec')
            # hopefully nbsec<(seconds per day)=86400<100000
            fctau0 = nbdate * 100000 + nbsec
            fctau0time = cdtime.abstime(nbdate, "day as %Y%m%d")
            fctau0time = fctau0time.add(
                nbsec, cdtime.Seconds)  # fctau0 as type comptime
            fc_time_attrs = []

        varnames = list(f.variables.keys())

        # Try to force all axes to be included, but only small ones, length<100.
        # This section was motivated by a need to preserve the cloud axes isccp_prs,isccp_tau.
        # If we ever need to preserve longer axes as well, we could create one
        # variable per axis...
        crude_var_axes = [[ax[0] for ax in var.getDomain()]
                          for var in list(f.variables.values())]
        var_axes = set().union(*crude_var_axes)
        other_axes = list(set(f.axes.values()) - var_axes)
        if len(other_axes) > 0:
            other_axes = [ax for ax in other_axes if len(ax) < 100]
            other_axes.sort(key=(lambda ax: ax.id))
            axisvar = cdms2.createVariable(numpy.ones([len(ax) for ax in other_axes]),
                                           axes=other_axes, id='allaxesdummy')
            # all CdmsObj objects have this attribute, but for unknown
            axisvar.autoApiInfo = None
            # reasons datasetnode.dump() fails trying to dump this attribute's
            # default value (jfp)
            varnames.append(axisvar.id)
        # ...try to force all axes to be considered

        varnames = sorted(varnames)
        for varname in varnames:

            # If --var-locate is specified for the variable, match the basename
            # before processing
            if varLocate is not None and varname in varLocate:
                varpattern = varLocate[varname]
                base = os.path.basename(path)
                mobj = re.match(varpattern, base)
                if mobj is None:
                    continue

            # was var = f.variables[varname]
            if varname == 'allaxesdummy':
                var = axisvar
            else:
                var = f.variables[varname]

            # Reset the variable ID to any specified alias
            if aliasMapFile is not None:
                varalias = aliasMap.get(var.id)
                if varalias is not None:
                    var.name_in_file = var.id
                    var.id = varalias
                    varname = varalias

            # Append a suffix to the variable ID, if applicable
            if varsuffix is not None:
                if not hasattr(var, 'name_in_file'):
                    var.name_in_file = var.id
                var.id += varsuffix
                varname += varsuffix

            # [timestart, timeend, levstart, levend, path, timeid, levid, calendar, fctau0]
            varentry = [None] * 9
            varentry[4] = basepath
            varentry[8] = fctau0

            # Generate a temporary domain entry, and
            # create axis dictionary entries.
            domain = var.getDomain()
            if forecast:
                tempdomain = ['fctau0']
            else:
                # List of axis names and/or objects (if not partitioned)
                tempdomain = []
            for axis, start, length, truelen in domain:
                if (splitOnTime and (axis.isTime() or axis.id == timeid)) or \
                   (splitOnLevel and (axis.isLevel() or axis.id == levelid)):
                    tempdomain.append(axis.id)
                elif forecast and (axis.isTime() or axis.id == timeid):
                    # time axis isn't split but needs special treatment for
                    # forecasts
                    tempdomain.append(axis.id)
                    fc_time_attrs.append(axis.attributes)
                else:
                    axis = cloneWithLatCheck(axis)  # Transient copy
                    if axis.id in axisdict:
                        currentaxis = axisdict[axis.id]

                        # Check that the axis has the same length and values as the saved value. If not,
                        # create an unambiguous name in the axis dictionary.
                        if compareaxes(axis, currentaxis):
                            sepname = disambig(
                                axis.id, axisdict, len(axis), compareaxes, axis)
                            axis.name_in_file = axis.id
                            axis.id = sepname

                            # Fix boundary variable names if using suffixes.
                            if varsuffix is not None and hasattr(
                                    axis, 'bounds'):
                                axis.bounds += varsuffix
                            axisdict[sepname] = axis
                        else:
                            axis = currentaxis
                    else:
                        # Fix boundary variable names if using suffixes.
                        if varsuffix is not None and hasattr(axis, 'bounds'):
                            axis.bounds += varsuffix
                        axisdict[axis.id] = axis
                    tempdomain.append(axis)

            # Create a dictionary entry for the variable if not already there.
            if var.id in boundsids:
                boundsattrs = copyDict(var.attributes)
                boundsdict[var.id] = [tempdomain, boundsattrs, var.typecode()]
                continue                # Don't set a filemap entry until axes are sorted out
            elif var.id not in vardict:
                varattrs = copyDict(var.attributes)
                if varsuffix is not None or aliasMapFile is not None:
                    varattrs['name_in_file'] = var.name_in_file
                vardict[var.id] = [tempdomain, varattrs, var.typecode()]
            else:
                currentdomain, attrs, tcode = vardict[var.id]
                if comparedomains(currentdomain, tempdomain):
                    sepname = disambig(
                        var.id, vardict, var.size(), compareVarDictValues, (tempdomain, None))
                    saveid = var.id
                    varname = var.id = sepname
                    varattrs = copyDict(var.attributes)
                    var.name_in_file = varattrs['name_in_file'] = saveid
                    vardict[sepname] = [tempdomain, varattrs, var.typecode()]

            # Create a filemap entry for this variable/file, if split on time
            # or forecast
            axisids = [x[0].id for x in var.getDomain()]
            if splitOnTime or forecast:
                vartime = None
                if timeid is not None:
                    if timeid in axisids:
                        vartime = f.axes.get(timeid)
                    else:
                        if verbose:
                            print(
                                'Warning, time axis %s not found, -t option ignored' %
                                timeid, file=sys.stderr)
                if vartime is None:
                    vartime = var.getTime()
                if vartime is not None:
                    if not overrideCalendar:
                        calendar = vartime.getCalendar()
                    if referenceTime is None:
                        referenceTime = vartime.units
                    if verbose and not forecast:
                        print('Setting reference time units to', referenceTime)
                    if timeIsLinear is None and timeIsVector is None:
                        timeIsLinear = (referenceTime[0].lower().split() in
                                        ['hour',
                                         'hours',
                                         'minute',
                                         'minutes',
                                         'second',
                                         'seconds'])
                        if timeIsLinear and verbose:
                            print('Setting time representation to "linear"')  # '
                    if timeIsLinear and referenceDelta is None:
                        if len(vartime) > 1:
                            time1 = timeindex(
                                vartime[1], vartime.units, referenceTime, None, calendar)
                            time0 = timeindex(
                                vartime[0], vartime.units, referenceTime, None, calendar)
                            referenceDelta = time1 - time0
                        else:
                            referenceDelta = 1
                        if verbose:
                            print('Setting time delta to', referenceDelta)

#                    starttime = vartime[0]
#                    endtime = vartime[-1]
                    startindex = timeindex(
                        vartime[0],
                        vartime.units,
                        referenceTime,
                        referenceDelta,
                        calendar)
                    endindex = timeindex(
                        vartime[-1], vartime.units, referenceTime, referenceDelta, calendar)
                    if forecast:
                        # split on forecast, hence no split on time
                        varentry[0] = None
                        varentry[1] = None
                        referenceTime = None
                    else:
                        varentry[0] = startindex
                        varentry[1] = endindex
                    varentry[5] = vartime.id
                    varentry[7] = calendar

                    if (basepath, vartime.id) not in timedict:
                        values = vartime[:]
                        timedict[(basepath, vartime.id)] = (
                            values, vartime.units, calendar)

            if splitOnLevel:
                varlev = None
                if (levelid is not None) and (levelid in axisids):
                    varlev = f.axes.get(levelid)
                if varlev is None:
                    varlev = var.getLevel()
                if varlev is not None:
                    startlev = varlev[0]
                    if isinstance(startlev, numpy.ndarray):
                        startlev = startlev[0]
                    endlev = varlev[-1]
                    if isinstance(endlev, numpy.ndarray):
                        endlev = endlev[0]
                    varentry[2] = startlev
                    varentry[3] = endlev
                    varentry[6] = varlev.id

                    if (basepath, varlev.id, None) not in levdict:
                        values = varlev[:]
                        levdict[(basepath, varlev.id)] = (
                            values, varlev.units, None)

            if forecast:
                if (basepath, 'fctau0') not in fcdict:
                    fcdict[(basepath, 'fctau0')] = ([fctau0], "", None)

            if varname in filemap:
                filemap[varname].append(tuple(varentry))
            else:
                filemap[varname] = [tuple(varentry)]

            # Set boundsmap : varid => timebounds_id
            # was if splitOnTime and vartime is not None and hasattr(vartime,
            # "bounds") and not boundsmap.has_key(varname):
            if splitOnTime and vartime is not None and (vartime.getBounds() is not None) and\
                    varname not in boundsmap:
                boundsmap[varname] = vartime.bounds

            # End of loop "for varname in varnames"

        f.close()
        # End of loop "for path in fileargs"

    # ------------------------------------------------------------------------

    # Generate varindex, by combining variable names with
    # identical varentry values.
    varindex = []
#    varnames = sorted(filemap.keys())
    varnames = list(OrderedDict(sorted(filemap.items())).keys())

    for varname in varnames:
        varentry = sorted(filemap[varname])

        for varindexname, varindexvalue in varindex:
            if varentry == varindexvalue:
                varindexname.append(varname)
                break
        else:
            varindex.append(([varname], varentry))

    # If a variable is not a function of one of the partitioned dimensions,
    # no indexing is necessary: just read from the first file containing it.
    for varlist, slicelist in varindex:
        slice0 = slicelist[0]
        a, b, c, d, path0, timename, levname, calen, fctau0 = slice0
        if (a, b, c, d, fctau0) == (None, None, None, None, None):
            del slicelist[1:]

    # Change times to constant units
    sameCalendars = 1                   # True iff all time calendars are the same
    prevcal = None
    if forecast:
        # The data files' time axis is interpreted to be tau time, i.e. the forecast_period.
        # Find the axis, and remember it in timedict.
        for key in timedict.keys():
            values, units, calendar = timedict[key]
            if prevcal is not None and calendar != prevcal:
                sameCalendars = 0
            prevcal = calendar
            if units.find(" as ") == -1:
                time0 = cdtime.reltime(values[0], units)
            else:
                time0 = cdtime.abstime(values[0], units)
            offset = time0.torel(units, calendar).value  # normally will be 0
            values = values + offset - values[0]
            # Switch units from "normal" time such as "days since 2001-06-01"
            # to "basic" time such as "days", which makes sense for a
            # forecast_period.
            baslen = time0.units.find(' since ')
            basic_units = time0.units[0:baslen]  # e.g. 'days'
            fc_units = basic_units
            timedict[key] = (values, fc_units, calendar)
    else:       # splitOnTime is true
        for key in timedict.keys():
            values, units, calendar = timedict[key]
            if prevcal is not None and calendar != prevcal:
                sameCalendars = 0
            prevcal = calendar
            if units.find(" as ") == -1:
                time0 = cdtime.reltime(values[0], units)
            else:
                time0 = cdtime.abstime(values[0], units)
            offset = time0.torel(referenceTime, calendar).value
            values = values + offset - values[0]
            timedict[key] = (values, referenceTime, calendar)

    if sameCalendars and prevcal is not None:
        calenkey = reverseCalendarMap[calendar]

    if forecast:
        # For forecasts, make sure that the above has made all timedict values the same.
        # >>> It's conceivable that different forecasts will have different time (really, tau)
        # >>> axes.  If so, at this point we'll want to merge and mask all the time values, so
        # >>> that all variables can have the same time axis..  For now, just raise an error
        # >>> if there are time axis differences at this point.
        values0, units0, calendar0 = timedict[list(timedict.keys())[0]]

        timedict_same = all([((values0 == values1).all() and units0 == units1 and calendar0 == calendar1)
                             for (values1, units1, calendar1) in list(timedict.values())])
        if not timedict_same:
            raise CDMSError(
                'cdscan is confused about times for a forecast set')
        # Earlier we had saved all the time axis attributes.  Keep whatever
        # they have in common.
        fc_time_attr = fc_time_attrs[0]
        # go through all time attributes (each a dictionary)
        for fcta in fc_time_attrs:
            for attrn in fc_time_attr.keys():
                if attrn not in fcta:
                    # key attrn isn't in all time attributes
                    del fc_time_attr[attrn]
                elif fcta[attrn] != fc_time_attr[attrn]:
                    # not all time attributes have the same value for attrn
                    del fc_time_attr[attrn]
        # At this point fc_time_attr is the dictionary of those time attributes which are common to
        # all time axes encountered (in the context of a forecast dataset).
        # Finally, add the appropriate standard_name to it, if we haven't already gotten one from
        # the data file.  If the file has anything other than 'forecast_period', it's wrong, but
        # we'll stick with it anyway.
        if 'standard_name' not in fc_time_attr.keys():
            fc_time_attr['standard_name'] = 'forecast_period'

    # Create partitioned axes
    axes = []
    # varkey => (timeCoordToInd, levCoordToInd)
    masterCoordToInd = {}
    errorOccurred = 0
    for varlist, varentry in varindex:

        # Project time, level indices
        timeproj = {}
        levproj = {}
        fctproj = {}
        for time0, time1, lev0, lev1, path, timename, levname, calendar, fctau0 in varentry:
            if timename is not None:
                timeproj[(time0, time1)] = (path, timename)
            if levname is not None:
                try:
                    levproj[(lev0, lev1)] = (path, levname)
                except BaseException:
                    print('Cannot hash level %s range (%f,%f)' % (
                        levname, lev0, lev1), file=sys.stderr)
                    print(type(lev0), file=sys.stderr)
                    raise
            if fctau0 is not None:
                fctproj[(fctau0, fctau0)] = (path, 'fctau0')

        # and combine the projected indices into axes
        timeCoordToInd = None
        timelinCoordToInd = None
        if splitOnTime and timename is not None:
            fullaxis, name, partition, timeCoordToInd, units, opartition, timelinCoordToInd, errflag = \
                combineKeys(timeproj, timedict, timeIsLinear, referenceDelta)
            axes.append((varlist, fullaxis, name, partition, timeCoordToInd, units, opartition,
                         timelinCoordToInd, calendar))
            if errflag:
                errorOccurred = 1
        levCoordToInd = None
        if splitOnLevel and levname is not None:
            fullaxis, name, partition, levCoordToInd, units, opartition, levlinCoordToInd, errflag = \
                combineKeys(levproj, levdict)
            axes.append(
                (varlist,
                 fullaxis,
                 name,
                 partition,
                 levCoordToInd,
                 units,
                 opartition,
                 levlinCoordToInd,
                 None))
            if errflag:
                errorOccurred = 1
        fcCoordToInd = None
        if forecast:
            fullaxis, name, partition, fcCoordToInd, units, opartition, fclinCoordToInd, errflag = \
                combineKeys(fctproj, fcdict, forecast=forecast)
            axes.append(
                (varlist,
                 fullaxis,
                 name,
                 partition,
                 fcCoordToInd,
                 units,
                 opartition,
                 fclinCoordToInd,
                 None))
            if errflag:
                errorOccurred = 1
            # i.e., if time is in this variable's domain.
            if len(timeproj) > 0:
                # The useKeys call is like combineKeys, except that it's for a variable not partitioned
                # among files.  It just sets up axis data and (in the context of this variable loop)
                # propagates what's in timedict to every variable with time in
                # its domain.
                fullaxis, name, partition, timeCoordToInd, units, opartition, timelinCoordToInd, errflag = \
                    useKeys(timeproj, timedict, timeIsLinear, referenceDelta)
                axes.append((varlist, fullaxis, name, partition, timeCoordToInd, units, opartition,
                             timelinCoordToInd, calendar))
                if errflag:
                    errorOccurred = 1

        masterCoordToInd[varlist[0]] = (
            timeCoordToInd, levCoordToInd, timelinCoordToInd, fcCoordToInd)

    if errorOccurred:
        raise RuntimeError(
            'Error(s) determining axis values - see previous message(s)')

    # Eliminate duplicate axes
    axes2 = []
    for vlist1, axis1, name1, partition1, coordToInd1, units1, opartition1, linCoordToInd1, calen1 in axes:
        for vlist2, axis2, name2, partition2, coordToInd2, units2, opartition2, linCoordToInd2, calen2 in axes2:
            if len(axis1) == len(axis2) and name1 == name2 and partition1 == partition2 and units1 == units2 and \
                    numpy.ma.allclose(axis1, axis2) == 1 and calen1 == calen2:
                vlist2.extend(vlist1)
                break
        else:
            axes2.append((copy.copy(vlist1), axis1, name1, partition1, coordToInd1, units1, opartition1,
                          linCoordToInd1, calen1))

    # For each axis described by axis2, disambiguate its name, create the axis
    # object, etc.
    assignedBounds = {}
    for vlist, axis, name, partition, coordToInd, units, opartition, linCoordToInd, calendar in axes2:
        # print vlist, coordToInd
        uniqname = disambig(name, axisdict, len(axis), compareaxes, axis)
        axisobj = cdms2.createAxis(axis)
        axisobj.name_in_file = name
        axisobj.id = uniqname
        axisobj.units = units
        # For forecasts, give the time axis some saved attributes.
        if forecast and axisobj.isTime():
            for attr in fc_time_attr.keys():
                if not hasattr(axisobj, attr):
                    setattr(axisobj, attr, fc_time_attr[attr])
        if timeIsLinear and axisobj.isTime():
            axisobj.partition = numpy.ma.ravel(numpy.ma.array(opartition))
            axisobj.length = axisobj.partition[-1] - axisobj.partition[0]
            mopartition = numpy.array(opartition)
            partition_length = numpy.sum(mopartition[:, 1] - mopartition[:, 0])
            if partition_length < axisobj.length:
                axisobj.partition_length = partition_length
        elif partition is not None:
            axisobj.partition = numpy.ma.ravel(numpy.ma.array(partition))
        if axisobj.isTime():
            axisobj.calendar = reverseCalendarMap[calendar]
        # axisobj.reference_partition = str(numpy.ma.ravel(numpy.ma.array(opartition)))
        if uniqname not in axisdict:
            axisdict[uniqname] = axisobj
        for varname in vlist:
            domain, attributes, tcode = vardict[varname]
            for i in range(len(domain)):
                item = domain[i]
                if isinstance(item, string_types) and item == name:
                    domain[i] = axisobj

        # Add bounds variables to vardict, varindex
        if axisobj.isTime():
            # 'Representative' variable having this time axis
            reprVar = vlist[0]
            if reprVar in boundsmap:
                boundsname = boundsmap[reprVar]
                boundsinfo = boundsdict[boundsname]
                boundsattrs = boundsinfo[1]
                if uniqname != name:
                    boundsattrs['name_in_file'] = boundsname
                    boundsname = uniqname + '_bnds'
                if boundsname not in assignedBounds:
                    axisobj.bounds = boundsname
                    for varids, ranges in varindex:
                        if reprVar in varids:
                            varids.append(boundsname)
                    tmpdom = boundsinfo[0]
                    if isinstance(tmpdom[1], string_types):
                        bndsobj = tmpdom[0]
                        boundsdomain = (bndsobj, axisobj)
                    else:
                        bndsobj = tmpdom[1]
                        boundsdomain = (axisobj, bndsobj)
                    vardict[boundsname] = (
                        boundsdomain, boundsinfo[1], boundsinfo[2])
                    assignedBounds[boundsname] = 1

    # Collapse like indices in filemap. For example, transform
    # [x,[[0,10,-,-,file1], [0,10,-,-,file2]]] into
    # [x,[[0,10,-,-,file1]]]
    # This occurs for variables such as time boundaries, which are
    # often duplicated in different files.
    cdms_filemap_list = []
    duplicatevars = {}
    for varindexname, varindexvalue in varindex:
        timeCoordToInd, levCoordToInd, linCoordToInd, fcCoordToInd = masterCoordToInd[
            varindexname[0]]
        newslicedict = {}
        for time0, time1, lev0, lev1, path, timename, levname, calendar, fctau0 in varindexvalue:
            if timeCoordToInd is not None:
                if timeIsLinear:
                    i0, i1 = linCoordToInd[(time0, time1)]
                else:
                    i0, i1 = timeCoordToInd[(time0, time1)]
            else:
                i0 = i1 = None
            if levCoordToInd is not None:
                j0, j1 = levCoordToInd[(lev0, lev1)]
            else:
                j0 = j1 = None
            if (i0, i1, j0, j1, fctau0) in newslicedict:
                currentpath = newslicedict[(i0, i1, j0, j1, fctau0)]
                if tuple(varindexname) not in duplicatevars:
                    duplicatevars[tuple(varindexname)] = (currentpath, path)
            else:
                newslicedict[(i0, i1, j0, j1, fctau0)] = path
#        keys = sorted(newslicedict.keys())
        keys = list(OrderedDict(sorted(newslicedict.items())).keys())

        newslicelist = []
        for i0, i1, j0, j1, fctau0 in keys:
            path = newslicedict[(i0, i1, j0, j1, fctau0)]
            newslicelist.append([i0, i1, j0, j1, fctau0, path])
        cdms_filemap_list.append([varindexname, newslicelist])

    # Check if any duplicated variables are a function of longitude or latitude.
    # Raise an exception if so.
    illegalvars = []
    for varlist in duplicatevars.keys():
        for varname in varlist:
            if (excludeList is not None) and (varname in excludeList):
                continue
            dom, attrs, tcode = vardict[varname]
            for axisobj in dom:
                if axisobj.isLatitude() or axisobj.isLongitude():
                    path1, path2 = duplicatevars[varlist]
                    illegalvars.append((varname, path1, path2))
    if len(illegalvars) > 0:
        raise RuntimeError(
            "Variable '%s' is duplicated, and is a function of lat or lon: files %s, %s" %
            illegalvars[0])

    if verbose and len(list(duplicatevars.values())) > 0:
        print('Duplicate variables:', file=sys.stderr)
        for varlist in duplicatevars.keys():
            path1, path2 = duplicatevars[varlist]
            print('\t', varlist, '\t', path1, '\t', path2, file=sys.stderr)

    # Generate the cdms_filemap attribute
    cdms_filemap = str(cdms_filemap_list)
    cdms_filemap = cdms_filemap.replace(' ', '')
    cdms_filemap = cdms_filemap.replace('None', '-')
    cdms_filemap = cdms_filemap.replace('"', '')  # "
    cdms_filemap = cdms_filemap.replace("'", '')

    # Dump to XML
    datasetnode = cdmsNode.DatasetNode(datasetid)
    global_attrs['cdms_filemap'] = cdms_filemap
    global_attrs['directory'] = directory
    if sameCalendars and calenkey is not None:
        global_attrs['calendar'] = calenkey
    elif 'calendar' in global_attrs:
        del global_attrs['calendar']
    cleanupAttrs(global_attrs)
    # datasetnode.setExternalDict(global_attrs)
    setNodeDict(datasetnode, global_attrs)
    validateAttrs(datasetnode)

    timeWasOverridden = 0
    keys = list(OrderedDict(sorted(axisdict.items())).keys())
    for key in keys:
        axis = axisdict[key]
        tcode = axis.typecode()
        if tcode in [numpy.float32, numpy.float, numpy.int16,
                     numpy.int32, numpy.int, numpy.intc, numpy.int8]:
            tcode = numpy.sctype2char(tcode)
        cdtype = cdmsNode.NumericToCdType[tcode]
        node = cdmsNode.AxisNode(axis.id, len(axis), cdtype)

        # Override the time axis as a linear axis
        if axis.isTime() and (overrideTimeLinear is not None):
            ttzero = overrideTimeLinear[0]
            ttdelta = overrideTimeLinear[1]
            axis.units = overrideTimeLinear[2]
            if overrideTimeLinear[3] is None:
                axis.calendar = reverseCalendarMap[axis.getCalendar()]
            else:
                axis.calendar = overrideTimeLinear[3]
            linearnode = cdmsNode.LinearDataNode(ttzero, ttdelta, len(axis))
            node.setLinearData(linearnode)
            if verbose:
                if timeWasOverridden == 0:
                    print("Overriding values for axis '%s'" % axis.id)
                else:
                    print(
                        'Warning, overriding more than one time axis (%s)' %
                        axis.id, file=sys.stderr)
            timeWasOverridden = 1

        # Represent time as linear axis using time values in the file
        elif axis.isTime() and timeIsLinear:
            reference_length = axis.partition[-1] - axis.partition[0]
            linearnode = cdmsNode.LinearDataNode(
                axis[0], referenceDelta, reference_length)
            node.setLinearData(linearnode)
        else:
            try:
                node.setData(axis[:])
            except cdms2.cdmsNode.NotMonotonicError:
                if verbose:
                    print('Warning: Axis values for axis %s are not monotonic:' % axis.id, axis[
                        :], file=sys.stderr)
                    print('Warning: Resetting axis %s values to:' % axis.id, numpy.arange(
                        len(axis)), file=sys.stderr)
                node.setData(numpy.arange(len(axis)))
        axisattrs = copyDict(axis.attributes)

        # Make sure that new axes have attribute mods
        if key in extraDict:
            for eattr, evalue in extraDict[key]:
                axisattrs[eattr] = evalue
        cleanupAttrs(axisattrs)
        # node.setExternalDict(axisattrs)
        setNodeDict(node, axisattrs)
        validateAttrs(node)
        datasetnode.addId(axis.id, node)

#    keys = sorted(vardict.keys())
#    keys = vardict.keys()
#    keys = list(keys).sort()
    keys = list(OrderedDict(sorted(vardict.items())).keys())
    for key in keys:
        if (includeList is not None) and (key not in includeList):
            continue
        if (excludeList is not None) and (key in excludeList):
            continue
        domain, attrs, tcode = vardict[key]
        if tcode in [numpy.float32, numpy.float, numpy.int16,
                     numpy.int32, numpy.int, numpy.intc, numpy.int8]:
            tcode = numpy.sctype2char(tcode)
        domainNode = cdmsNode.DomainNode()
        cdtype = cdmsNode.NumericToCdType[tcode]
        node = cdmsNode.VariableNode(key, cdtype, domainNode)
        cleanupAttrs(attrs)
        # node.setExternalDict(attrs)
        setNodeDict(node, attrs)
        validateAttrs(node)
        for axis in domain:
            if hasattr(axis, 'length'):
                length = axis.length
            else:
                length = len(axis)
            try:
                elemnode = cdmsNode.DomElemNode(axis.id, 0, length)
            except AttributeError:
                print('Axis %s for variable %s does not have attribute "id"' % (
                    repr(axis), key), file=sys.stderr)
            if hasattr(axis, 'partition_length'):
                elemnode.setExternalAttr(
                    'partition_length', axis.partition_length)
            domainNode.add(elemnode)
        datasetnode.addId(key, node)

    # Add the Conventions attribute if not present
    conventions = datasetnode.getExternalAttr('Conventions')
    if conventions is None:
        datasetnode.setExternalAttr('Conventions', '')
    if templatestr is not None:
        datasetnode.setExternalAttr('template', templatestr)

    # Add/modify history
    history = datasetnode.getExternalAttr('history')
    if history is None:
        history = ""
    stringargv = reduce(lambda x, y: x + ' ' + y, argv)
    stringtime = "\n[%s] " % timestamp()
    if len(stringargv) <= 256:
        history += stringtime + stringargv
    else:
        history += stringtime + stringargv[:256] + " ..."
    datasetnode.setExternalAttr('history', history)

    # datasetnode.validate()
    if writeToStdout:
        datasetnode.dump()
    else:
        datasetnode.dump(xmlpath)
        if verbose:
            print(xmlpath, 'written')


# ------------------------------------------------------------------------
if __name__ == '__main__':
    main(sys.argv)
    mpi_disabled = getenv_bool("CDMS_NO_MPI", "False")
    try:
        # skip trying to load mpi4py module
        if mpi_disabled:
            raise Exception()

        from mpi4py import MPI
        comm = MPI.Comm.Get_parent()
        comm.send('done', dest=0)
    except BaseException:
        pass