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.insert_axis(1) # axis = domain.new_axis_identifier() # print aux_coord, aux_coord.ndim, axis 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': # }, {'name': 'forecast_ref_time', 'key': # 'aux0', '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 '' % (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 >>> m.find_coordrefs('dim0', dim_coord) >>> aux_coord >>> 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 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 [, , , , , ] >>> g = cf.aggregate(f) >>> g [] >>> g[0].source 'Model A' >>> g = cf.aggregate(f, dimension=('source',)) [] >>> 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