Source code for cdms2.avariable

# Automatically adapted for numpy.oldnumeric Aug 01, 2007 by
# Further modified to be pure new numpy June 24th 2008

"CDMS Variable objects, abstract interface"
import numpy
import string
import re
import warnings
from .cdmsobj import CdmsObj
import cdms2
from .slabinterface import Slab
from .sliceut import splitSliceExt, splitSlice
from .error import CDMSError
from .axis import axisMatchIndex, axisMatchAxis, axisMatches, unspecified, CdtimeTypes, AbstractAxis
from . import selectors
import copy
from .mvCdmsRegrid import CdmsRegrid, getBoundList, _getCoordList
from regrid2.mvGenericRegrid import guessPeriodicity
# import PropertiedClasses
from .convention import CF1
from .grid import AbstractRectGrid
# import internattr
from six import string_types

InvalidRegion = "Invalid region: "
OutOfRange = "Coordinate interval is out of range or intersection has no data: "
NotImplemented = "Child of AbstractVariable failed to implement: "
# Backward compatibility with numpy behavior
_numeric_compatibility = False
# False: return scalars from 0-D slices
#        MV axis=None by default
# True:  return 0-D arrays
#        MV axis=1 by default


[docs]def getMinHorizontalMask(var): """ Get the minimum mask associated with 'x' and 'y' (i.e. with the min number of ones) across all axes Parameters ---------- var : CDMS variable with a mask N/A : None Returns ------- mask array or None : if order 'x' and 'y' were not found """ from distarray import MultiArrayIter if not hasattr(var, 'mask'): return None shp = var.shape ndims = len(shp) order = var.getOrder() # e.g. 'zxty-', ndims = 5 # run a few checks numX = order.count('x') numY = order.count('y') num_ = order.count('-') hasXY = (numX == 1) and (numY == 1) if numX + numY + num_ < 2: msg = """ Not able to locate the horizontal (y, x) axes for order = %s in getMinHorizontalMask """ % str(order) raise CDMSError(msg) ps = [] # index position of x/y, e.g. [1,3] es = [] # end indices, sizes of x/y axes nonHorizShape = [] found = False for i in range(ndims - 1, -1, -1): # iterate backwards because the horizontal # axes are more likely to be last o = order[i] # curvilinear coordinates have '-' in place of # x or y, also check for '-' but exit if we think # we found the x and y coords if not found and (o in 'xy') or (not hasXY and o == '-'): ps = [i, ] + ps es = [shp[i], ] + es if len(ps) == 2: found = True else: nonHorizShape = [shp[i], ] + nonHorizShape if len(ps) == 2: # found all the horizontal axes, start with mask # set to invalid everywhere mask = numpy.ones(es, numpy.bool8) # iterate over all non-horizontal axes, there can be as # many as you want... for it in MultiArrayIter(nonHorizShape): inds = it.getIndices() # (i0, i1, i2) # build the slice operator, there are three parts # (head, middle, and tail), some parts may be # missing # slce = 'i0,' + ':,' + 'i1,' + ':,' + 'i2,' slce = ('%d,' * ps[0]) % tuple(inds[:ps[0]]) + ':,' \ + ('%d,' * (ps[1] - ps[0] - 1)) % tuple(inds[ps[0]:ps[1] - 1]) \ + ':,' + ('%d,' * (ndims - ps[1] - 1) ) % tuple(inds[ps[1] - 1:]) # evaluate the slice for this time, level.... mask &= eval('var.mask[' + slce + ']') return mask else: msg = """ Could not find all the horizontal axes for order = %s in getMinHorizontalMask """ % str(order) raise CDMSError(msg) return None
[docs]def setNumericCompatibility(mode): global _numeric_compatibility if mode is True or mode == 'on': _numeric_compatibility = True elif mode is False or mode == 'off': _numeric_compatibility = False
[docs]def getNumericCompatibility(): return _numeric_compatibility
[docs]class AbstractVariable(CdmsObj, Slab): """Not to be called by users. Parameters ---------- variableNode is the variable tree node, if any. parent is the containing dataset instance. """
[docs] def info(self, flag=None, device=None): Slab.info(self, flag, device)
[docs] def __init__(self, parent=None, variableNode=None): if variableNode is not None and variableNode.tag != 'variable': raise CDMSError('Node is not a variable node') CdmsObj.__init__(self, variableNode) val = self.__cdms_internals__ + ['id', 'domain'] self.___cdms_internals__ = val Slab.__init__(self) self.id = None # Transient variables key on this to create a default ID self.parent = parent self._grid_ = None # Variable grid, if any if not hasattr(self, 'missing_value'): self.missing_value = None else: if isinstance(self.missing_value, bytes): self.missing_value = None elif isinstance(self.missing_value, string_types): self.missing_value = None elif numpy.isnan(self.missing_value): self.missing_value = None
# Reminder: children to define self.shape and set self.id
[docs] def __array__(self, t=None, context=None): # Numeric, ufuncs call this return numpy.ma.filled(self.getValue(squeeze=0))
[docs] def __call__(self, *args, **kwargs): """ Selection of a subregion using selectors. Parameters ---------- raw : if set to 1, return numpy.ma only squeeze : if set to 1, eliminate any dimension of length 1 grid : if given, result is regridded ont this grid order : if given, result is permuted into this order Returns ------- Subregion selected """ # separate options from selector specs d = kwargs.copy() raw = d.setdefault('raw', 0) squeeze = d.setdefault('squeeze', 0) grid = d.setdefault('grid', None) order = d.setdefault('order', None) del d['squeeze'], d['grid'], d['order'], d['raw'] # make the selector s = selectors.Selector(*args, **d) # get the selection return s.unmodified_select(self, raw=raw, squeeze=squeeze, order=order, grid=grid)
select = __call__
[docs] def rank(self): return len(self.shape)
[docs] def _returnArray(self, ar, squeeze, singles=None): # ar is a Numeric array, numpy.ma, or scalar, possibly numpy.ma.masked. # job is to make sure we return an numpy.ma or a scalar. # If singles is not None, squeeze dimension indices in singles inf = 1.8e308 if isinstance(ar, cdms2.tvariable.TransientVariable): result = numpy.ma.array(ar._data, mask=ar.mask) # already numpy.ma, only need squeeze. elif numpy.ma.isMaskedArray(ar): result = ar elif isinstance(ar, numpy.ndarray): missing = self.getMissing() if missing is None: result = numpy.ma.masked_array(ar) elif missing == inf or missing != missing: # (x!=x) ==> x is NaN result = numpy.ma.masked_object(ar, missing, copy=0) elif ar.dtype.char == 'c' or ar.dtype.char == 'S': # umath.equal is not implemented resultmask = (ar == missing) if not resultmask.any(): resultmask = numpy.ma.nomask result = numpy.ma.masked_array( ar, mask=resultmask, fill_value=missing).astype(str) else: result = numpy.ma.masked_values(ar, missing, copy=0) elif ar is numpy.ma.masked: return ar else: # scalar, but it might be the missing value missing = self.getMissing() if missing is None: return ar # scalar else: result = numpy.ma.masked_values(ar, missing, copy=0) squoze = 0 if squeeze: n = 1 newshape = [] for s in result.shape: if s == 1: squoze = 1 continue else: n = n * s newshape.append(s) elif singles is not None: n = 1 newshape = [] oldshape = result.shape for i in range(len(oldshape)): if i in singles: squoze = 1 continue else: s = oldshape[i] n = n * s newshape.append(s) else: n = numpy.ma.size(result) if n == 1 and squeeze: return numpy.ma.ravel(result)[0] # scalar or masked if squoze: result.shape = newshape return result
[docs] def generateGridkey(self, convention, vardict): """Determine if the variable is gridded. Parameters ---------- convention : Metadata convention class vardict : Variable metedata Returns ------- ((latname, lonname, order, maskname, class), lat, lon) if gridded (None, None, None) if not gridded """ lat, nlat = convention.getVarLatId(self, vardict) lon, nlon = convention.getVarLonId(self, vardict) if (lat is not None) and (lat is lon): raise CDMSError( "Axis %s is both a latitude and longitude axis! Check standard_name and/or axis attributes." % lat.id) # Check for 2D grid if (lat is None) or (lon is None): return None, lat, lon # Check for a rectilinear grid if isinstance(lat, AbstractAxis) and isinstance( lon, AbstractAxis) and (lat.rank() == lon.rank() == 1): return self.generateRectGridkey(lat, lon), lat, lon # Check for a curvilinear grid: if lat.rank() == lon.rank() == 2: # check that they are defined on the same indices as self vardomain = self.getAxisIds() allok = 1 for axisid in lat.getAxisIds(): if axisid not in vardomain: allok = 0 break if allok: for axisid in lon.getAxisIds(): if axisid not in vardomain: allok = 0 break # It's a curvilinear grid if allok: if hasattr(lat, 'maskid'): maskid = lat.maskid else: maskid = '' return (lat.id, lon.id, 'yx', maskid, 'curveGrid'), lat, lon # Check for a generic grid: if lat.rank() == lon.rank() == 1: # check that they are defined on the same indices as self vardomain = self.getAxisIds() allok = 1 for axisid in lat.getAxisIds(): if axisid not in vardomain: allok = 0 break if allok: for axisid in lon.getAxisIds(): if axisid not in vardomain: allok = 0 break # It's a generic grid if allok: if hasattr(lat, 'maskid'): maskid = lat.maskid else: maskid = '' return (lat.id, lon.id, 'yx', maskid, 'genericGrid'), lat, lon return None, lat, lon
[docs] def generateRectGridkey(self, lat, lon): """Determine if the variable is gridded, rectilinear. Parameters ---------- lat : latitude axis lon : longitude axis Returns ------- (latname, lonname, order, maskname, class) if gridded, None if not gridded.""" ilat = ilon = -1 k = 0 for axis in self.getAxisList(): if axis is lon: ilon = k elif axis is lat: ilat = k k += 1 if ilat == -1: raise CDMSError( "Cannot find latitude axis; check standard_name and/or axis attributes") if ilon == -1: raise CDMSError( "Cannot find longitude axis; check standard_name and/or axis attributes") if ilat < ilon: order = "yx" else: order = "xy" gridkey = (lat.id, lon.id, order, '', 'rectGrid') return gridkey
[docs] def isAbstractCoordinate(self): return 0
# Set the variable grid
[docs] def setGrid(self, grid): if grid is None: gridok = 1 else: alist = [d[0] for d in self.getDomain()] gridok = grid.checkAxes(alist) if not gridok: raise CDMSError( "grid does not match axes for variable %s" % self.id) self._grid_ = grid
[docs] def getDomain(self): "Get the list of axes" raise CDMSError("getDomain not overriden in child")
[docs] def getConvention(self): "Get the metadata convention associated with this object." if hasattr(self, 'parent') and self.parent is not None: result = self.parent._convention_ else: result = CF1 return result
# A child class may want to override this
[docs] def getAxis(self, n): """Get the n-th axis. Parameters ---------- n : Axis number Returns ------- if n < 0: n = n + self.rank() self.getDomain()[n][0]""" if n < 0: n = n + self.rank() return self.getDomain()[n][0]
[docs] def getAxisIndex(self, axis_spec): """Get the index of the axis specificed by axis_spec. Parameters ---------- axis_spec : Returns ------- the axis index or -1 if no match is found. """ for i in range(self.rank()): if axisMatches(self.getAxis(i), axis_spec): return i return -1
[docs] def hasCellData(self): """ If any of the variable's axis has explicit bounds, we have cell data otherwise we have point data. Returns ------- True or False if axis has cell data. """ for axis in self.getAxisList(): if (axis.getExplicitBounds() is not None): return True return False
[docs] def getAxisListIndex(self, axes=None, omit=None, order=None): """Get Axis List Index Returns ------- a list of indices of axis objects Notes ----- If axes is **not** `None`, include only certain axes. less the ones specified in omit. If axes is `None`, use all axes of this variable. Other specificiations are as for axisMatchIndex. """ return axisMatchIndex(self.getAxisList(), axes, omit, order)
[docs] def getAxisList(self, axes=None, omit=None, order=None): """Get the list of axis objects Notes If axes is **not** `None`, include only certain axes. If omit is **not** `None`, omit those specified by omit. Arguments omit or axes may be as specified in axisMatchAxis order is an optional string determining the output order """ alist = [d[0] for d in self.getDomain()] return axisMatchAxis(alist, axes, omit, order)
[docs] def getAxisIds(self): """Get a list of axis identifiers. Returns ------- array list of axis ids""" return [x[0].id for x in self.getDomain()]
# Return the grid
[docs] def getGrid(self): return self._grid_
[docs] def getMissing(self, asarray=0): """Get Missing Parameters ---------- asarray : '0' : scalar '1' : numpy array Returns ------- the missing value as a scalar, or as a numpy array if asarray==1""" if hasattr(self, 'missing_value'): try: mv = self.missing_value.item() except BaseException: mv = self.missing_value if mv is None and hasattr(self, '_FillValue'): mv = self._FillValue if asarray == 0 and isinstance(mv, numpy.ndarray): mv = mv[0] if isinstance(mv, string_types) and self.dtype.char not in [ '?', 'c', 'O', 'S']: try: mv = float(mv) except BaseException: if hasattr(self, '_FillValue'): try: mv = float(self._FillValue) except BaseException: mv = None else: mv = None return mv
[docs] def _setmissing(self, name, value): self.setMissing(value)
[docs] def setMissing(self, value): """Set the missing value. Parameters ---------- value : scalar, a single-valued numpy array, or None. Note : The value is cast to the same type as the variable.""" # Check for None first, so that constructors can # set missing_value before typecode() is initialized. if value is None: self._basic_set('missing_value', value) return selftype = self.typecode() valuetype = type(value) if valuetype is numpy.ndarray: value = value.astype(selftype).item() elif isinstance(value, numpy.floating) or isinstance(value, numpy.integer): value = numpy.array([value], selftype) elif valuetype in [float, int, int, complex]: try: value = numpy.array([value], selftype) # Set fill value when ar[i:j] returns a masked value except BaseException: value = numpy.array( [numpy.ma.default_fill_value(self)], selftype) # '?' for Boolean and object elif isinstance(value, (str, numpy.string_, numpy.str, numpy.string0, numpy.str_)) and selftype in ['?', 'c', 'O', 'S']: pass else: raise CDMSError('Invalid missing value %s' % repr(value)) self.missing_value = value
[docs] def getTime(self): """Get the first time dimension. Returns ------- First Time dimension axis or `None`. """ for k in range(self.rank()): axis = self.getAxis(k) if axis.isTime(): return axis break else: return None
[docs] def getForecastTime(self): """Get the first forecast time dimension. Returns ------- First forecast time dimension axis or `None`. """ for k in range(self.rank()): axis = self.getAxis(k) if axis.isForecast(): return axis break else: return None
[docs] def getForecast(self): return self.getForecastTime()
[docs] def getLevel(self): """Get the first vertical level dimension in the domain. Returns ------- First vertical level dimension axis or `None`. """ for k in range(self.rank()): axis = self.getAxis(k) if axis.isLevel(): return axis break else: return None
[docs] def getLatitude(self): """Get the first latitude dimension. Returns ------- First latitude dimension axis or `None`. """ grid = self.getGrid() if grid is not None: result = grid.getLatitude() else: result = None if result is None: for k in range(self.rank()): result = self.getAxis(k) if result.isLatitude(): break else: result = None return result
[docs] def getLongitude(self): """Get the first longitude dimension. Returns ------- First longitude dimension axis or `None`. """ grid = self.getGrid() if grid is not None: result = grid.getLongitude() else: result = None if result is None: for k in range(self.rank()): result = self.getAxis(k) if result.isLongitude(): break else: result = None return result
# Get an order string, such as "tzyx"
[docs] def getOrder(self, ids=0): """ Get Order Parameters ---------- id : 0 or 1 Returns ------- the order string, such as t, z, y, x (time, level, lat, lon). Notes * if ids == 0 (the default) for an axis that is not t,z,x,y the order string will contain a (-) character in that location. The result string will be of the same length as the number of axes. This makes it easy to loop over the dimensions. * if ids == 1 those axes will be represented in the order string as (id) where id is that axis' id. The result will be suitable for passing to order2index to get the corresponding axes, and to orderparse for dividing up into components. """ order = "" for k in range(self.rank()): axis = self.getAxis(k) if axis.isLatitude(): order = order + "y" elif axis.isLongitude(): order = order + "x" elif axis.isLevel(): order = order + "z" elif axis.isTime(): order = order + "t" elif ids: order = order + '(' + axis.id + ')' else: order = order + "-" return order
[docs] def subSlice(self, *specs, **keys): speclist = self._process_specs(specs, keys) numericSqueeze = keys.get('numericSqueeze', 0) # Get a list of single-index specs if numericSqueeze: singles = self._single_specs(specs) else: singles = None slicelist = self.specs2slices(speclist, force=1) d = self.expertSlice(slicelist) squeeze = keys.get('squeeze', 0) raw = keys.get('raw', 0) order = keys.get('order', None) grid = keys.get('grid', None) # Force result to have these axes forceaxes = keys.get('forceaxes', None) raweasy = raw == 1 and order is None and grid is None if not raweasy: if forceaxes is None: axes = [] allaxes = [None] * self.rank() for i in range(self.rank()): slice = slicelist[i] if squeeze and numpy.ma.size(d, i) == 1: continue elif numericSqueeze and i in singles: continue # Don't wrap square-bracket slices axis = self.getAxis(i).subaxis( slice.start, slice.stop, slice.step, wrap=( numericSqueeze == 0)) axes.append(axis) allaxes[i] = axis else: axes = forceaxes # Slice the grid, if non-rectilinear. Don't carry rectilinear grids, since # they can be inferred from the domain. selfgrid = self.getGrid() if selfgrid is None or isinstance(selfgrid, AbstractRectGrid): resultgrid = None else: alist = [item[0] for item in self.getDomain()] gridslices, newaxes = selfgrid.getGridSlices( alist, allaxes, slicelist) # If one of the grid axes was squeezed, the result grid is None if None in newaxes: resultgrid = None else: resultgrid = selfgrid.subSlice( *gridslices, **{'forceaxes': newaxes}) resultArray = self._returnArray(d, squeeze, singles=singles) if self.isEncoded(): resultArray = self.decode(resultArray) newmissing = resultArray.fill_value else: newmissing = self.getMissing() if raweasy: return resultArray elif len(axes) > 0: # If forcing use of input axes, make sure they are not copied. # Same if the grid is not rectilinear - this is when forceaxes is # set. copyaxes = (forceaxes is None) and (resultgrid is None) result = TransientVariable(resultArray, copy=0, fill_value=newmissing, axes=axes, copyaxes=copyaxes, grid=resultgrid, attributes=self.attributes, id=self.id) if grid is not None: order2 = grid.getOrder() if order is None: order = order2 elif order != order2: raise CDMSError('grid, order options not compatible.') result = result.reorder(order).regrid(grid) if raw == 0: return result else: return result.getSlice(squeeze=0, raw=1) else: # Return numpy.ma for zero rank, so that __cmp__ works. return resultArray
[docs] def getSlice(self, *specs, **keys): """getSlice takes arguments of the following forms and produces a return array. Parameters ---------- raw : if set to 1, return numpy.ma only squeeze : if set to 1, eliminate any dimension of length 1 grid : if given, result is regridded ont this grid. order : if given, result is permuted into this order numericSqueeze : if index slice is given, eliminate that dimension. isitem : if given, result is return as a scaler for 0-D data Notes There can be zero or more positional arguments, each of the form: #. a single integer n, meaning slice(n, n+1) #. an instance of the slice class #. a tuple, which will be used as arguments to create a slice #. `None` or `:`, which means a slice covering that entire dimension #. Ellipsis (...), which means to fill the slice list with `:` leaving only enough room at the end for the remaining positional arguments There can be keyword arguments of the form key = value, where key can be one of the names `time`, `level`, `latitude`, or `longitude`. The corresponding value can be any of (1)-(5) above. There must be no conflict between the positional arguments and the keywords. In (1)-(5) negative numbers are treated as offsets from the end of that dimension, as in normal Python indexing. """ # Turn on squeeze and raw options by default. keys['numericSqueeze'] = keys.get('numericSqueeze', 0) keys['squeeze'] = keys.get('squeeze', 1 - keys['numericSqueeze']) keys['raw'] = keys.get('raw', 1) keys['order'] = keys.get('order', None) keys['grid'] = keys.get('grid', None) isitem = keys.get('isitem', 0) result = self.subSlice(*specs, **keys) # return a scalar for 0-D slices if isitem and result.size == 1 and ( not _numeric_compatibility) and not result.mask.item(): result = result.item() return result
[docs] def expertSlice(self, slicelist): raise CDMSError(NotImplemented + 'expertSlice')
[docs] def getRegion(self, *specs, **keys): """ Read a region of data. A region is an n-dimensional rectangular region specified in coordinate space. Parameters ---------- slice : is an argument list, each item of which has one of the following forms: * x, where x is a scalar * Map the scalar to the index of the closest coordinate value. * (x, y) * Map the half-open coordinate interval [x,y) to index interval. * (x, y, 'cc') * Map the closed interval [x,y] to index interval. Other options are 'oo' (open), 'oc' (open on the left), and 'co' (open on the right, the default). * (x, y, 'co', cycle) * Map the coordinate interval with wraparound. If no cycle is specified, wraparound will occur iff axis.isCircular() is true. Ellipsis : Represents the full range of all dimensions bracketed by non-Ellipsis items. None, colon : Represents the full range of one dimension. Notes Only one dimension may be wrapped. Example Suppose the variable domain is `(time, level, lat, lon)`. Then >>> getRegion((10, 20), 850, Ellipsis,(-180, 180)) retrieves : * all times t such that 10.<=t<20. * level 850. * all values of all dimensions between level and lon (namely, lat). * longitudes x such that -180<=x<180. This will be wrapped unless lon.topology=='linear'. """ # By default, squeeze and raw options are on keys['squeeze'] = keys.get('squeeze', 1) keys['raw'] = keys.get('raw', 1) keys['order'] = keys.get('order', None) keys['grid'] = keys.get('grid', None) return self.subRegion(*specs, **keys)
[docs] def subRegion(self, *specs, **keys): speclist = self._process_specs(specs, keys) slicelist = self.reg_specs2slices(speclist) squeeze = keys.get('squeeze', 0) raw = keys.get('raw', 0) order = keys.get('order', None) grid = keys.get('grid', None) raweasy = raw == 1 and order is None and grid is None if grid is not None and order is None: order = grid.getOrder() # Check if any slice wraps around. wrapdim = -1 axes = [] circulardim = None for idim in range(len(slicelist)): item = slicelist[idim] axis = self.getAxis(idim) axislen = len(axis) if(axis.isCircular()): circulardim = idim wraptest1 = (axis.isCircular() and speclist[idim] != unspecified) start, stop = item.start, item.stop wraptest2 = not ( (start is None or ( 0 <= start < axislen)) and ( stop is None or ( 0 <= stop <= axislen))) if (wraptest1 and wraptest2): if wrapdim >= 0: raise CDMSError("Too many dimensions wrap around.") wrapdim = idim break else: # No wraparound, just read the data # redo the speclist -> slice if passed circular test but not # wrapped test if(circulardim is not None): slicelist = self.reg_specs2slices(speclist, force=circulardim) d = {'raw': raw, 'squeeze': squeeze, 'order': order, 'grid': grid, } return self.subSlice(*slicelist, **d) # # get the general wrap slice (indices that are neg -> pos and vica versa) # wrapslice = slicelist[wrapdim] wrapaxis = self.getAxis(wrapdim) length = len(wrapaxis) # # shift the wrap slice to the positive side and calc number of cycles shifted # wb = wrapslice.start we = wrapslice.stop ws = wrapslice.step size = length cycle = self.getAxis(wrapdim).getModulo() # # ncycle: # number of cycles for slicing purposes and # resetting world coordinate from the positive only direction # # ncyclesrev: # resetting the world coordinate for reversed direction # ncycles = 0 ncyclesrev = 0 if(ws > 0): if(wb > 0): ncycles = 1 while(wb >= 0): wb = wb - size we = we - size ncycles = ncycles - 1 else: ncycles = 0 while(wb < 0): wb = wb + size we = we + size ncycles = ncycles + 1 if(wb < 0): wb = wb + size we = we + size # reversed direction else: # do the ncycles for resetting world coordinate wbrev = wb werev = we werevNoneTest = 0 if(werev is None): werev = 0 werevNoneTest = 1 ncycleRevStart = 1 if(wbrev > 0): ncyclesrev = ncycleRevStart while(wbrev >= 0): wbrev = wbrev - size werev = werev - size ncyclesrev = ncyclesrev - 1 else: ncyclesrev = 0 while(wbrev < 0): wbrev = wbrev + size werev = werev + size ncyclesrev = ncyclesrev + 1 while(werev < 0): wbrev = wbrev + size werev = werev + size # number of cycles to make the slice positive while(we < 0 and we is not None): wb = wb + size we = we + size ncycles = ncycles + 1 wb = wbrev we = werev if(werevNoneTest): we = None wrapslice = slice(wb, we, ws) # # calc the actual positive slices and create data array # donew = 1 if(donew): wraps = splitSliceExt(wrapslice, length) for kk in range(0, len(wraps)): sl = wraps[kk] slicelist[wrapdim] = sl if(kk == 0): ar1 = self.getSlice(squeeze=0, *slicelist) result = ar1 else: ar2 = self.getSlice(squeeze=0, *slicelist) result = numpy.ma.concatenate((result, ar2), axis=wrapdim) else: wrap1, wrap2 = splitSlice(wrapslice, length) slicelist[wrapdim] = wrap1 ar1 = self.getSlice(squeeze=0, *slicelist) slicelist[wrapdim] = wrap2 ar2 = self.getSlice(squeeze=0, *slicelist) result = numpy.ma.concatenate((ar1, ar2), axis=wrapdim) if raweasy: return self._returnArray(result, squeeze) # ---------------------------------------------------------------------- # # set ALL the axes (transient variable) # # ---------------------------------------------------------------------- axes = [] for i in range(self.rank()): if squeeze and numpy.ma.size(result, i) == 1: continue sl = slicelist[i] if i == wrapdim: axis = self.getAxis(i).subAxis(wb, we, ws) if(ws > 0): delta_beg_wrap_dimvalue = ncycles * cycle else: delta_beg_wrap_dimvalue = ncyclesrev * cycle isGeneric = [False] b = axis.getBounds(isGeneric) - delta_beg_wrap_dimvalue axis.setBounds(b, isGeneric=isGeneric[0]) axis[:] = ( axis[:] - delta_beg_wrap_dimvalue).astype( axis.typecode()) else: axis = self.getAxis(i).subaxis(sl.start, sl.stop, sl.step) axes.append(axis) result = self._returnArray(result, squeeze) result = TransientVariable(result, copy=0, fill_value=self.missing_value, axes=axes, attributes=self.attributes, id=self.id) if grid is not None: order2 = grid.getOrder() if order is None: order = order2 elif order != order2: raise CDMSError('grid, order options not compatible.') result = result.reorder(order).regrid(grid) if raw == 0: return result else: return result.getSlice(squeeze=0, raw=1)
[docs] def getValue(self, squeeze=1): """Get the entire set of values. Returns ------- All values and elimite the 1-D dimension. """ return self.getSlice(Ellipsis, squeeze=squeeze)
[docs] def assignValue(self, data): raise CDMSError(NotImplemented + 'assignValue')
[docs] def reorder(self, order): """Reorder per the specification order. Parameters ---------- order : string can be "tzyx" with all possible axes permutation. Returns ------- New reordered variable. """ if order is None: return self axes = self.getAxisList() permutation = order2index(axes, order) if permutation == list(range(len(axes))): return self return MV.transpose(self, permutation)
[docs] def regrid(self, togrid, missing=None, order=None, mask=None, **keywords): """return self regridded to the new grid. One can use the regrid2.Regridder optional arguments as well. Example >>> new_cdmsVar = cdmsVar.regrid(newGrid) # uses libcf >>> new_cdmsVar = cdmsVar.regrid(newGrid, regridMethod = 'conserve', coordSys = 'cart') Parameters ---------- togrid : togrid destination grid. CDMS grid missing : Optional missing missing values order : Optional order axis order mask : Optional mask grid/data mask **keyords keywords optional keyword arguments dependent on regridTool Returns ------- regridded variable """ # there is a circular dependency between cdms2 and regrid2. In # principle, cdms2 files should not import regrid2, we're bending # rules here... import regrid2 from regrid2 import Horizontal if togrid is None: return self else: fromgrid = self.getGrid() # this returns the horizontal grid only # default w/o bounds regridTool = 'libcf' regridMethod = 'linear' if 'topology' in self.getAxis(-1).attributes: if self.getAxis(-1).attributes['topology'] == 'circular': # for the ESMF regridders if periodicity is not set. if 'periodicity' not in keywords.keys(): keywords['periodicity'] = guessPeriodicity( getBoundList(_getCoordList(self.getGrid()))) keywords['mkCyclic'] = 1 # for LibCF regridder # check if there are bounds and we have esmf if fromgrid.getBounds() is not None and hasattr(regrid2, "ESMFRegrid"): regridTool = 'esmf' regridMethod = 'linear' # Hum ok if only 1 longitude regrid fails, let's check if len(togrid.getLongitude()) == 1: # esmf can't deal with this regridTool = "regrid2" # let user override userSpecifiesMethod = False for rm in 'rm', 'method', 'regridmethod', 'regrid_method', 'regridMethod': if rm in keywords: regridMethod = keywords[rm] del keywords[rm] userSpecifiesMethod = True userSpecifiesTool = False for rt in 'rt', 'tool', 'regridtool', 'regrid_tool', 'regridTool': if rt in keywords: regridTool = keywords[rt] del keywords[rt] userSpecifiesTool = True # the method determines the tool if re.search('conserv', regridMethod, re.I) or \ re.search('patch', regridMethod, re.I): # only esmf can do conservative and patch regridTool = 'esmf' # make sure the tool can do it if re.search('^regrid', regridTool, re.I) is not None and \ (len(fromgrid.getLatitude().shape) > 1 or len(togrid.getLatitude().shape) > 1): message = """ avariable.regrid: regrid2 cannot do curvilinear, will switch to esmf..." """ warnings.warn(message, Warning) regridTool = 'esmf' if re.search('esmf', regridTool, re.I): # make sure source grids have bounds haveBounds = True for g in fromgrid, : for c in g.getLatitude(), g.getLongitude(): haveBounds &= (c.getBounds() is not None) if not haveBounds: message = """ avariable.regrid: regridTool = 'esmf' requires bounds for source grid, will switch to regridTool = 'libcf' """ warnings.warn(message, Warning) regridTool = 'libcf' regridMethod = 'linear' if not hasattr(regrid2, "ESMFRegrid"): message = """ avariable.regrid: regridTool = 'esmf' but your version does not seems to be built with esmf, will switch to regridTool = 'libcf' """ warnings.warn(message, Warning) regridTool = 'libcf' regridMethod = 'linear' if re.search('conserv', regridMethod, re.I): # make sure destination grid has bounds haveBounds = True for g in togrid, : for c in g.getLatitude(), g.getLongitude(): haveBounds &= (c.getBounds() is not None) if not haveBounds: message = """ avariable.regrid: regridMethod = 'conserve' requires bounds for destination grid, will switch to regridMethod = 'linear' """ warnings.warn(message, Warning) regridMethod = 'linear' if not userSpecifiesTool and re.search( '^regrid', regridTool, re.I) is None: message = """ avariable.regrid: We chose regridTool = %s for you among the following choices: Tools -> 'regrid2' (old behavior) 'esmf' (conserve, patch, linear) or 'libcf' (linear)""" % regridTool warnings.warn(message, Warning) if not userSpecifiesMethod and re.search( '^regrid', regridTool, re.I) is None: message = """ avariable.regrid: We chose regridMethod = %s for you among the following choices: 'conserve' or 'linear' or 'patch'""" % regridMethod warnings.warn(message, Warning) if re.search('^regrid', regridTool, re.I): if 'diag' in keywords and \ isinstance(keywords['diag'], dict): keywords['diag']['regridTool'] = 'regrid' # the original cdms2 regridder regridf = Horizontal(fromgrid, togrid) return regridf(self, missing=missing, order=order, mask=mask, **keywords) # emsf or libcf... srcGridMask = None # set the source mask if a mask is defined with the source data # if numpy.any(self.mask == True): # srcGridMask = getMinHorizontalMask(self) # compute the interpolation weights ro = CdmsRegrid(fromgrid, togrid, dtype=self.dtype, regridMethod=regridMethod, regridTool=regridTool, srcGridMask=srcGridMask, srcGridAreas=None, dstGridMask=None, dstGridAreas=None, **keywords) # now interpolate return ro(self, **keywords)
[docs] def pressureRegrid(self, newLevel, missing=None, order=None, method="log"): """Return the variable regridded to new pressure levels. The variable should be a function of lat, lon, pressure, and (optionally) time. Parameters ---------- newLevel : is an axis of the result pressure levels. method : is optional, either `log` to interpolate in the log of pressure (default), or `linear` for linear interpolation. missing and order : are as for regrid.PressureRegridder. """ from regrid2 import PressureRegridder fromlevel = self.getLevel() if fromlevel is None: raise CDMSError('No pressure level') pregridf = PressureRegridder(fromlevel, newLevel) result = pregridf(self, missing=missing, order=order, method=method) return result
[docs] def crossSectionRegrid(self, newLevel, newLatitude, missing=None, order=None, method="log"): """Return the variable regridded to new pressure levels and latitudes. The variable should be a function of lat, level, and (optionally) time. Parameters ---------- newLevel : is an axis of the result pressure levels. newLatitude : is an axis of latitude values. method : Optional either "log" to interpolate in the log of pressure (default), or "linear" for linear interpolation. missing and order : are as for regrid.CrossSectionRegridder. """ from regrid2 import CrossSectionRegridder fromlevel = self.getLevel() fromlat = self.getLatitude() if fromlevel is None: raise CDMSError('No pressure level') if fromlat is None: raise CDMSError('No latitude level') xregridf = CrossSectionRegridder( fromlat, newLatitude, fromlevel, newLevel) result = xregridf(self, missing=missing, order=order, method=method) return result
[docs] def _process_specs(self, specs, keys): """Process the arguments for a getSlice, getRegion, etc. time, level, latitude, longitude keywords handled here Returns ------- An array of specifications for all dimensions. Any Ellipsis has been eliminated. """ myrank = self.rank() nsupplied = len(specs) # numpy will broadcast if we have a new axis in specs # --------------------------------------------------- if [x for x in specs if numpy.array_equal(x, numpy.newaxis)] == [ numpy.newaxis]: nnewaxis = 1 else: nnewaxis = 0 if [x for x in specs if numpy.array_equal(x, Ellipsis)] == [Ellipsis]: nellipses = 1 else: nellipses = 0 if nsupplied - nellipses - nnewaxis > myrank: raise CDMSError(InvalidRegion + "too many dimensions: %d, for variable %s" % (len(specs), self.id)) speclist = [unspecified] * myrank i = 0 j = 0 while i < nsupplied: if (specs[i] is Ellipsis) or (specs[i] is numpy.newaxis): j = myrank - (nsupplied - (i + 1)) else: speclist[j] = specs[i] j = j + 1 i = i + 1 for k, v in list(keys.items()): if k in ['squeeze', 'raw', 'grid', 'order']: continue i = self.getAxisIndex(k) if i >= 0: if speclist[i] is not unspecified: raise CDMSError( 'Conflict between specifier %s and %s' % (repr( speclist[i]), repr(keys))) speclist[i] = v # Replace remaining unspecified for slice(None, ...) for i in range(myrank): if speclist[i] is unspecified: speclist[i] = slice(None, None, None) return speclist
[docs] def _single_specs(self, specs): """Return a list of dimension indices where the spec is an index.""" myrank = self.rank() nsupplied = len(specs) i = 0 j = 0 singles = [] while i < nsupplied: if specs[i] is Ellipsis: j = myrank - (nsupplied - (i + 1)) else: if isinstance(specs[i], int): singles.append(j) j = j + 1 i = i + 1 return singles
[docs] def specs2slices(self, speclist, force=None): """Create an equivalent list of slices from an index specification. An index specification is a list of acceptable items, which are * an integer * a slice instance (slice(start, stop, stride)) * the object "unspecified" * the object None * a colon The size of the speclist must be self.rank() """ if len(speclist) != self.rank(): raise CDMSError("Incorrect length of speclist in specs2slices.") slicelist = [] for i in range(self.rank()): key = speclist[i] if isinstance(key, int): # x[i] slicelist.append(slice(key, key + 1)) elif isinstance(key, slice): # x[i:j:k] slicelist.append(key) elif key is unspecified or key is None or key == ':': slicelist.append(slice(0, len(self.getAxis(i)))) elif key is Ellipsis: raise CDMSError("Misuse of ellipsis in specification.") elif isinstance(key, tuple): slicelist.append(slice(*key)) else: raise CDMSError('invalid index: %s' % str(key)) # Change default or negative start, stop to positive for i in range(self.rank()): axis = self.getAxis(i) length = len(axis) start = slicelist[i].start stop = slicelist[i].stop step = slicelist[i].step # # allow negative indices in a wrapped (isCircular() = 1) axis # circular = (axis.isCircular() and force is None) altered = 0 if step is None: altered = 1 step = 1 if ((start is None or stop is None or start < 0 or stop < 0) and (circular == 0)): altered = 1 adjustit = 1 if step > 0: if start is None: start = 0 if stop is None: stop = length if start == -1 and stop == 0: stop = length else: if start is None: start = length - 1 if stop is None: # stop=-1 adjustit = 0 if (start is None or start < 0): start = start % length if (stop is None or stop < 0) and adjustit: stop = stop % length if altered: slicelist[i] = slice(start, stop, step) return slicelist
[docs] def reg_specs2slices(self, initspeclist, force=None): # Don't use input to store return value speclist = copy.copy(initspeclist) for i in range(self.rank()): item = speclist[i] if isinstance(item, slice): newitem = item elif item == ':' or item is None or item is unspecified: axis = self.getAxis(i) newitem = slice(0, len(axis)) elif isinstance(item, list) or \ isinstance(item, tuple): axis = self.getAxis(i) if len(item) == 2: # (start,end) indexInterval = axis.mapIntervalExt(item) elif len(item) == 3: # (start,end,'xxx') coordInterval = (item[0], item[1]) indexInterval = axis.mapIntervalExt(coordInterval, item[2]) elif len(item) == 4: coordInterval = (item[0], item[1]) indexInterval = axis.mapIntervalExt( coordInterval, item[2], item[3]) elif len(item) == 5: coordInterval = (item[0], item[1]) indexInterval = axis.mapIntervalExt( coordInterval, item[2], item[3], item[4]) elif len(item) == 6: coordInterval = (item[0], item[1]) indexInterval = axis.mapIntervalExt( coordInterval, item[2], item[3], item[4], item[5]) else: raise CDMSError( InvalidRegion + "invalid format for coordinate interval: %s" % str(item)) if indexInterval is None: raise CDMSError(OutOfRange + str(item)) newitem = slice( indexInterval[0], indexInterval[1], indexInterval[2]) elif isinstance(item, numpy.floating) or \ isinstance(item, float) or \ isinstance(item, numpy.integer) or \ isinstance(item, int) or \ isinstance(item, int) or \ isinstance(item, string_types) or \ type(item) in CdtimeTypes: axis = self.getAxis(i) # # default is 'ccn' in axis.mapIntervalExt # indexInterval = axis.mapIntervalExt((item, item)) if indexInterval is None: raise CDMSError(OutOfRange + str(item)) newitem = slice( indexInterval[0], indexInterval[1], indexInterval[2]) else: raise CDMSError( InvalidRegion + "invalid format for coordinate interval: %s" % str(item)) speclist[i] = newitem slicelist = self.specs2slices(speclist, force) return slicelist
[docs] def _decodedType(self): "The datatype after decoding." if hasattr(self, 'scale_factor') and isinstance( self.scale_factor, numpy.ndarray): result = self.scale_factor.dtype.char elif hasattr(self, 'add_offset') and isinstance(self.add_offset, numpy.ndarray): result = self.add_offset.dtype.char else: result = numpy.float32 return result
[docs] def isEncoded(self): "True if self is represented as packed data." return (hasattr(self, "scale_factor") or hasattr(self, "add_offset"))
[docs] def decode(self, ar): """Decode compressed data. Parameters ---------- ar : is a masked array, scalar, or numpy.ma.masked. """ resulttype = self._decodedType() if hasattr(self, 'scale_factor'): scale_factor = self.scale_factor else: scale_factor = numpy.array([1.0], resulttype) if hasattr(self, 'add_offset'): add_offset = self.add_offset else: add_offset = numpy.array([0.0], resulttype) if ar is not numpy.ma.masked: result = scale_factor * ar.data[:] + add_offset if isinstance(ar, numpy.ma.MaskedArray): result = numpy.ma.masked_array( result, mask=ar.mask).astype(resulttype) numpy.ma.set_fill_value( result, numpy.ma.default_fill_value(0.)) else: tmp = numpy.array(result) result = tmp.astype(resulttype)[0] return result else: return ar
[docs] def getGridIndices(self): """ Get Grid Indices Returns ------- a tuple of indices corresponding to the variable grid.""" grid = self.getGrid() result = [] if grid is not None: gridaxes = grid.getAxisList() varaxes = [d[0] for d in self.getDomain()] for i in range(len(gridaxes)): for j in range(len(varaxes)): if gridaxes[i] is varaxes[j]: result.append(j) break else: raise CDMSError( 'Variable and grid do not share common dimensions: %s' % self.id) return tuple(result)
# numpy.ma overrides
[docs] def squeeze(self): return(MV.squeeze(self))
[docs] def __getitem__(self, key): if isinstance(key, tuple): speclist = self._process_specs(key, {}) else: if isinstance(key, int) and key >= len(self): raise IndexError("Index too large: %d" % key) speclist = self._process_specs([key], {}) if [x for x in speclist if (isinstance( x, numpy.ndarray) or isinstance(x, list))] != []: index = [ x for x in speclist if ( isinstance( x, numpy.ndarray) or isinstance( x, list))] return self.data.take(index) # Note: raw=0 ensures that a TransientVariable is returned return self.getSlice(numericSqueeze=1, raw=0, isitem=1, *speclist)
[docs] def __getslice__(self, low, high): # Note: raw=0 ensures that a TransientVariable is returned return self.getSlice(slice(low, high), numericSqueeze=1, raw=0)
[docs] def typecode(self): raise CDMSError(NotImplemented + 'typecode')
[docs] def __abs__(self): return MV.absolute(self)
[docs] def __neg__(self): return MV.negative(self)
[docs] def __add__(self, other): return MV.add(self, other)
[docs] def __copy__(self): return MV.array(self, copy=1)
__radd__ = __add__
[docs] def __lshift__(self, n): return MV.left_shift(self, n)
[docs] def __rshift__(self, n): return MV.right_shift(self, n)
[docs] def __sub__(self, other): return MV.subtract(self, other)
[docs] def __rsub__(self, other): return MV.subtract(other, self)
[docs] def __mul__(self, other): return MV.multiply(self, other)
__rmul__ = __mul__
[docs] def __floordiv__(self, other): return MV.floor_divide(self, other)
[docs] def __truediv__(self, other): return MV.true_divide(self, other)
[docs] def __div__(self, other): return MV.divide(self, other)
[docs] def __rdiv__(self, other): return MV.divide(other, self)
[docs] def __pow__(self, other, third=None): return MV.power(self, other, third)
[docs] def __iadd__(self, other): "Add other to self in place." return MV.add(self, other) # See if we can improve these later.
[docs] def __isub__(self, other): "Subtract other from self in place." return MV.subtract(self, other) # See if we can improve these later.
[docs] def __imul__(self, other): "Multiply self by other in place." return MV.multiply(self, other) # See if we can improve these later.
[docs] def __idiv__(self, other): "Divide self by other in place." return MV.divide(self, other) # See if we can improve these later.
[docs] def __eq__(self, other): return MV.equal(self, other)
[docs] def __ne__(self, other): return MV.not_equal(self, other)
[docs] def __lt__(self, other): return MV.less(self, other)
[docs] def __le__(self, other): return MV.less_equal(self, other)
[docs] def __gt__(self, other): return MV.greater(self, other)
[docs] def __ge__(self, other): return MV.greater_equal(self, other)
[docs] def __sqrt__(self): return MV.sqrt(self)
[docs] def astype(self, tc): "return self as array of given type." return self.subSlice().astype(tc)
# internattr.add_internal_attribute(AbstractVariable, 'id', 'parent') # PropertiedClasses.set_property(AbstractVariable, 'missing_value', acts=AbstractVariable._setmissing, nodelete=1) __rp = r'\s*([-txyz0-9]{1,1}|\(\s*\w+\s*\)|[.]{3,3})\s*' __crp = re.compile(__rp)
[docs]def orderparse(order): """Parse an order string. Returns a list of axes specifiers. Note: Order elements can be: * Letters t, x, y, z meaning time, longitude, latitude, level * Numbers 0-9 representing position in axes * The letter - meaning insert the next available axis here. * The ellipsis ... meaning fill these positions with any remaining axes. * (name) meaning an axis whose id is name """ if not isinstance(order, string_types): raise CDMSError('order arguments must be strings.') pos = 0 result = [] order = order.strip() while pos < len(order): m = __crp.match(order, pos) if m is None: break r = m.group(1) if r[0] == '(': pass elif r == '...': r = Ellipsis elif len(r) == 1: if r in string.digits: r = int(r) result.append(r) pos = m.end(0) if pos != len(order): raise CDMSError('Order string "' + order + '" malformed, index ' + str(pos)) return result
[docs]def order2index(axes, order): """Find the index permutation of axes to match order. The argument order is a string. Note: Order elements can be: * Letters t, x, y, z meaning time, longitude, latitude, level. * Numbers 0-9 representing position in axes * The letter - meaning insert the next available axis here. * The ellipsis ... meaning fill these positions with any remaining axes. * (name) meaning an axis whose id is name. """ if isinstance(order, string_types): result = orderparse(order) elif isinstance(order, list): result = order else: raise CDMSError( 'order2index, order specified of bad type:' + str(type(order))) n = len(axes) permutation = [None] * n j = 0 pos = 0 while j < len(result): item = result[j] if isinstance(item, string_types): if item == 't': spec = 'time' elif item == 'x': spec = 'longitude' elif item == 'y': spec = 'latitude' elif item == 'z': spec = 'level' elif item == '-': pos += 1 j += 1 continue else: spec = item[1:-1] for k in range(n): if axisMatches(axes[k], spec): if k in permutation: raise CDMSError('Duplicate item in order %s' % order) permutation[pos] = k pos += 1 break else: raise CDMSError('No axis matching order spec %s' % str(item)) elif isinstance(item, int): if item in permutation: raise CDMSError('Duplicate item in order %s' % order) if item >= n: raise CDMSError('Index %d out of range in order %s' % (item, order)) permutation[pos] = item pos += 1 elif item is Ellipsis: nleft = len(result) - j - 1 pos = n - nleft else: raise CDMSError( 'List specified for order contains bad item: ' + repr(item)) j += 1 for i in range(n): if i not in permutation: for j in range(n): if permutation[j] is None: permutation[j] = i break return permutation
from .tvariable import TransientVariable # noqa from . import MV2 as MV # noqa