from numpy import argsort as numpy_argsort
from numpy import dtype as numpy_dtype
from numpy import sort as numpy_sort
from operator import attrgetter, itemgetter
from itertools import izip
#from .ancillaryvariables import AncillaryVariables
#from .comparison import gt
from .coordinate import AuxiliaryCoordinate
from .coordinatereference import CoordinateReference
from .field import Field, FieldList
#from .fieldlist import FieldList
from .query import gt
from .units import Units
from .functions import (flat, RTOL, ATOL, equals, hash_array, allclose)
from .functions import inspect as cf_inspect
from .data.data import Data
from .data.filearray import FileArray
_dtype_float = numpy_dtype(float)
## --------------------------------------------------------------------
## Global properties, as defined in Appendix A of the CF conventions.
## --------------------------------------------------------------------
#_global_properties = set(('comment',
# 'Conventions',
# 'history',
# 'institution',
# 'references',
# 'source',
# 'title',
# ))
# --------------------------------------------------------------------
# Data variable properties, as defined in Appendix A of the CF
# conventions, without those which are not simple. And less
# 'long_name'.
# --------------------------------------------------------------------
_signature_properties = set(('add_offset',
'calendar',
'cell_methods',
'_FillValue',
'flag_masks',
'flag_meanings',
'flag_values',
'missing_value',
'scale_factor',
'standard_error_multiplier',
'standard_name',
'units',
'valid_max',
'valid_min',
'valid_range',
))
#_standard_properties = _data_properties.union(_global_properties)
_no_units = Units()
class _HFLCache(object):
'''
A cache for coordinate and cell measure hashes, first and last values
and first and last cell bounds
'''
def __init__(self):
self.hash = {}
self.fl = {}
self.flb = {}
#--- End: def
def inspect(self):
'''
Inspect the object for debugging.
.. seealso:: `cf.inspect`
:Returns:
None
:Examples:
>>> f.inspect()
'''
print cf_inspect(self)
#--- End: def
#--- End: class
class _Meta(object):
'''
A summary of a field.
This object contains everything you need to know in order to aggregate
the field.
'''
#
_canonical_units = {}
#
_canonical_cell_methods = []
def __init__(self, f,
rtol=None, atol=None,
info=0,
relaxed_units=False,
allow_no_identity=False,
respect_valid=False,
equal_all=False,
exist_all=False,
equal=None,
exist=None,
ignore=None,
dimension=(),
relaxed_identities=False,
ncvar_identities=False):
'''
**initialization**
:Parameters:
f : cf.Field
info : int, optional
See the `aggregate` function for details.
relaxed_units : bool, optional
See the `aggregate` function for details.
allow_no_identity : bool, optional
See the `aggregate` function for details.
rtol : float, optional
See the `aggregate` function for details.
atol : float, optional
See the `aggregate` function for details.
dimension : (sequence of) str, optional
See the `aggregate` function for details.
:Examples:
'''
self._nonzero = False
self.info = info
self.sort_indices = {}
self.sort_keys = {}
self.cell_values = False
self.message = ''
strict_identities = not (relaxed_identities or ncvar_identities)
self.strict_identities = strict_identities
self.ncvar_identities = ncvar_identities
# Initialize the flag which indicates whether or not this
# field has already been aggregated
self.aggregated_field = False
# ------------------------------------------------------------
# Field
# ------------------------------------------------------------
self.field = f
self._hasData = f._hasData
self.identity = f.name(identity=strict_identities,
ncvar=ncvar_identities)
# ------------------------------------------------------------
#
# ------------------------------------------------------------
signature_override = getattr(f, 'aggregate', None)
if signature_override is not None:
self.signature = signature_override
self._nonzero = True
return
if self.identity is None:
if not allow_no_identity and self._hasData:
if info:
self.message = \
"no identity; consider setting relaxed_identities"
return
elif not self._hasData:
if info:
self.message = \
"no data array"
return
#--- End: if
domain = f.domain
items = domain.items
# ------------------------------------------------------------
# Promote selected properties to 1-d, size 1 auxiliary
# coordinates
# ------------------------------------------------------------
for property in dimension:
value = f.getprop(property, None)
if value is None:
continue
aux_coord = AuxiliaryCoordinate(properties={'long_name': property},
attributes={'id': property},
data=Data([value], units=''),
copy=False)
axis = domain.new_axis_identifier()
domain.insert_aux(aux_coord, axes=[axis], copy=False)
f.delprop(property) ### dch COPY issue?
#--- End: for
self.units = self.canonical_units(f, self.identity,
relaxed_units=relaxed_units)
self.cell_methods = self.canonical_cell_methods(rtol=rtol, atol=atol)
# ------------------------------------------------------------
# Formula_terms
# ------------------------------------------------------------
coordrefs = items(role='r')
if not coordrefs:
self.coordrefs = ()
self.coordref_signatures = ()
else:
self.coordrefs = coordrefs.values()
self.coordref_signatures = sorted([cr.structural_signature()
for cr in self.coordrefs])
for s in self.coordref_signatures:
if s[0] is None:
if info:
self.messsage = \
"%r field can't be aggregated due to it having an unidentifiable coordinate reference" % \
f.name('')
return
#--- End: if
#--- End: if
# ------------------------------------------------------------
# Ancillary variables
# ------------------------------------------------------------
if not self.set_ancillary_variables():
return
# ------------------------------------------------------------
# Coordinate and cell measure arrays
# ------------------------------------------------------------
self.hash_values = {}
self.first_values = {}
self.last_values = {}
self.first_bounds = {}
self.last_bounds = {}
# Map axis canonical identities to their domain identifiers
#
# For example: {'time': 'dim2'}
self.id_to_axis = {}
# Map domain axis identifiers to their canonical identities
#
# For example: {'dim2': 'time'}
self.axis_to_id = {}
# Dictionaries mapping domain auxiliary coordinate identifiers
# to their auxiliary coordiante objects
aux_1d = items(role='a', ndim=1)
# A set containing the identity of each domain coordinate
#
# For example: set(['time', 'height', 'latitude',
# 'longitude'])
self.all_coord_identities = set()
self.axis = {}
for axis in domain._axes_sizes:
# List some information about each 1-d coordinate which
# spans this axis. The order of elements is arbitrary, as
# ultimately it will get sorted by each element's 'name'
# key values.
#
# For example: [{'name': 'time', 'key': 'dim0', 'units':
# <CF Units: ...>}, {'name': 'forecast_ref_time', 'key':
# 'aux0', 'units': <CF Units: ...>}]
info_dim = []
dim_coord = domain.item(axis)
if dim_coord is not None:
# ----------------------------------------------------
# 1-d dimension coordinate
# ----------------------------------------------------
identity = self.coord_has_identity_and_data(dim_coord)
if identity is None:
return
# Find the canonical units for this dimension
# coordinate
units = self.canonical_units(dim_coord, identity,
relaxed_units=relaxed_units)
info_dim.append(
{'identity' : identity,
'key' : axis,
'units' : units,
'hasbounds': dim_coord._hasbounds,
'coordrefs': self.find_coordrefs(axis, dim_coord)})
#--- End: if
# Find the 1-d auxiliary coordinates which span this axis
aux_coords = {}
for aux in aux_1d.keys():
if axis in domain.item_axes(aux): #dimensions[aux]:
aux_coords[aux] = aux_1d.pop(aux)
#--- End: for
info_aux = []
for aux, aux_coord in aux_coords.iteritems():
# ----------------------------------------------------
# 1-d auxiliary coordinate
# ----------------------------------------------------
identity = self.coord_has_identity_and_data(aux_coord)
if identity is None:
return
# Find the canonical units for this 1-d auxiliary
# coordinate
units = self.canonical_units(aux_coord, identity,
relaxed_units=relaxed_units)
info_aux.append(
{'identity' : identity,
'key' : aux,
'units' : units,
'hasbounds': aux_coord._hasbounds,
'coordrefs': self.find_coordrefs(aux, aux_coord)})
#--- End: for
# Sort the 1-d auxiliary coordinate information
info_aux.sort(key=itemgetter('identity'))
# Prepend the dimension coordinate information to the
# auxiliary coordinate information
info_1d_coord = info_dim + info_aux
if not info_1d_coord:
if info:
self.message ="\
axis has no one dimensional or scalar coordinates"
#"% field can't be aggregated due to an axis having no 1-d coordinates" %
#f.name(''))
return
#--- End: if
# Find the canonical identity for this axis
identity = info_1d_coord[0]['identity']
self.axis[identity] = \
{'ids' : tuple([i['identity'] for i in info_1d_coord]),
'keys' : tuple([i['key'] for i in info_1d_coord]),
'units' : tuple([i['units'] for i in info_1d_coord]),
'hasbounds': tuple([i['hasbounds'] for i in info_1d_coord]),
'coordrefs': tuple([i['coordrefs'] for i in info_1d_coord])}
if info_dim:
self.axis[identity]['dim_coord_index'] = 0
else:
self.axis[identity]['dim_coord_index'] = None
self.id_to_axis[identity] = axis
self.axis_to_id[axis] = identity
#--- End: for
# Create a sorted list of the axes' canonical identities
#
# For example: ['latitude', 'longitude', 'time']
self.axis_ids = sorted(self.axis)
# ------------------------------------------------------------
# N-d auxiliary coordinates
# ------------------------------------------------------------
self.nd_aux = {}
for aux, nd_aux_coord in items(role='a', ndim=gt(1)).iteritems():
# Find this N-d auxiliary coordinate's identity
identity = self.coord_has_identity_and_data(nd_aux_coord)
if identity is None:
return
# Find the canonical units
units = self.canonical_units(nd_aux_coord, identity,
relaxed_units=relaxed_units)
# Find axes' canonical identities
axes = [self.axis_to_id[axis] for axis in domain.item_axes(aux)]
axes = tuple(sorted(axes))
self.nd_aux[identity] = {
'key' : aux,
'units' : units,
'axes' : axes,
'hasbounds': nd_aux_coord._hasbounds,
'coordrefs': self.find_coordrefs(aux, nd_aux_coord)}
#--- End: for
# ------------------------------------------------------------
# Cell measures
# ------------------------------------------------------------
self.msr = {}
info_msr = {}
for key, msr in items(role='m').iteritems():
if not self.cell_measure_has_data_and_units(msr):
return
# Find the canonical units for this cell measure
units = self.canonical_units(msr,
msr.name(identity=strict_identities,
ncvar=ncvar_identities),
relaxed_units=relaxed_units)
# Find axes' canonical identities
axes = [self.axis_to_id[axis] for axis in domain.item_axes(key)]
axes = tuple(sorted(axes))
if units in info_msr:
# Check for ambiguous cell measures, i.e. those which
# have the same units and span the same axes.
for value in info_msr[units]:
if axes == value['axes']:
if info:
self.message = \
"duplicate %r cell measures" % msr.name('')
return
else:
info_msr[units] = []
#--- End: if
info_msr[units].append({'key' : key,
'axes': axes})
#--- End: for
# For each cell measure's canonical units, sort the
# information by axis identities.
for units, value in info_msr.iteritems():
value.sort(key=itemgetter('axes'))
self.msr[units] = {'keys': tuple([v['key'] for v in value]),
'axes': tuple([v['axes'] for v in value])}
#--- End: for
# ------------------------------------------------------------
# Properties and attributes
# ------------------------------------------------------------
if not (equal or exist or equal_all or exist_all):
self.properties = ()
else:
properties = f.properties
for p in ignore:
properties.pop(p, None)
if equal:
eq = dict([(p, properties[p]) for p in equal
if p in properties])
else:
eq = {}
if exist:
ex = [p for p in exist if p in properties]
else:
ex = []
eq_all = {}
ex_all = []
if equal_all or exist_all:
if equal_all:
if not equal and not exist:
eq_all = properties
elif equal and exist:
eq_all = dict([(p, properties[p]) for p in properties
if p not in ex and p not in eq])
elif equal:
eq_all = dict([(p, properties[p]) for p in properties
if p not in eq])
elif exist:
eq_all = dict([(p, properties[p]) for p in properties
if p not in ex])
elif exist_all:
if not equal and not exist:
ex_all = list(properties)
elif equal and exist:
ex_all = [p for p in properties
if p not in ex and p not in eq]
elif equal:
ex_all = [p for p in properties if p not in eq]
elif exist:
ex_all = [p for p in properties if p not in ex]
#--- End: if
self.properties = tuple(sorted(ex_all + ex +
eq_all.items() + eq.items()))
#--- End: if
# Attributes
self.attributes = set(('file',))
# ----------------------------------------------------------------
# Still here? Then create the structural signature.
# ----------------------------------------------------------------
self.respect_valid = respect_valid
self.structural_signature()
# Initialize the flag which indicates whether or not this
# field has already been aggregated
self.aggregated_field = False
self.sort_indices = {}
self.sort_keys = {}
# Finally, set the object to True
self._nonzero = True
#--- End: def
def __nonzero__(self):
'''
x.__nonzero__() <==> bool(x)
'''
return self._nonzero
#--- End: if
def __repr__(self):
'''
x.__repr__() <==> repr(x)
'''
return '<CF %s: %r>' % (self.__class__.__name__,
getattr(self, 'field', None))
#--- End: def
def __str__(self):
'''
x.__str__() <==> str(x)
'''
strings = []
for attr in sorted(self.__dict__):
strings.append('%s.%s = %r' % (self.__class__.__name__, attr,
getattr(self, attr)))
return '\n'.join(strings)
#--- End: def
def coordinate_values(self):
'''
'''
string = ['First cell: '+str(self.first_values)]
string.append('Last cell: '+str(self.last_values))
string.append('First bounds: ' +str(self.first_bounds))
string.append('Last bounds: ' +str(self.last_bounds))
return '\n'.join(string)
#--- End: def
def copy(self):
'''
'''
new = _Meta.__new__(_Meta)
new.__dict__ = self.__dict__.copy()
new.field = new.field.copy()
return new
def canonical_units(self, variable, identity, relaxed_units=False):
'''
Updates the `_canonical_units` attribute.
:Parameters:
variable : cf.Variable
identity : str
relaxed_units : bool
See the `cf.aggregate` for details.
:Returns:
out : cf.Units or None
:Examples:
'''
var_units = variable.Units
_canonical_units = self._canonical_units
if identity in _canonical_units:
if var_units:
for u in _canonical_units[identity]:
if var_units.equivalent(u):
return u
#--- End: for
# Still here?
_canonical_units[identity].append(var_units)
elif relaxed_units or variable.dtype.kind == 'S':
var_units = _no_units
else:
if var_units:
_canonical_units[identity] = [var_units]
elif relaxed_units or variable.dtype.kind == 'S':
var_units = _no_units
#--- End: if
# Still here?
return var_units
#--- End: def
def canonical_cell_methods(self, rtol=None, atol=None):
'''
Updates the `_canonical_cell_methods` attribute.
:Parameters:
atol : float
rtol : float
:Returns:
out : cf.CellMethods or None
:Examples:
'''
cell_methods = getattr(self.field, 'cell_methods', None)
if cell_methods is None:
return None
_canonical_cell_methods = self._canonical_cell_methods
for cm in _canonical_cell_methods:
if cell_methods.equivalent(cm, rtol=rtol, atol=atol):
return cm
#--- End: for
# Still here?
_canonical_cell_methods.append(cell_methods)
return cell_methods
#--- End: def
def cell_measure_has_data_and_units(self, msr):
'''
:Parameters:
msr : cf.CellMeasure
:Returns:
out : bool
:Examples:
'''
if not msr.Units:
if self.info:
self.message = \
"%r cell measure has no units" % msr.name('')
return
if not msr._hasData:
if self.info:
self.message = \
"%r cell measure has no data" % msr.name('')
return
return True
#--- End: def
def coord_has_identity_and_data(self, coord):
'''
:Parameters:
coord : cf.Coordinate
:Returns:
out : str or None
The coordinate object's identity, or None if there is no
identity and/or no data.
:Examples:
'''
identity = coord.name(identity=self.strict_identities,
ncvar=self.ncvar_identities)
if identity is None:
# Coordinate has no identity, but it may have a recognised
# axis.
for ctype in ('T', 'X', 'Y', 'Z'):
if getattr(coord, ctype):
identity = ctype
break
#--- End: if
if identity is not None:
all_coord_identities = self.all_coord_identities
if identity in all_coord_identities:
if self.info:
self.message = \
"multiple %r coordinates" % identity
return None
#--- End: if
if coord._hasData or (coord._hasbounds and coord.bounds._hasData):
all_coord_identities.add(identity)
return identity
#--- End: if
if self.info:
self.message = \
"%r coordinate has no identity or no data" % coord.name('')
return None
#--- End: def
def set_ancillary_variables(self):
'''
:Returns:
out : dict or None
:Examples:
'''
f_ancillary_variables = getattr(self.field, 'ancillary_variables',
None)
if f_ancillary_variables is None:
self.ancillary_variables = {}
return True
ancillary_variables = {}
for av in f_ancillary_variables:
identity = av.name(identity=self.strict_identities,
ncvar=self.ncvar_identities)
if identity in ancillary_variables:
if self.info:
self.message = \
"multiple %r ancillary variables" % av.name('')
return None
#--- End: if
ancillary_variables[identity] = av
#--- End: for
self.ancillary_variables = ancillary_variables
return True
#--- End: def
def print_info(self, info, signature=True):
'''
:Parameters:
m : _Meta
info : int
'''
if info >= 2:
if signature:
print 'STRUCTURAL SIGNATURE:\n', self.string_structural_signature()
if self.cell_values:
print 'CANONICAL COORDINATES:\n', self.coordinate_values()
if info >= 3:
print 'COMPLETE AGGREGATION METADATA:\n', self
#--- End: def
def string_structural_signature(self):
'''
'''
keys = ('Identity',
'Units',
'Cell methods',
'Data',
'Properties',
'standard error multiplier',
'valid_min',
'valid_max',
'valid_range',
'Flags',
'Ancillary variables',
'Coordinate reference systems',
'1-d coordinates',
'Dimension coordinates',
'N-d coordinates',
'Cell measures',
)
string = []
# d = {}
for key, value in zip(keys[:], self.signature[:]):
if not (value == () or value is None):
string.append('%s: %r' % (key, value))
# return '{'+'\n'.join(string)+'}'
return '\n'.join(string)
#--- End: def
def structural_signature(self):
'''
:Returns:
out : tuple
:Examples:
'''
f = self.field
# Initialize the structual signature with:
#
# * the identity
# * the canonical units
# * the canonical cell methods
# * whether or not there is a data array
signature = [self.identity, self.units, self.cell_methods, self._hasData]
signature_append = signature.append
# Properties
signature_append(self.properties)
# standard_error_multiplier
signature_append(f.getprop('standard_error_multiplier', None))
# valid_min, valid_max, valid_range
if self.respect_valid:
signature.extend((f.getprop('valid_min' , None),
f.getprop('valid_max' , None),
f.getprop('valid_range', None)))
else:
signature.extend((None, None, None))
# Flags
signature_append(getattr(f, 'Flags', None))
# Add ancillary variables
if self.ancillary_variables:
signature_append(tuple(sorted(self.ancillary_variables)))
else:
signature_append(None)
# Coordinate references
signature_append(tuple(self.coordref_signatures))
# 1-d coordinates for each axis. Note that self.axis_ids has
# already been sorted.
axis = self.axis
x = [(axis[identity]['ids'],
axis[identity]['units'],
axis[identity]['hasbounds'],
axis[identity]['coordrefs']) for identity in self.axis_ids]
signature_append(tuple(x))
# Whether or not each axis has a dimension coordinate
x = [False if axis[identity]['dim_coord_index'] is None else True
for identity in self.axis_ids]
signature_append(tuple(x))
# N-d auxiliary coordinates
nd_aux = self.nd_aux
x = [(identity,
nd_aux[identity]['units'],
nd_aux[identity]['axes'],
nd_aux[identity]['hasbounds'],
nd_aux[identity]['coordrefs']) for identity in sorted(nd_aux)]
signature_append(tuple(x))
# Cell measures
msr = self.msr
x = [(units,
msr[units]['axes']) for units in sorted(msr)]
signature_append(tuple(x))
self.signature = tuple(signature)
#--- End: def
def find_coordrefs(self, key, coord):
'''
:Parameters:
key : str
The domain identifier of the coordinate object
coord : cf.Coordinate
:Returns:
out : tuple or None
:Examples:
>>> dim_coord
<CF DimensionCoordinate: ....>
>>> m.find_coordrefs('dim0', dim_coord)
>>> aux_coord
<CF AuxiliaryCoordinate: ....>
>>> m.find_coordrefs('aux1', aux_coord)
'''
coordrefs = self.coordrefs
if not coordrefs:
return None
# Select the coordinate references which contain a pointer to
# this coordinate
names = [ref.name for ref in coordrefs if key in ref.coords]
if not names:
return None
return tuple(sorted(names))
#--- End: def
#--- End: class
[docs]def aggregate(fields,
info=0,
relaxed_units=False,
no_overlap=False,
contiguous=False,
relaxed_identities=False,
ncvar_identities=False,
respect_valid=False,
equal_all=False,
exist_all=False,
equal=None,
exist=None,
ignore=None,
exclude=False,
dimension=(),
concatenate=True,
copy=True,
axes=None,
donotchecknonaggregatingaxes=False,
allow_no_identity=False,
shared_nc_domain=False
):
'''
Aggregate fields into as few fields as possible.
The aggregation of fields may be thought of as the combination fields
into each other to create a new field that occupies a larger domain.
Using the CF aggregation rules, input fields are separated into
aggregatable groups and each group (which may contain just one field)
is then aggregated to a single field. These aggregated fields are
returned in a field list.
**Identities**
In order for aggregation to be possible, fields and their components
need to be unambiguously identifiable. By default, these identities
are taken from `!standard_name` CF properties or else `!id`
attributes. If both of these identifiers are absent then `!long_name`
CF properties or else `!ncvar` attributes may be used if the
*relaxed_identities* parameter is True.
:Parameters:
fields : (sequence of) cf.Field or cf.FieldList
The field or fields to aggregated.
info : int, optional
Print information about the aggregation process. If *info* is
0 then no information is displayed. If *info* is 1 or more
then display information on which fields are unaggregatable,
and why. If *info* is 2 or more then display the structural
signatures of the fields and, when there is more than one
field with the same structural signature, their canonical
first and last coordinate values. If *info* is 3 or more then
display the fields' complete aggregation metadata. By default
*info* is 0 and no information is displayed.
no_overlap : bool, optional
If True then require that aggregated fields have adjacent
dimension coordinate object cells which do not overlap (but
they may share common boundary values). Ignored if the
dimension coordinates objects do not have bounds. See the
*contiguous* parameter.
contiguous : bool, optional
If True then require that aggregated fields have adjacent
dimension coordinate object cells which partially overlap or
share common boundary values. Ignored if the dimension
coordinate objects do not have bounds. See the *no_overlap*
parameter.
relaxed_units : bool, optional
If True then assume that fields or domain items (such as
coordinate objects) with the same identity (as returned by
their `!identity` methods) but missing units all have
equivalent but unspecified units, so that aggregation may
occur. By default such fields are not aggregatable.
allow_no_identity : bool, optional
If True then treat fields with data arrays but with no
identities (see the above notes) as having equal but
unspecified identities, so that aggregation may occur. By
default such fields are not aggregatable.
relaxed_identities : bool, optional
If True then allow fields and their components to be
identified by their `!long_name` CF properties or else
`!ncvar` attributes if their `!standard_name` CF properties or
`!id` attributes are missing.
ncvar_identities : bool, optional
If True then Force fields and their components (such as
coordinates) to be identified by their netCDF file variable
names.
shared_nc_domain : bool, optional
If True then match axes between a field and its contained
ancillary variable and coordinate reference fields via their
netCDF dimension names and not via their domains.
equal_all : bool, optional
If True then require that aggregated fields have the same set
of non-standard CF properties (including
`~cf.Field.long_name`), with the same values. See the
*concatenate* parameter.
equal_ignore : (sequence of) str, optional
Specify CF properties to omit from any properties specified by
or implied by the *equal_all* and *equal* parameters.
equal : (sequence of) str, optional
Specify CF properties for which it is required that aggregated
fields all contain the properties, with the same values. See
the *concatenate* parameter.
exist_all : bool, optional
If True then require that aggregated fields have the same set
of non-standard CF properties (including, in this case,
long_name), but not requiring the values to be the same. See
the *concatenate* parameter.
exist_ignore : (sequence of) str, optional
Specify CF properties to omit from the properties specified by
or implied by the *exist_all* and *exist* parameters.
exist : (sequence of) str, optional
Specify CF properties for which it is required that aggregated
fields all contain the properties, but not requiring the
values to be the same. See the *concatenate* parameter.
respect_valid : bool, optional
If True then the CF properties `~cf.Field.valid_min`,
`~cf.Field.valid_max` and `~cf.Field.valid_range` are taken
into account during aggregation. By default these CF
properties are ignored and are not set in the output fields.
dimension : (sequence of) str, optional
Create new axes for each input field which has one or more of
the given properties. For each CF property name specified, if
an input field has the property then, prior to aggregation, a
new axis is created with an auxiliary coordinate whose datum
is the property's value and the property itself is deleted
from that field.
concatenate : bool, optional
If False then a CF property is omitted from an aggregated
field if the property has unequal values across constituent
fields or is missing from at least one constituent field. By
default a CF property in an aggregated field is the
concatenated collection of the distinct values from the
constituent fields, delimited with the string
``' :AGGREGATED: '``.
copy : bool, optional
If False then do not copy fields prior to aggregation.
Setting this option to False may change input fields in place,
and the output fields may not be independent of the
inputs. However, if it is known that the input fields are
never to accessed again (such as in this case: ``f =
cf.aggregate(f)``) then setting *copy* to False can reduce the
time taken for aggregation.
axes : (sequence of) str, optional
Select axes to aggregate over. Aggregation will only occur
over as large a subset as possible of these axes. Each axis is
identified by the exact identity of a one dimensional
coordinate object, as returned by its `!identity`
method. Aggregations over more than one axis will occur in the
order given. By default, aggregation will be over as many axes
as possible.
donotchecknonaggregatingaxes : bool, optional
If True, and *axes* is set, then checks for consistent data
array values will only be made for one dimensional coordinate
objects which span the any of the given aggregating axes. This
can reduce the time taken for aggregation, but if any those
checks would have failed then this clearly allows the
possibility of an incorrect result. Therefore, this option
should only be used in cases for which it is known that the
non-aggregating axes are in fact already entirely consistent.
:Returns:
out : cf.FieldList
The aggregated fields.
:Examples:
The following six fields comprise eastward wind at two different times
and for three different atmospheric heights for each time:
>>> f
[<CF Field: eastward_wind(latitude(73), longitude(96)>,
<CF Field: eastward_wind(latitude(73), longitude(96)>,
<CF Field: eastward_wind(latitude(73), longitude(96)>,
<CF Field: eastward_wind(latitude(73), longitude(96)>,
<CF Field: eastward_wind(latitude(73), longitude(96)>,
<CF Field: eastward_wind(latitude(73), longitude(96)>]
>>> g = cf.aggregate(f)
>>> g
[<CF Field: eastward_wind(height(3), time(2), latitude(73), longitude(96)>]
>>> g[0].source
'Model A'
>>> g = cf.aggregate(f, dimension=('source',))
[<CF Field: eastward_wind(source(1), height(3), time(2), latitude(73), longitude(96)>]
>>> g[0].source
AttributeError: 'Field' object has no attribute 'source'
'''
# Initialise the cache for coordinate and cell measure hashes,
# first and last values and first and last cell bounds
hfl_cache = _HFLCache()
output_fields = FieldList()
output_fields_append = output_fields.append
if exclude:
exclude = ' NOT'
else:
exclude = ''
atol = ATOL()
rtol = RTOL()
if axes is not None and isinstance(axes, basestring):
axes = (axes,)
# Parse parameters
strict_identities = not (relaxed_identities or ncvar_identities)
if exist_all and equal_all:
raise AttributeError("asdasdas jnf0____")
if equal or exist or ignore:
properties = {'equal' : equal,
'exist' : exist,
'ignore': ignore}
for key, value in properties.iteritems():
if not value:
continue
if isinstance(equal, basestring):
# If it is a string then convert to a single element
# sequence
properties[key] = (value,)
else:
try:
value[0]
except TypeError:
raise TypeError("Bad type of %r parameter: %r" %
(key, type(value)))
#--- End: for
equal = properties['equal']
exist = properties['exist']
ignore = properties['ignore']
if equal and exist:
if set(equal).intersection(exist):
raise AttributeError("888888888888888 asdasdas jnf0____")
if ignore:
ignore = _signature_properties.union(ignore)
else:
ignore = _signature_properties
#--- End: if
unaggregatable = False
status = 0
# ================================================================
# Group together fields with the same structural signature
# ================================================================
signatures = {}
for f in flat(fields):
# ------------------------------------------------------------
# Create the metadata summary, including the structural
# signature
# ------------------------------------------------------------
meta = _Meta(f,
info=info, rtol=rtol, atol=atol,
relaxed_units=relaxed_units,
allow_no_identity=allow_no_identity,
equal_all=equal_all,
exist_all=exist_all,
equal=equal,
exist=exist,
ignore=ignore,
dimension=dimension,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
respect_valid=respect_valid)
if not meta:
unaggregatable = True
status = 1
if info:
print(
"Unaggregatable %r field has%s been output: %s" %
(f.name(''), exclude, meta.message))
if not exclude:
# This field does not have a structural signature, so
# it can't be aggregated. Put it straight into the
# output list and move on to the next input field.
if not copy:
output_fields_append(f)
else:
output_fields_append(f.copy())
#--- End: if
continue
#--- End: if
# ------------------------------------------------------------
# This field has a structural signature, so append it to the
# list of fields with the same structural signature.
# ------------------------------------------------------------
signatures.setdefault(meta.signature, []).append(meta)
#--- End: for
# ================================================================
# Within each group of fields with the same structural signature,
# aggregate as many fields as possible. Sort the signatures so
# that independent aggregations of the same set of input fields
# return fields in the same order.
# ================================================================
for signature in sorted(signatures):
meta = signatures[signature]
if info >= 2:
# Print useful information
meta[0].print_info(info)
print ''
#--- End: if
if len(meta) == 1:
# --------------------------------------------------------
# There's only one field with this signature, so we can
# add it straight to the output list and move on to the
# next signature.
# --------------------------------------------------------
if not copy:
output_fields_append(meta[0].field)
else:
output_fields_append(meta[0].field.copy())
# if info >= 2:
# meta[0].print_info(info)
# if info:
# print(
#"%r field can't be aggregated due to a unique structural signature" %
#meta[0].field.name(''))
continue
#--- End: if
# ------------------------------------------------------------
# Still here? Then there are 2 or more fields with this
# signature which may be aggregatable. These fields need to be
# passed through until no more aggregations are possible. With
# each pass, the number of fields in the group will reduce by
# one for each aggregation that occurs. Each pass represents
# an aggregation in another axis.
# ------------------------------------------------------------
# ------------------------------------------------------------
# For each axis's 1-d coordinates, create the canonical hash
# value and the first and last cell values.
# ------------------------------------------------------------
if axes is None:
# Aggregation will be over as many axes as possible
aggregating_axes = meta[0].axis_ids
_create_hash_and_first_values(meta, None, False, hfl_cache)
#def _create_hash_and_first_values(meta, axes, donotchecknonaggregatingaxes,
# hfl_cache):
else:
# Specific aggregation axes have been selected
aggregating_axes = []
axis_items = meta[0].axis.items()
for axis in axes:
coord = meta[0].field.coord(axis, exact=True)
if coord is None:
continue
coord_identity = coord.name(identity=strict_identities,
ncvar=ncvar_identities)
for identity, value in axis_items:
if (identity not in aggregating_axes and
coord_identity in value['ids']):
aggregating_axes.append(identity)
break
#--- End: for
_create_hash_and_first_values(meta, aggregating_axes,
donotchecknonaggregatingaxes,
hfl_cache)
#--- End: if
if info >= 2:
# Print useful information
for m in meta:
m.print_info(info, signature=False)
print ''
#--- End: if
# Take a shallow copy in case we abandon and want to output
# the original, unaggregated fields.
meta0 = meta[:]
unaggregatable = False
for axis in aggregating_axes:
number_of_fields = len(meta)
if number_of_fields == 1:
break
# --------------------------------------------------------
# Separate the fields with the same structural signature
# into groups such that either within each group the
# fields' domains differ only long the axis or each group
# contains only one field.
#
# Note that the 'a_identity' attribute is set in the
# _group_fields function.
# --------------------------------------------------------
grouped_meta = _group_fields(meta, axis)
if not grouped_meta:
if info:
print(
"Unaggregatable %r fields have%s been output: %s" %
(meta[0].field.name(''), exclude, meta[0].message))
unaggregatable = True
break
#--- End: if
if len(grouped_meta) == number_of_fields:
if info >= 3:
print(
"%r fields can't be aggregated along their %r axis" %
(meta[0].field.name(''), axis))
continue
# --------------------------------------------------------
# Within each group, aggregate as many fields as possible.
# --------------------------------------------------------
for m in grouped_meta:
if len(m) == 1:
continue
# ----------------------------------------------------
# Still here? The sort the fields in place by the
# canonical first values of their 1-d coordinates for
# the aggregating axis.
# ----------------------------------------------------
_sorted_by_first_values(m, axis)
# ----------------------------------------------------
# Check that the aggregating axis's 1-d coordinates
# don't overlap, and don't aggregate anything in this
# group if any do.
# ----------------------------------------------------
if not _ok_coordinate_arrays(m, axis, no_overlap, contiguous,
info):
if info:
print(
"Unaggregatable %r fields have%s been output: %s" %
(m[0].field.name(''), exclude, m[0].message))
unaggregatable = True
break
#--- End: if
# ----------------------------------------------------
# Still here? Then pass through the fields
# ----------------------------------------------------
m0 = m[0].copy()
for m1 in m[1:]:
m0 = _aggregate_2_fields(m0, m1,
rtol=rtol, atol=atol,
respect_valid=respect_valid,
contiguous=contiguous,
no_overlap=no_overlap,
relaxed_units=relaxed_units,
info=info,
concatenate=concatenate,
copy=(copy or not exclude),
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
shared_nc_domain=shared_nc_domain)
if not m0:
# Couldn't aggregate these two fields, so
# abandon all aggregations on the fields with
# this structural signature, including those
# already done.
if info:
print(
"Unaggregatable %r fields have%s been output: %s" %
(m1.field.name(''), exclude, m1.message))
unaggregatable = True
break
#--- End: while
m[:] = [m0]
#--- End: for
if unaggregatable:
break
# --------------------------------------------------------
# Still here? Then the aggregation along this axis was
# completely successful for each sub-group, so reassemble
# the aggregated fields as a single list ready for
# aggregation along the next axis.
# --------------------------------------------------------
meta = [m for gm in grouped_meta for m in gm]
#--- End: for
# Add fields to the output list
if unaggregatable:
# info > 0:
# print ''
status = 1
if not exclude:
if copy:
output_fields.extend((m.field.copy() for m in meta0))
else:
output_fields.extend((m.field for m in meta0))
else:
output_fields.extend((m.field for m in meta))
#--- End: for
aggregate.status = status
if status and info > 0:
print ''
if len(output_fields) == 1:
return output_fields[0]
else:
return output_fields
#--- End: def
# --------------------------------------------------------------------
# Initialise the status
# --------------------------------------------------------------------
aggregate.status = 0
def _create_hash_and_first_values(meta, axes, donotchecknonaggregatingaxes,
hfl_cache):
'''
Updates each field's _Meta object.
:Parameters:
meta : list of _Meta
axes : None or list
donotchecknonaggregatingaxes : bool
:Returns:
None
'''
for m in meta:
domain = m.field.domain
domain_dimensions = domain._axes
m_sort_keys = m.sort_keys
m_sort_indices = m.sort_indices
m_hash_values = m.hash_values
m_first_values = m.first_values
m_last_values = m.last_values
m_id_to_axis = m.id_to_axis
# --------------------------------------------------------
# Create a hash value for each metadata array
# --------------------------------------------------------
# --------------------------------------------------------
# 1-d coordinates
# --------------------------------------------------------
for identity in m.axis_ids:
# print 'identity=', identity #dch
if (axes is not None and donotchecknonaggregatingaxes and
identity not in axes):
x = [None] * len(m.axis[identity]['keys'])
m_hash_values[identity] = x
m_first_values[identity] = x[:]
m_last_values[identity] = x[:]
continue
# Still here?
m_axis_identity = m.axis[identity]
axis = m_id_to_axis[identity]
dim_coord = domain.get(axis, None)
# Find the sort indices for this axis ...
if dim_coord is not None:
# ... which has a dimension coordinate
m_sort_keys[axis] = axis
if not domain.direction(axis):
# Axis is decreasing
sort_indices = slice(None, None, -1)
null_sort = False
else:
# Axis is increasing
sort_indices = slice(None)
null_sort = True
else:
# ... which doesn't have a dimension coordinate but
# does have one or more 1-d auxiliary coordinates
aux = m_axis_identity['keys'][0]
sort_indices = numpy_argsort(domain.get(aux).unsafe_array)
m_sort_keys[axis] = aux
null_sort = False
#-- End: if
m_sort_indices[axis] = sort_indices
hash_values = []
first_values = []
last_values = []
for key, canonical_units in izip(m_axis_identity['keys'],
m_axis_identity['units']):
coord = domain.get(key)
# print repr(coord) #dch
# Get the hash of the data array and its first and
# last values
h, f, l = _get_hfl(coord, canonical_units,
sort_indices, null_sort,
True, False, hfl_cache)
# print h, f, l #dch
first_values.append(f)
last_values.append(l)
if coord._hasbounds:
if coord.isdimension:
# Get the hash of the dimension coordinate
# bounds data array and its first and last
# cell values
hb, fb, lb = _get_hfl(coord.bounds, canonical_units,
sort_indices, null_sort,
False, True, hfl_cache)
m.first_bounds[identity] = fb
m.last_bounds[identity] = lb
else:
# Get the hash of the auxiliary coordinate
# bounds data array
hb = _get_hfl(coord.bounds, canonical_units,
sort_indices, null_sort,
False, False, hfl_cache)
#--- End: if
h = (h, hb)
#--- End: if
hash_values.append(h)
## else:
## coord_units = coord.Units
##
## # Change the coordinate data type if required
## if coord.dtype.char not in ('d', 'S'):
## coord = coord.copy(_only_Data=True)
## coord.dtype = _dtype_float
##
## # Change the coordinate's units to the canonical ones
## coord.Units = canonical_units
##
## # Get the coordinate's data array
## if null_sort:
## array = coord.Data.unsafe_array
## else:
## array = coord.Data.array[sort_indices]
##
## hash_value = hash_array(array)
##
## first_values.append(array.item(0)) #[0])
## last_values.append(array.item(-1)) #[-1])
##
## if coord._hasbounds:
## if null_sort:
## array = coord.bounds.Data.unsafe_array
## else:
## array = coord.bounds.Data.array[sort_indices, ...]
##
## hash_value = (hash_value, hash_array(array))
##
## if key[:3] == 'dim': # can do better than this! DCH
## # Record the bounds of the first and last
## # (sorted) cells of a dimension coordinate
## # (don't need to do this for an auxiliary
## # coordinate).
## array0 = array[0, ...].copy()
## array0.sort()
## m.first_bounds[identity] = array0
##
## array0 = array[-1, ...].copy()
## array0.sort()
## m.last_bounds[identity] = array0
## #--- End: if
##
## hash_values.append(hash_value)
##
## # Reinstate the coordinate's original units
## coord.Units = coord_units
#--- End: for
m_hash_values[identity] = hash_values
m_first_values[identity] = first_values
m_last_values[identity] = last_values
#--- End: for
# ------------------------------------------------------------
# N-d auxiliary coordinates
# ------------------------------------------------------------
if donotchecknonaggregatingaxes:
for aux in m.nd_aux.itervalues():
aux['hash_value'] = None
else:
for aux in m.nd_aux.itervalues():
key = aux['key']
canonical_units = aux['units']
coord = domain.get(key)
axes = [m_id_to_axis[identity] for identity in aux['axes']]
domain_axes = domain_dimensions[key]
if axes != domain_axes:
coord = coord.copy(_only_Data=True)
iaxes = [domain_axes.index(axis) for axis in axes]
coord.transpose(iaxes, i=True)
#--- End: if
sort_indices = tuple([m_sort_indices[axis] for axis in axes])
# Get the hash of the data array
h = _get_hfl(coord, canonical_units, sort_indices,
False, False, False, hfl_cache)
if coord._hasbounds:
# Get the hash of the bounds data array
hb = _get_hfl(coord.bounds, canonical_units,
sort_indices,
False, False, False, hfl_cache)
h = (h, hb)
#--- End: if
aux['hash_value'] = h
## else:
## coord_units = coord.Units
##
## # Change the coordinate data type if required
## if coord.dtype.char not in ('d', 'S'):
## coord = coord.copy(_only_Data=True)
## coord.dtype = _dtype_float
## copied = True
## else:
## copied = False
##
## # Change the coordinate's units to the canonical ones
## coord.Units = aux['units'] #canonical_units
##
## # Get the coordinate's data array
## array = coord.Data.array[sort_indices]
##
## hash_value = hash_array(array)
##
## if coord._hasbounds:
## sort_indices.append(Ellipsis)
## array = coord.bounds.Data.array[sort_indices]
## hash_value = (hash_value, hash_array(array))
##
## aux['hash_value'] = hash_value
##
## # Reinstate the coordinate's original units
## coord.Units = coord_units
#--- End: for
#--- End: if
# ------------------------------------------------------------
# Cell measures
# ------------------------------------------------------------
if donotchecknonaggregatingaxes:
for msr in m.msr.itervalues():
msr['hash_values'] = [None] * len(msr['keys'])
else:
for canonical_units, msr in m.msr.iteritems():
hash_values = []
for key, axes in izip(msr['keys'], msr['axes']):
coord = domain.get(key)
axes = [m_id_to_axis[identity] for identity in axes]
domain_axes = domain_dimensions[key]
if axes != domain_axes:
coord = coord.copy(_only_Data=True)
iaxes = [domain_axes.index(axis) for axis in axes]
coord.transpose(iaxes, i=True)
#--- End: if
sort_indices = [m_sort_indices[axis] for axis in axes]
## if qwerty:
# Get the hash of the data array
h = _get_hfl(coord, canonical_units,
tuple(sort_indices),
False, False, False, hfl_cache)
hash_values.append(h)
## else:
## coord_units = coord.Units
##
## # Change the coordinate data type if required
## if coord.dtype.char not in ('d', 'S'):
## coord = coord.copy(_only_Data=True)
## coord.dtype = _dtype_float
## copied = True
## else:
## copied = False
##
## # Change the coordinate's units to the canonical ones
## coord.Units = canonical_units
##
## array = coord.Data.array[tuple(sort_indices)]
##
## hash_values.append(hash_array(array))
##
## # Reinstate the coordinate's original units
## coord.Units = coord_units
#--- End: for
msr['hash_values'] = hash_values
#--- End: for
#--- End: if
# m.calculate_hash_values = set()
m.cell_values = True
#--- End: for
#--- End: def
def _get_hfl(v, canonical_units, sort_indices, null_sort,
first_and_last_values, first_and_last_bounds,
hfl_cache):
'''
Return the hash value, and optionally first and last values (or cell
bounds)
'''
create_hash = True
create_fl = first_and_last_values
create_flb = first_and_last_bounds
key = None
d = v.Data
if d._pmsize == 1:
partition = d.partitions.matrix.item()
if not partition.part:
key = getattr(partition.subarray, 'file_pointer', None)
if key is not None:
hash_value = hfl_cache.hash.get(key, None)
create_hash = hash_value is None
if first_and_last_values:
first, last = hfl_cache.fl.get(key, (None, None))
create_fl = first is None
if first_and_last_bounds:
first, last = hfl_cache.flb.get(key, (None, None))
create_flb = first is None
#--- End: if
if create_hash or create_fl or create_flb:
# Change the data type if required
if d.dtype.char not in ('d', 'S'):
d = d.copy()
d.dtype = _dtype_float
# Change the units to the canonical ones
units = d.Units
d.Units = canonical_units
# Get the data array
if null_sort:
array = d.unsafe_array
else:
array = d.array[sort_indices]
# Reinstate the original units
d.Units = units
if create_hash:
# if v.standard_name=='latitude':
# print repr(array)
# print array.dtype
hash_value = hash_array(array)
hfl_cache.hash[key] = hash_value
if create_fl:
first = array.item(0)
last = array.item(-1)
hfl_cache.fl[key] = (first, last)
if create_flb:
# Record the bounds of the first and last (sorted) cells
first = numpy_sort(array[0, ...])
last = numpy_sort(array[-1, ...])
hfl_cache.flb[key] = (first, last)
#--- End: if
if first_and_last_values or first_and_last_bounds:
return hash_value, first, last
else:
return hash_value
#--- End: def
def _group_fields(meta, axis):
'''
:Parameters:
meta : list of _Meta
axis : str
The name of the axis to group for aggregation.
:Returns:
out : list of cf.FieldList
'''
axes = meta[0].axis_ids
if axes:
if axis in axes:
# Move axis to the end of the axes list
axes = axes[:]
axes.remove(axis)
axes.append(axis)
#--- End: if
sort_by_axis_ids = itemgetter(*axes)
def _hash_values(m):
return sort_by_axis_ids(m.hash_values)
meta.sort(key=_hash_values)
#--- End: if
# Create a new group of potentially aggregatable fields (which
# contains the first field in the sorted list)
m0 = meta[0]
groups_of_fields = [[m0]]
hash0 = m0.hash_values
for m0, m1 in izip(meta[:-1], meta[1:]):
#-------------------------------------------------------------
# Count the number of axes which are different between the two
# fields
# -------------------------------------------------------------
count = 0
hash1 = m1.hash_values
for identity, value in hash0.iteritems():
if value != hash1[identity]:
count += 1
a_identity = identity
#--- End: for
hash0 = hash1
if count == 1:
# --------------------------------------------------------
# Exactly one axis has different 1-d coordinate values
# --------------------------------------------------------
if a_identity != axis:
# But it's not the axis that we're trying currently to
# aggregate over
groups_of_fields.append([m1])
continue
# Still here? Then it is the axis that we're trying
# currently to aggregate over.
ok = True
# Check the N-d auxiliary coordinates
for identity, aux0 in m0.nd_aux.iteritems():
if (a_identity not in aux0['axes'] and
aux0['hash_value'] != m1.nd_aux[identity]['hash_value']):
# This matching pair of N-d auxiliary coordinates
# does not span the aggregating axis and they have
# different data array values
ok = False
break
#--- End: for
if not ok:
groups_of_fields.append([m1])
continue
# Still here? Then check the cell measures
msr0 = m0.msr
for units in msr0:
for axes, hash_value0, hash_value1 in izip(
msr0[units]['axes'],
msr0[units]['hash_values'],
m1.msr[units]['hash_values']):
if a_identity not in axes and hash_value0 != hash_value1:
# There is a matching pair of cell measures
# with these units which does not span the
# aggregating axis and they have different
# data array values
ok = False
break
#--- End: for
if not ok:
groups_of_fields.append([m1])
continue
# Still here? Then set the identity of the aggregating
# axis
m0.a_identity = a_identity
m1.a_identity = a_identity
# Append field1 to this group of potentially aggregatable
# fields
groups_of_fields[-1].append(m1)
elif not count:
# --------------------------------------------------------
# Zero axes have different 1-d coordinate values, so don't
# aggregate anything in this entire group.
# --------------------------------------------------------
meta[0].message = \
"indistinguishable coordinates or other domain information"
return ()
else:
# --------------------------------------------------------
# Two or more axes have different 1-d coordinate values,
# so create a new sub-group of potentially aggregatable
# fields which contains field1.
# --------------------------------------------------------
groups_of_fields.append([m1])
#--- End: if
#--- End: for
return groups_of_fields
#--- End: def
def _sorted_by_first_values(meta, axis):
'''
Sort fields inplace
:Parameters:
meta : list of _Meta
axis : str
:Returns:
None
'''
sort_by_axis_ids = itemgetter(axis)
def _first_values(m):
return sort_by_axis_ids(m.first_values)
#--- End: def
meta.sort(key=_first_values)
#--- End: def
def _ok_coordinate_arrays(meta, axis, no_overlap, contiguous, info):
'''
Return True if the aggregating axis's 1-d coordinates are all
aggregatable.
It is assumed that the input metadata objects have already been sorted
by the canonical first values of their 1-d coordinates.
:Parameters:
meta : list of _Meta
axis : str
Find the canonical identity of the aggregating axis.
no_overlap : bool
See the `aggregate` function for details.
contiguous : bool
See the `aggregate` function for details.
NOT : str
:Returns:
out : bool
:Examples:
>>> if not _ok_coordinate_arrays(meta, True, False)
... print "Don't aggregate"
'''
m = meta[0]
dim_coord_index = m.axis[axis]['dim_coord_index']
if dim_coord_index is not None:
# ------------------------------------------------------------
# The aggregating axis has a dimension coordinate
# ------------------------------------------------------------
# Check for overlapping dimension coordinate cell centres
dim_coord_index0 = dim_coord_index
for m0, m1 in izip(meta[:-1], meta[1:]):
dim_coord_index1 = m1.axis[axis]['dim_coord_index']
if (m0.last_values[axis][dim_coord_index0] >=
m1.first_values[axis][dim_coord_index1]):
# Found overlap
if info:
meta[0].message = \
"%r dimension coordinate values overlap (%s >= %s)" % \
(m.axis[axis]['ids'][dim_coord_index],
m0.last_values[axis][dim_coord_index0],
m1.first_values[axis][dim_coord_index1])
#
#
#"%r fields can't be aggregated due to their %r dimension coordinate values over#lapping (%s >= %s)" %
#(m.field.name(''),
# m.axis[axis]['ids'][dim_coord_index],
# m0.last_values[axis][dim_coord_index0],
# m1.first_values[axis][dim_coord_index1]))
return
dim_coord_index0 = dim_coord_index1
#--- End: for
if axis in m.first_bounds:
# --------------------------------------------------------
# The dimension coordinates have bounds
# --------------------------------------------------------
if no_overlap:
for m0, m1 in izip(meta[:-1], meta[1:]):
if (m1.first_bounds[axis][0] <
m0.last_bounds[axis][1]):
# Do not aggregate anything in this group
# because overlapping has been disallowed and
# the first cell from field1 overlaps with the
# last cell from field0.
if info:
meta[0].message = \
"%r dimension coordinate bounds values overlap (%s < %s)" % \
(m.axis[axis]['ids'][dim_coord_index],
m1.first_bounds[axis][0],
m0.last_bounds[axis][1])
# print(
#"%r fields can't be aggregated due to their %r dimension coordinate bounds valu#es overlapping (%s < %s)" %
#(m.field.name(''),
# m.axis[axis]['ids'][dim_coord_index],
# m1.first_bounds[axis][0],
# m0.last_bounds[axis][1]
# ))
return
#--- End: for
else:
for m0, m1 in izip(meta[:-1], meta[1:]):
m0_last_bounds = m0.last_bounds[axis]
m1_first_bounds = m1.first_bounds[axis]
if (m1_first_bounds[0] <= m0_last_bounds[0] or
m1_first_bounds[1] <= m0_last_bounds[1]):
# Do not aggregate anything in this group
# because, even though overlapping has been
# allowed, the first cell from field1 overlaps
# in an unreasonable way with the last cell
# from field0.
if info:
meta[0].message = \
"%r dimension coordinate bounds values overlap by too much (%s <= %s)" % \
(m.axis[axis]['ids'][dim_coord_index],
m1_first_bounds[0], m0_last_bounds[0],
m1_first_bounds[1], m0_last_bounds[1])
# print(
#"%r fields can't be aggregated due to their %r dimension coordinate bounds valu#es overlapping by too much (%s <= %s)" %
#(m.field.name(''),
# m.axis[axis]['ids'][dim_coord_index],
# m1_first_bounds[0], m0_last_bounds[0],
# m1_first_bounds[1], m0_last_bounds[1]
# ))
return
#--- End: for
#--- End: if
if contiguous:
for m0, m1 in izip(meta[:-1], meta[1:]):
if (m0.last_bounds[axis][1] <
m1.first_bounds[axis][0]):
# Do not aggregate anything in this group
# because contiguous coordinates have been
# specified and the first cell from field1 is
# not contiguous with the last cell from
# field0.
if info:
meta[0].message = \
"%r dimension coordinate cells are not contiguous (%s < %s)" % \
(m.axis[axis]['ids'][dim_coord_index],
m0.last_bounds[axis][1],
m1.first_bounds[axis][0])
# print(
#"%r fields can't be aggregated due to their %r dimension coordinate cells not b#eing contiguous (%s < %s)" %
#(m.field.name(''),
# m.axis[axis]['ids'][dim_coord_index],
# m0.last_bounds[axis][1],
# m1.first_bounds[axis][0]
# ))
return
#--- End: for
#--- End: if
#--- End: if
else:
# ------------------------------------------------------------
# The aggregating axis does not have a dimension coordinate,
# but it does have at least one 1-d auxiliary coordinate.
# ------------------------------------------------------------
# Check for duplicate auxiliary coordinate values
for i, identity in enumerate(meta[0].axis[axis]['ids']):
set_of_1d_aux_coord_values = set()
number_of_1d_aux_coord_values = 0
for m in meta:
aux = m.axis[axis]['keys'][i]
array = m.field.domain.get(aux).array
set_of_1d_aux_coord_values.update(array)
number_of_1d_aux_coord_values += array.size
if len(set_of_1d_aux_coord_values) != number_of_1d_aux_coord_values:
if info:
meta[0].message = \
"no %r dimension coordinates and %r auxiliary coordinates have duplicate values" % \
(identity, identity)
# print(
#"%r fields can't be aggregated due to their %r axes having no dimension coordin#ates and their %r auxiliary coordinates have duplicate values" %
#(m.field.name(''),
# identity,
# identity))
return
#--- End: for
#--- End: for
#--- End: if
# ----------------------------------------------------------------
# Still here? Then the aggregating axis does not overlap between
# any of the fields.
# ----------------------------------------------------------------
return True
#--- End: def
def _aggregate_2_fields(m0, m1,
rtol=None, atol=None,
info=0,
respect_valid=False,
relaxed_units=False,
no_overlap=False,
contiguous=False,
concatenate=True,
copy=True,
relaxed_identities=False,
ncvar_identities=False,
shared_nc_domain=False):
'''
:Parameters:
m0 : _Meta
m1 : _Meta
contiguous : bool, optional
See the `aggregate` function for details.
rtol : float, optional
See the `aggregate` function for details.
atol : float, optional
See the `aggregate` function for details.
info : int, optional
See the `aggregate` function for details.
no_overlap : bool, optional
See the `aggregate` function for details.
relaxed_units : bool, optional
See the `aggregate` function for details.
relaxed_identities : bool, optional
See the `aggregate` function for details.
ncvar_identities : bool, optional
See the `aggregate` function for details.
:Returns:
out : _Meta or bool
'''
# if copy and not m0.aggregated_field:
# m0.field = m0.field.copy()
a_identity = m0.a_identity
# ----------------------------------------------------------------
# Aggregate coordinate references
# ----------------------------------------------------------------
if m0.coordref_signatures:
t = _aggregate_coordrefs(m0, m1,
axis=a_identity,
rtol=rtol, atol=atol,
respect_valid=respect_valid,
relaxed_units=relaxed_units,
no_overlap=no_overlap, info=info,
contiguous=contiguous,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
shared_nc_domain=shared_nc_domain)
if not t:
return
else:
t = None
# ----------------------------------------------------------------
# Aggregate ancillary variables
# ----------------------------------------------------------------
if m0.ancillary_variables:
av = _aggregate_ancillary_variables(m0, m1,
axis=a_identity,
rtol=rtol, atol=atol,
respect_valid=respect_valid,
relaxed_units=relaxed_units,
no_overlap=no_overlap,
info=info,
contiguous=contiguous,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
shared_nc_domain=shared_nc_domain)
if not av:
return
else:
av = None
# Still here?
field0 = m0.field
field1 = m1.field
if copy:
field1 = field1.copy()
domain0 = field0.domain
domain1 = field1.domain
if t:
# ------------------------------------------------------------
# Update coordinate references
# ------------------------------------------------------------
for key, ref in t.iteritems():
domain0.insert_ref(ref, key=key, copy=False, replace=True)
#--- End: if
if av:
# ------------------------------------------------------------
# Update ancillary variables
# ------------------------------------------------------------
field0.ancillary_variables = av
# ----------------------------------------------------------------
# Map the axes of field1 to those of field0
# ----------------------------------------------------------------
dim1_name_map = {}
for identity in m0.axis_ids:
dim1_name_map[m1.id_to_axis[identity]] = m0.id_to_axis[identity]
dim0_name_map = {}
for axis1, axis0 in dim1_name_map.iteritems():
dim0_name_map[axis0] = axis1
# ----------------------------------------------------------------
# In each field, find the identifier of the aggregating axis.
# ----------------------------------------------------------------
adim0 = m0.id_to_axis[a_identity]
adim1 = m1.id_to_axis[a_identity]
# ----------------------------------------------------------------
# Make sure that, along the aggregating axis, field1 runs in the
# same direction as field0
# ----------------------------------------------------------------
direction0 = domain0.direction(adim0)
if domain1.direction(adim1) != direction0:
field1.flip(adim1, i=True)
# ----------------------------------------------------------------
# Find matching pairs of coordinates and cell measures which span
# the aggregating axis
# ----------------------------------------------------------------
# 1-d coordinates
spanning_variables = [(key0, key1, domain0.get(key0), domain1.get(key1))
for key0, key1 in izip(m0.axis[a_identity]['keys'],
m1.axis[a_identity]['keys'])]
hash_values0 = m0.hash_values[a_identity]
hash_values1 = m1.hash_values[a_identity]
for i, (hash0, hash1) in enumerate(izip(hash_values0, hash_values1)):
try:
hash_values0[i].append(hash_values1[i])
except AttributeError:
hash_values0[i] = [hash_values0[i], hash_values1[i]]
#--- End: for
# N-d auxiliary coordinates
for identity in m0.nd_aux:
aux0 = m0.nd_aux[identity]
aux1 = m1.nd_aux[identity]
if a_identity in aux0['axes']:
key0 = aux0['key']
key1 = aux1['key']
spanning_variables.append((key0, key1,
domain0.get(key0),
domain1.get(key1)))
hash_value0 = aux0['hash_value']
hash_value1 = aux1['hash_value']
try:
hash_value0.append(hash_value1)
except AttributeError:
aux0['hash_value'] = [hash_value0, hash_value1]
#--- End: for
# Cell measures
for units in m0.msr:
hash_values0 = m0.msr[units]['hash_values']
hash_values1 = m1.msr[units]['hash_values']
for i, (axes, key0, key1) in enumerate(izip(m0.msr[units]['axes'],
m0.msr[units]['keys'],
m1.msr[units]['keys'])):
if a_identity in axes:
spanning_variables.append((key0, key1,
domain0.get(key0),
domain1.get(key1)))
try:
hash_values0[i].append(hash_values1[i])
except AttributeError:
hash_values0[i] = [hash_values0[i], hash_values1[i]]
#--- End: for
# ----------------------------------------------------------------
# For each matching pair of coordinates and cell measures which
# span the aggregating axis, insert the one from field1 into the
# one from field0
# ----------------------------------------------------------------
domain_axes = domain0._axes
for key0, key1, item0, item1 in spanning_variables:
item_axes0 = domain0.item_axes(key0)
item_axes1 = domain1.item_axes(key1)
# Ensure that the axis orders are the same in both items
iaxes = [item_axes1.index(dim0_name_map[axis0]) for axis0 in item_axes0]
item1.transpose(iaxes, i=True)
# Find the position of the concatenating axis
axis = item_axes0.index(adim0)
if direction0:
# The fields are increasing along the aggregating axis
item0.Data = Data.concatenate((item0.Data, item1.Data), axis,
_preserve=False)
if item0._hasbounds:
item0.bounds.Data = Data.concatenate((item0.bounds.Data,
item1.bounds.Data),
axis, _preserve=False)
else:
# The fields are decreasing along the aggregating axis
item0.Data = Data.concatenate((item1.Data, item0.Data), axis,
_preserve=False)
if item0._hasbounds:
item0.bounds.Data = Data.concatenate((item1.bounds.Data,
item0.bounds.Data),
axis, _preserve=False)
#--- End: for
# ----------------------------------------------------------------
# Insert the data array from field1 into the data array of field0
# ----------------------------------------------------------------
if m0._hasData:
data_axes0 = domain0.data_axes()
data_axes1 = domain1.data_axes()
# Ensure that both data arrays span the same axes, including
# the aggregating axis.
for axis1 in data_axes1:
axis0 = dim1_name_map[axis1]
if axis0 not in data_axes0:
field0.expand_dims(0, axis0, i=True)
data_axes0.append(axis0)
for axis0 in data_axes0:
axis1 = dim0_name_map[axis0]
if axis1 not in data_axes1:
field1.expand_dims(0, axis1, i=True)
# Find the position of the concatenating axis
if adim0 not in data_axes0:
# Insert the aggregating axis at position 0 because is not
# already spanned by either data arrays
field0.expand_dims(0, adim0, i=True)
field1.expand_dims(0, adim1, i=True)
axis = 0
else:
axis = data_axes0.index(adim0)
# Ensure that the axis orders are the same in both fields
transpose_axes1 = [dim0_name_map[axis0] for axis0 in data_axes0]
if transpose_axes1 != data_axes1:
field1.transpose(transpose_axes1, i=True)
if direction0:
# The fields are increasing along the aggregating axis
field0.Data = Data.concatenate((field0.Data, field1.Data), axis,
_preserve=False)
else:
# The fields are decreasing along the aggregating axis
field0.Data = Data.concatenate((field1.Data, field0.Data), axis,
_preserve=False)
#--- End: if
# Update the size of the aggregating axis in field0
domain0._axes_sizes[adim0] += domain1._axes_sizes[adim1]
# Make sure that field0 has a standard_name, if possible.
if getattr(field0, 'id', None) is not None:
standard_name = field1.getprop('standard_name', None)
if standard_name is not None:
field0.standard_name = standard_name
del field0.id
#--- End: if
#-----------------------------------------------------------------
# Update the properties in field0
#-----------------------------------------------------------------
for prop in set(field0._simple_properties()) | set(field1._simple_properties()):
value0 = field0.getprop(prop, None)
value1 = field1.getprop(prop, None)
if prop in ('valid_min', 'valid_max', 'valid_range'):
if not m0.respect_valid:
try:
field0.delprop(prop)
except AttributeError:
pass
#--- End: if
continue
#--- End: if
if prop == '_FillValue' or prop == 'missing_value':
continue
# Still here?
if equals(value0, value1):
continue
if concatenate:
if value1 is not None:
if value0 is not None:
field0.setprop(prop, '%s :AGGREGATED: %s' % (value0, value1))
else:
field0.setprop(prop, ' :AGGREGATED: %s' % value1)
else:
if value0 is not None:
field0.delprop(prop)
#--- End: for
#-----------------------------------------------------------------
# Update the attributes in field0
#-----------------------------------------------------------------
for attr in m0.attributes | m1.attributes:
value0 = getattr(field0, attr, None)
value1 = getattr(field1, attr, None)
if equals(value0, value1):
continue
if concatenate:
if value1 is not None:
if value0 is not None:
setattr(field0, attr, '%s :AGGREGATED: %s' % (value0, value1))
else:
setattr(field0, attr, ' :AGGREGATED: %s' % value1)
else:
m0.attributes.discard(attr)
if value0 is not None:
delattr(field0, attr)
#--- End: for
# Note that the field in this _Meta object has already been
# aggregated
m0.aggregated_field = True
# ----------------------------------------------------------------
# Return the _Meta object containing the aggregated field
# ----------------------------------------------------------------
return m0
#--- End: def
def _aggregate_coordrefs(m0, m1,
axis=None,
rtol=None, atol=None,
respect_valid=False,
relaxed_units=False,
no_overlap=False,
info=0,
contiguous=False,
relaxed_identities=False,
ncvar_identities=False,
shared_nc_domain=False):
'''
Aggregate fields in coordinate references.
:Parameters:
m0 : _Meta
m1 : _Meta
no_overlap : bool, optional
See the `aggregate` function for details.
contiguous : bool, optional
See the `aggregate` function for details.
rtol : float, optional
See the `aggregate` function for details.
atol : float, optional
See the `aggregate` function for details.
info : int, optional
See the `aggregate` function for details.
relaxed_units : bool, optional
See the `aggregate` function for details.
:Returns:
out : dict
'''
# axis = m0.a_identity
field0 = m0.field
field1 = m1.field
out = {}
for signature in m0.coordref_signatures:
name = signature[0]
key, coordref0 = field0.refs(name, exact=True).popitem()
coordref1 = field1.ref(name, exact=True)
# Initialize the new coordinate reference
new_coordref = CoordinateReference(name=name)
for term in set(coordref0).union(coordref1):
value0 = coordref0.get(term, None)
value1 = coordref1.get(term, None)
if value1 is None and value0 is None:
# ----------------------------------------------------
# Both terms are undefined
# ----------------------------------------------------
continue
if value1 is None:
t, u, m, value = coordref0, coordref1, m0, value0
elif value0 is None:
t, u, m, value = coordref1, coordref0, m1, value1
else:
t = None
if t is not None:
# ----------------------------------------------------
# Exactly one term is undefined
# ----------------------------------------------------
if term in t.coord_terms:
# Term is a coordinate
value = m.field.item(t[term], exact=True)
if value is None:
continue
#--- End: if
default = t.default_value(term)
if default is None:
if info:
m1.message = \
"%r %s %r parameter has no default value" % (name, t.type, term)
return
if isinstance(value, Field):
x = _Meta(value, info=info,
relaxed_units=relaxed_units,
allow_no_identity=True,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
respect_valid=respect_valid)
if not x:
if info:
m1.message = \
"%r %s %r parameter is a field with no structural signature" % \
(name, t.type, term)
return
#--- End: if
if not allclose(value, default, rtol=rtol, atol=atol):
if info:
m1.message = \
"%r %s %r parameters have non-equivalent values" % (name, t.type, term)
return
# Update the new coordinate reference
if term in t.coord_terms:
new_coordref.setcoord(term, t[term])
else:
new_coordref[term] = value
continue
else:
t = coordref0
#--- End: if
coord0 = term in coordref0.coord_terms
coord1 = term in coordref1.coord_terms
if coord0 and coord1:
# ----------------------------------------------------
# Both terms are coordinates
# ----------------------------------------------------
coord0 = field0.item(value0, exact=True)
coord1 = field1.item(value1, exact=True)
coord0_name = coord0.name(identity=m0.strict_identities,
ncvar=m0.ncvar_identities)
coord1_name = coord1.name(identity=m1.strict_identities,
ncvar=m1.ncvar_identities)
if (coord0 is None or
coord1 is None or
coord0_name != coord1_name):
if info:
m1.message = \
"%r %s %r parameters are unaggregatable coordinates" % \
(name, t.type, term)
return
# Update the new coordinate reference
new_coordref.setcoord(term, value0)
continue
if coord0 or coord1:
# ----------------------------------------------------
# Exactly one term is a coordinate
# ----------------------------------------------------
if info:
m1.message = \
"%r %s %r parameters are not all coordinates" % (name, t.type, term)
return
is_field0 = isinstance(value0, Field)
is_field1 = isinstance(value1, Field)
if not is_field0 and not is_field1:
# ----------------------------------------------------
# Neither term is a field
# ----------------------------------------------------
if not allclose(value0, value1, rtol=rtol, atol=atol):
# The values are not equivalent
if info:
m1.message = \
"%r %s %r parameters have non-equivalent values" % (name, t.type, term)
return
# Update the new coordinate reference
new_coordref[term] = value0
continue
if is_field0 != is_field1:
# ----------------------------------------------------
# Exactly one term is a field
# ----------------------------------------------------
if info:
m1.message = \
"%r %s %r parameters are not all fields" % (name, t.type, term)
return
# --------------------------------------------------------
# Both terms are fields
# --------------------------------------------------------
if shared_nc_domain:
role = '%r %s %r' % (name, t.type, term)
value0, message0 = _share_nc_domain(value0, field0, role)
if message0:
if info:
m0.message = message0
return
#--- End: if
value1, message1 = _share_nc_domain(value1, field1, role)
if message1:
if info:
m1.message = message1
return
#--- End: if
#--- End: if
x0 = _Meta(value0, info=info,
relaxed_units=relaxed_units,
allow_no_identity=True,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
respect_valid=respect_valid)
x1 = _Meta(value1, info=info,
relaxed_units=relaxed_units,
allow_no_identity=True,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
respect_valid=respect_valid)
if not (x0 and x1):
# At least one field doesn't have a structual
# signature
if info:
m1.message = \
"%r %s %r parameter is a field with no structural signature" % \
(name, t.type, term)
return
#--- End: if
if axis not in x0.axis and axis not in x1.axis:
# Neither field spans the aggregating axis ...
if value0.equivalent_data(value1, rtol=rtol, atol=atol):
# ... and the fields have equivalent data
# arrays. Therefore we don't need to do any
# aggregation.
# Update the new coordinate reference
new_coordref[term] = value0
continue
else:
# ... and the fields do not have equivalent data
if info:
m1.message = \
"%r %s %r parameters are fields with non-equivalent values" % \
(name, t.type, term)
return
#--- End: if
if not (axis in x0.axis and axis in x1.axis):
# Only one of the fields spans the aggregating axis
if info:
m1.message = \
"%r %s %r parameters are unaggregatable fields" % (name, t.type, term)
# print(
#"%r fields can't be aggregated due to their %r %s %r parameters being fields wh#ich are not aggregatable" %
#(field0.name(''), name, t.type, term))
return
#--- End: if
# Both fields span the aggregating axis, so try to
# aggregate them.
new_value = aggregate((value0, value1),
info=info,
no_overlap=no_overlap,
contiguous=contiguous,
respect_valid=respect_valid,
relaxed_units=relaxed_units,
allow_no_identity=True,
axes=axis,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities)
if len(new_value) == 2:
# Couldn't aggregate them (because we got two fields
# back instead of one)
if info:
m1.message = \
"%r %s %r parameters are unaggregatable fields" % (name, t.type, term)
# print(
#"%r fields can't be aggregated due to their %r %s %r parameters being fields wh#ich are not aggregatable" %
#(field0.name(''), name, t.type, term))
return
#--- End: if
# Successfully aggregated the coordinate reference fields
coordref0[term] = new_value[0] # DCH: Why?????????
# Update the new coordinate reference
new_coordref[term] = new_value[0]
#---End: for
out[key] = new_coordref
#---End: for
return out
#--- End: def
def _aggregate_ancillary_variables(m0, m1,
axis=None,
rtol=None, atol=None,
respect_valid=False,
relaxed_units=False,
no_overlap=False,
info=0,
contiguous=False,
relaxed_identities=False,
ncvar_identities=False,
shared_nc_domain=False):
'''
Aggregate the ancillary variable fields.
:Parameters:
m0 : _Meta
m1 : _Meta
no_overlap : bool, optional
See the `aggregate` function for details.
contiguous : bool, optional
See the `aggregate` function for details.
rtol : float, optional
See the `aggregate` function for details.
atol : float, optional
See the `aggregate` function for details.
info : int, optional
See the `aggregate` function for details.
relaxed_units : bool, optional
See the `aggregate` function for details.
relaxed_identities : bool, optional
See the `aggregate` function for details.
ncvar_identities : bool, optional
See the `aggregate` function for details.
:Returns:
out : cf.FieldList or bool
'''
field0 = m0.field
field1 = m1.field
ancillary_variables = m0.ancillary_variables
# new_ancillary_variables = AncillaryVariables()
new_ancillary_variables = FieldList()
for identity, ancil0 in ancillary_variables.iteritems():
ancil1 = m1.ancillary_variables[identity]
if shared_nc_domain:
ancil0, message0 = _share_nc_domain(ancil0, field0, 'ancillary variable')
if message0:
if info:
m0.message = message0
return
#--- End: if
ancil1, message1 = _share_nc_domain(ancil1, field1, 'ancillary variable')
if message1:
if info:
m1.message = message1
return
#--- End: if
#--- End: if
x0 = _Meta(ancil0, info=info, relaxed_units=relaxed_units,
allow_no_identity=False,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
respect_valid=respect_valid)
x1 = _Meta(ancil1, info=info, relaxed_units=relaxed_units,
allow_no_identity=False,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities,
respect_valid=respect_valid)
if not (x0 and x1):
# At least one field doesn't have a structual signature
if info:
m1.message = \
"%r ancillary variable is a field with no structural signature" % \
ancil0.name('')
return
#--- End: if
if axis not in x0.axis and axis not in x1.axis:
# Neither field spans the aggregating axis ...
if ancil0.equivalent_data(ancil1, rtol=rtol, atol=atol,
traceback=False):
# ... and the fields are equivalent
new_ancillary_variables.append(ancil0)
continue
else:
# ... and the fields are not equivalent
if info:
m1.message = \
"2 %r ancillary variable fields have non-equivalent values" % ancil0.name('')
return
#--- End: if
# Still here?
if not (axis in x0.axis and axis in x1.axis):
if info:
m1.message = \
"3 %r ancillary variable fields are unaggregatable" % ancil0.name('')
return
#--- End: if
# Both fields span the aggregating axis
if (ancil0.axis_size(axis, exact=True) == 1 and
ancil1.axis_size(axis, exact=True) == 1 and
field0.axis_size(axis, exact=True) > 1 or
field1.axis_size(axis, exact=True) > 1):
# The aggregating axis has size 1 in both ancillary fields
# and size > 1 in at least one parent field
if ancil0.equivalent(ancil1, rtol=rtol, atol=atol):
ancillary_variables[identity] = ancil0 ### WHY??
new_ancillary_variables.append(ancil0)
continue
#--- End: if
# Still here? Then try to aggregate the ancillary fields.
new_value = aggregate((ancil0, ancil1), info=info,
no_overlap=no_overlap, contiguous=contiguous,
respect_valid=respect_valid,
relaxed_units=relaxed_units,
allow_no_identity=True,
axes=axis,
relaxed_identities=relaxed_identities,
ncvar_identities=ncvar_identities)
if len(new_value) == 2:
# We got two fields back instead of one, therefore they
# couldn't be aggregated.
if info:
m1.message = \
"4 %r ancillary variable fields are unaggregatable" % ancil0.name('')
return
#--- End: if
# Update the m0.ancillary_variable dictionary, because it
# needs to contain the aggregated field.
ancillary_variables[identity] = new_value[0]
new_ancillary_variables.append(new_value[0])
#---End: for
return new_ancillary_variables
#--- End: def
def _share_nc_domain(child, field, role):
'''
perhaps this should be `cf.Field.share_nc_domain`, as it may be useful
in general.
'''
child_axis_to_ncdim = getattr(child.domain, 'nc_dimensions', {})
if len(set(child_axis_to_ncdim.values())) != len(child_axis_to_ncdim):
message = \
"%s %s field can't share domain with its parent field (ambiguous netCDF dimension names)" % \
(role, child.name(''))
return None, message
#--- End: if
field_axis_to_ncdim = getattr(field.domain, 'nc_dimensions', {})
n_axes = len(field_axis_to_ncdim)
field_ncdim_to_axis = dict([(v, k) for k, v in field_axis_to_ncdim.iteritems()])
if len(field_ncdim_to_axis) != n_axes:
message = \
"%s %s field can't share domain with its parent field (ambiguous netCDF dimension names)" % \
(role, child.name(''))
return None, message
#--- End: if
# Remove all items from the childlary field
new_child = child.copy()
new_child.remove_items()
for c_axis in child.data_axes():
# Find the parent field axis (f_axis) which correposnds to the
# child field axis (c_axis)
c_ncdim = child_axis_to_ncdim.get(c_axis, None)
f_axis = field_ncdim_to_axis.get(c_ncdim, None)
if f_axis is None:
message = \
"%s %s field can't share domain with its parent field (axis has no netCDF dimension name)" % \
(role, child.name(''))
return None, message
#--- End: if
# Copy 1-d dimension and auxiliary coordinates from the parent
# field to the child field
for coord in field.coords(axes=f_axis, ndim=1).itervalues():
if coord.isdimension:
new_child.insert_dim(coord, axis=c_axis)
else:
new_child.insert_aux(coord, axes=(c_axis,))
#--- End: for
#--- End: for
return new_child, None
#--- End: def